D-Waveから量子ゲートマシンまでのRBMのボルツマン学習


はじめに

普段量子コンピュータの勉強会をしているのですが、D-Waveを使用した機械学習に関して興味がある方が多いのと、質問が多いのでまず基本的な学習の過程のおさらいとD-Waveを活用してまずは簡単な例題を解いて見たいと思います。


参考資料や記事

これまで様々な問題を検証して記事にしてきました。基本的な用語などは過去の記事から参照してしまっているので、そちらもご覧ください。

イジングモデルや量子アニーリングについては過去の記事

「量子アニーリング、イジングモデルとフレームワーク」

https://qiita.com/YuichiroMinato/items/16574cdf88659dbe9783

また、Yoshua Bengio先生の論文にもある通り、ソフトウェアシミュレーションが活用できます。自社でもGPUやスパコンを活用したシミュレーションを活用していましたので、実機が活用できない環境でもシミュレーションでの実装をお勧めします。


RBMとは

RBMは制限ボルツマンマシン、制限付きボルツマンマシンなどと呼ばれ、ボルツマンマシンと呼ばれるタイプのNNに層を制限して2層の可視層と隠れ層に制限したものです。

参考:

Restricted Boltzmann machine

https://en.wikipedia.org/wiki/Restricted_Boltzmann_machine


D-WaveマシンとRBMとの関係とはじまり

モントリオール大学のこの論文が始まりです。ディープラーニングでの有名な先生がD-Waveでの深層学習の実装可能性として下記の論文を発表し、D-Waveマシンへの実装が検討されました。

On the Challenges of Physical Implementations of RBMs

Vincent Dumoulin and Ian J. Goodfellow and Aaron Courville and Yoshua Bengio

https://arxiv.org/pdf/1312.5258.pdf

この際にはD-Waveを模したソフトウェアシミュレータでNLLとよばれる学習の際に必要なRBMの勾配の推定にサンプリングが活用できるという方針は模索されましたが、実機は使われませんでした。

かなりいい論文で、ボルツマンマシンの基本が書かれています。この論文とGeoffrey Hinton先生の下記の論文はRBMを学ぶ上で大変役立ちました。

A Practical Guide to Training Restricted Boltzmann Machines

Geoffrey Hinton

https://www.cs.toronto.edu/~hinton/absps/guideTR.pdf


量子ボルツマンサンプリングの提唱

次に、実際に上記の問題をD-Wave実機で行なった結果として、

Application of Quantum Annealing to Training of Deep Neural Networks

Steven H. Adachi, Maxwell P. Henderson

https://arxiv.org/abs/1510.06356

があります。これは量子アニーラをボルツマンサンプリングマシンとして、上記のNLLの推定に役立てようという具体的な方針が示され、実機で学習が実際にされています。

12.png

引用:https://arxiv.org/pdf/1510.06356.pdf


D-Waveクラウドシステムへの実装

上記のボルツマンサンプリングを活用した機能は実はD-Waveマシンに搭載がされています。論文中に出て来る逆温度パラメータβとともに、機械学習への応用ということで近年のD-Waveの主要機能としてD-Wave社は提供を開始し、活用を期待しているようです。


ボルツマンサンプリング、RBM、DBMを活用した強化学習

さらに、上記のサンプリング機能を使って、RBMを複数接続したDBMを学習し、簡単な強化学習を解くという論文も出ています。

Reinforcement Learning Using Quantum Boltzmann Machines

Daniel Crawford, Anna Levit, Navid Ghadermarzy, Jaspreet S. Oberoi, Pooya Ronagh

https://arxiv.org/abs/1612.05695

詳細なアルゴリズムも掲載されていますので、実装に関してはとても有用かと思います。

13.png

引用:https://arxiv.org/pdf/1612.05695.pdf


今回やること・学ぶこと

今回は上記2つの論文を確認しながら、自社で借りている実機のD-Waveを使いながら、RBMの基礎とサンプリング学習の実装をみていきたいと思います。また、D-Waveのマニュアルにも最近は正式にボルツマンマシンのドキュメントが搭載されていますので、それを参照しながら機械学習の例題を行なってみたいと思います。

Application of Quantum Annealing to Training of Deep Neural Networks

Steven H. Adachi, Maxwell P. Henderson

https://arxiv.org/abs/1510.06356

Reinforcement Learning Using Quantum Boltzmann Machines

Daniel Crawford, Anna Levit, Navid Ghadermarzy, Jaspreet S. Oberoi, Pooya Ronagh

https://arxiv.org/abs/1612.05695


使用するツールとマシンについて

ツール類は今回はpythonを使います。マシンはD-Waveを使用します。

14.jpeg

引用元:https://www.dwavesys.com/resources/media-resources

15.jpeg

引用元:https://www.dwavesys.com/resources/media-resources

詳細はカナダのホームページを参照してくださいませ。

https://www.dwavesys.com/

日本語サイトもあるようです。

http://dwavejapan.com/


モデルと結合について

実は最初に前提としてモデルを構成する結合数に関しての考察がされているので、それは遵守したいと思います。

On the Challenges of Physical Implementations of RBMs

Vincent Dumoulin and Ian J. Goodfellow and Aaron Courville and Yoshua Bengio

https://arxiv.org/pdf/1312.5258.pdf

この論文中で、D-Waveの結合数とRBMの結合に関して、RBMの結合を簡略化して学習ができないかどうかという試みがされました。実際にはD-Waveはキメラグラフと呼ばれる結合を持っており、1量子ビットからの接続数が少ないので、RBMの結合を簡略化してD-Waveに合わせて学習ができないかということが検討されました。

しかし、上記の論文での結論は結合数を減らすこと(モデルを簡略化すること)では学習がうまくいかないという結論でしたので、以降はRBMの結合数をきちんと守った形で結合を遵守するという方針があります。

16.png

下記にRBMとchimera restricted RBMの模式図を描いてみました。

17.png

8量子ビットある場合、上記の完全二部グラフは結合数が、N*Nありますが、下のキメラグラフは元々の結合が少ないので、二部グラフがつくりにくいという事情があります。

ただ、下のChimera Restricted RBMはうまくはいかないので、上記のRBMモデルを忠実に再現するということを、D-Waveのキメラグラフを使いながら実現していく必要があります。

早速仕組みを少しずつみていきたいと思います。


モデル概略

まず、モデルは可視層と呼ばれる層と隠れ層と呼ばれる2層からなるモデルです。結合に方向性がなく無向グラフです。

18.png

引用:https://arxiv.org/pdf/1510.06356.pdf

結合分布がギブス分布、ボルツマン分布に従うようです。

p(v,h) = \frac{1}{Z}exp(-E(v,h))


結合分布

確率論において、同時分布(どうじぶんぷ)または結合分布(けつごうぶんぷ, joint distribution)とは、確率変数が複数個ある場合に、複数の確率変数がとる値の組に対して、その発生の度合いを確率を用いて記述するもので、確率分布の一種である。

引用:https://ja.wikipedia.org/wiki/%E5%90%8C%E6%99%82%E5%88%86%E5%B8%83


ボルツマン分布

ボルツマン分布(ボルツマンぶんぷ、英語: Boltzmann distribution)は、一つのエネルギー準位にある粒子の数(占有数)の分布を与える理論式の一つである。ギブス分布とも呼ばれる。気体分子の速度の分布を与えるマクスウェル分布をより一般化したものに相当する。

量子統計力学においては、占有数の分布がフェルミ分布に従うフェルミ粒子と、ボース分布に従うボース粒子の二種類の粒子に大別できる。ボルツマン分布はこの二種類の粒子の違いが現れないような条件におけるフェルミ分布とボーズ分布の近似形(古典近似)である。ボルツマン分布に従う粒子は古典的粒子とも呼ばれる。

引用:https://ja.wikipedia.org/wiki/%E5%90%8C%E6%99%82%E5%88%86%E5%B8%83


RBMの結合分布

ここで、この結合分布はエネルギー関数で規定され、可視層のノード数nと隠れ層のノード数mとで下記のように示されます。

E(v,h) = -\sum_{i=1}^nb_iv_i-\sum_{j=1}^mc_jh_j-\sum_{i=1}^n\sum_{j=1}^mW_{ij}v_ih_j\\ 

ここで、v_i,h_j\in\{0,1\}

また、正規化のための規格化定数、分配関数は下記の通りになります。

Z = \sum_{v_k}\sum_{h_l}exp\left(\sum_kb_kv_k + \sum_lc_lh_l + \sum_{k,l}W_{kl}v_kh_l  \right)

また、完全二部グラフより、条件付き確率分布はそれぞれ、vとhについてシグモイド関数を用いて下記の通りになります。

P(h_j=1|v) = sigm(c_j+\sum_iW_{ij}v_i)\\ 

P(v_i=1|h) = sigm(b_i+\sum_jW_{ij}h_j)


学習について

次に上記の確率分布からなるNNのモデルの学習方法を確認したいと思います。複層のDBN(DBM)でも、RBMの形にして学習をします。これらの学習はlogPを最大にするように、トレーニングデータとモデルの誤差計算をします。結合係数やバイアスの勾配計算はlogPを用いて下記のようになります。

\frac{\partial logP}{\partial W_{ij}} = < v_ih_j > _{data}- < v_ih_j > _{model}\\ 

\frac{\partial logP}{\partial b_i} = < v_i > _{data}- < v_i > _{model}\\ \frac{\partial logP}{\partial c_j} = < h_j > _{data}- < h_j > _{model}

ここで、上記の結合係数の勾配計算で、下記モデルの期待値の計算はあまり効率的な計算がありません。

< v_ih_j > _{model} = \frac{1}{Z}\sum_{\{v_k\}}\sum_{\{h_l\}}v_ih_j exp\left(\sum_kb_kv_k + \sum_lc_lh_l+\sum_{kl}W_{kl}v_kh_l\right)

実際には、この値を直接求めるのが大変なので、近似計算でCD法など、Gibbsサンプリングを応用した手法で、順番に隠れ層と可視層の計算を行なって値を取るという方法が行われます。

参考記事:https://qiita.com/t_Signull/items/f776aecb4909b7c5c116

なので、CD法などでは計算量を節約する意味で、かなり近似的な計算が行われ、これを精度改善しようとすると計算量と時間がかかります。

\Delta W_{ij} = \epsilon[ < v_ih_j > _{data}- < v_ih_j > _{recon}]

上記は結合係数の更新ですが、バイアスの更新でも同様です。


量子ボルツマンサンプリングを活用したパラメータ更新

ここで、このCD法の代わりに量子アニーリングマシンを使って元の勾配計算どおりしまおうというのが実機を使ったボルツマンサンプリング学習の基本です。

12.png

引用:https://arxiv.org/pdf/1510.06356.pdf

ボルツマンサンプリングを使った際の結合係数やバイアスの更新式は、下記の通りです。αはモーメンタムでϵは学習率です。

W_{ij}^{(t+1)} = \alpha W_{ij}^{(t)} + \epsilon[ < v_i h_j > _{data} – < v_ih_j > _{model}]\\ 

b_{i}^{(t+1)} = \alpha b_i^{(t)} + \epsilon[ < v_i > _{data} – < v_i > _{model}]\\
c_{j}^{(t+1)} = \alpha c_j^{(t)} + \epsilon[ < h_j > _{data} – < h_j > _{model}]

RBMを学習した後は古典計算機でバックプロパゲーションで仕上げの学習をするようです。


量子アニーラを活用したサンプリング手法

D-Waveマシンは量子アニーリングという理論をベースに作られており、基本的には最適化問題を解くために最小値問題を求めるような設計になっています。しかし実際には、外部環境の影響やその他の理由で、最適解に落ち着かない場合が多いです。そのため、この局所解へ落ちる性質を利用して、分布を求めるサンプリングマシンとして活用するというのが、サンプリング学習の基本です。

この励起状態のばらつきがボルツマン分布として仮定できると、下記式に近似できこれと初期に出てきた分布の式を対応させることで、サンプリング手法を更新式に導入できそうです。

P(x) = \frac{1}{Z}exp\left(-\beta_{eff}H_f(x)\right)

Hfは最終的に求めるハミルトニアン(エネルギーコスト関数)で、βeffはサンプリングの分布を調整する変数です。

これを使用することで、一番計算量のかかる部分を下記のように近似します。

\overline{v_ih_j} = \frac{1}{N}\sum_{n=1}^N v_i^{(n)}h_j^{(n)}

これをmodelに適用します。他の可視層や隠れ層のモデル期待値も同様です。


D-Waveマシンを想定した結合方法

D-Waveはキメラグラフと呼ばれる結合方式を採用しているので、RBMの完全二部グラフを再現するには少し工夫が必要です。

具体的には、量子ビット同士を繋げてコピー操作ができるので、それを活用してRBMを実現します。繋ぎ方は下記のようです。実際には完全二部グラフはとてもわかりやすい、縦横の交差になります。この交差部分に結合係数の値を設定していきます。

18.png

上の図では左側はハードウェアを想定した結合、実際にハードは量子ビットの形状は細長くなっています。右側は論理的に書いた図です。


QUBOの作成と、逆温度パラメータ

D-Waveで特には、QUBOmatrixというものを作ります。QUBOmatrixは問題をバイナリ値でかんがえた時に作る行列で、これを少し変形してD-Waveにかけます。

77.png

引用:https://arxiv.org/pdf/1510.06356.pdf

結果として、対角項にバイアスbとcがはいり、非対角項に結合係数wが入ります。バイアスは局所磁場とか磁場項とよび、結合係数は相互作用項と呼んだりします。


QUBOからイジングへ

次に、逆温度のパラメータはここで見ると、行列の外にだせるので全体のスケールを調整するパラメータになってます。QUBOからイジングに直す時にはまだかんがえなくていいので、イジングmatrixに直します。


あとはスケーリングパラメータを決めてD-Waveでサンプリングを実行

あとはβeffを決めて実行です。どうやら論文によると、βeffのあたいは実行するサンプルによって依存するようです。MNISTの画像認識の学習では、画像サイズが大きくなるにつれて、βeffは小さくするのが適切のようでした。

また、D-Waveの量子ビット数は多くないので、完全二部グラフを作って学習をさせる際には、MNISTをダウンサンプリングして実装しています。

MNISTの28*28pixelの画像を6*6まで下げています。

20.png

引用:https://arxiv.org/pdf/1510.06356.pdf

実行結果などの学習具合は詳しくは論文を参考にしてもらって、だいたい下記のように収束が良いようです。

22.png

引用:https://arxiv.org/pdf/1510.06356.pdf


続いて実機を使ってこれを実装する

D-Waveには前処理を変更して、逆温度パラメータを設定する項目があるので、これを使用してみます。

111.png


簡単な例題で実行する

簡単な例題で上記のRBMのサンプリングに関して学んでみたいと思います。

まずはD-Waveの例題でも載っている、簡単な問題で。

Q = \left[ 

\begin{array}{rr}
-1 & 2 \\
& -1 \\
\end{array}
\right]

このQUBOmatrixの場合、エネルギーコスト関数は、

E(x) = -x_1-x_2+2x_1x_2

となります。xの場合の数でエネルギーを求めて見ると、

x1.png.pagespeed.ic.Uon8jzb8v3.png

解析的にかんがえて見ると、ここではすべての場合の数が出せているので、まず規格化定数Zを求めます。

Z = \sum exp(E(x)) = exp(0) + exp(1) + exp(1) + exp(0) = 1+2.718+2.718+1 = 7.44

そして、

P = \frac{1}{Z}exp(E)

より、確率は0.13と0.37がでます。

11.png

これを手作業で出すのは大変なので、サンプリングでかんがえてみます。


D-Waveでサンプリングする

まずこれをQUBOからイジングに直します。

Q = \left[ 

\begin{array}{rr}
0 & 0.5 \\
& 0 \\
\end{array}
\right]

こうなりましたので、これをいれます。

係数のJijに-0.5いれて、1000回のサンプリングでかけると、、、

111-1.png

結構衝撃的です。。。大体の確率でサンプリングできてる。。。まじですか、、、

1111.png

解を確認したところ大丈夫そうです。

12.png

速度も速い。

ということで、サンプリングを行なった結果、モデルの分布が取れています。先ほどは解析的にときましたが、サンプリングから逆算してモデルの期待値が取れますので、高速に誤差計算できそうです。サンプリングした結果は古典計算機で値の更新をします。

大成功でした。


networkx

描画について少しみたいと思います。

RBMを量子アニーリングでやる際には無向グラフで十分ですので、networkxとmatplotlibを使ってやってみたいと思います。

import matplotlib.pyplot as plt

import networkx as nx

#無向グラフ
G = nx.Graph()

#RBMで使うエッジを指定、ノードは自動的に指定される
G.add_edges_from([(0,4),(0,5),(0,6),(0,7),(1,4),(1,5),(1,6),(1,7),(2,4),(2,5),(2,6),(2,7),(3,4),(3,5),(3,6),(3,7)])

#綺麗に並べたいので場所指定
pos = {0:[0.1,0.4],1:[0.1,0.3],2:[0.1,0.2],3:[0.1,0.1],4:[0.2,0.4],5:[0.2,0.3],6:[0.2,0.2],7:[0.2,0.1]}

nx.draw_networkx(G,pos)

ぼちぼちいい感じではないでしょうか。

13.png

次に、量子ビットを指定し、バイアスと結合荷重を初期化してみます。使用するのは量子ビットを表すリストのqとバイアスと結合荷重が上三角行列のmatrixになったJです。

class rbm:

def __init__(self):
self.J = []
self.v = 0
self.h = 0
self.n = 0

def setRBM(self,v=4,h=4):
n = v+h
self.v = v
self.h = h
self.n = n
self.J = np.diag([np.random.randint(99)/100 for i in range(n)])
for i in range(v):
for j in range(v,n):
self.J[i][j] = np.random.randint(99)/100
return self

def network(self):
G= nx.Graph()
edges = []
for i in range(self.v):
edges += [(i,j) for j in range(self.v,self.n)]
G.add_edges_from(edges)
pos = {}
for i in range(self.n):
if i<self.v:
pos[i] = [0.1,0.1*(self.v-i)]
else:
pos[i] = [0.2,0.1*(self.h+self.v-i)]
nx.draw_networkx(G,pos)
plt.show()

ついでにネットワークまで設定してくれるようにしました。こんな感じです。

import matplotlib.pyplot as plt

import networkx as nx
import numpy as np

a = rbm().setRBM(3,6).network()

14.png

まぁ、こんなもんでしょう。


アニーリング計算

この状態でアニーリングをかけてみます。アニーリングは任意の状態からスタートして、バイアスと結合荷重で結果が決まります。また、RBMの場合には近い値での局所解がたくさんありそうです。

from blueqat.opt import Opt

a = rbm().setRBM(4,4)
a.J

#=>
array([[0. , 0. , 0. , 0.52, 0.87, 0.09],
[0. , 0.08, 0. , 0.21, 0.35, 0.22],
[0. , 0. , 0.31, 0.51, 0.62, 0.3 ],
[0. , 0. , 0. , 0.7 , 0. , 0. ],
[0. , 0. , 0. , 0. , 0.55, 0. ],
[0. , 0. , 0. , 0. , 0. , 0.81]])

#これをアニーリングにかけてみます。
#複数回を使うカウンター
import collections
def counter(narr):
dis = []
for i in range(len(narr)):
dis.append(''.join([str(x) for x in narr[i]]))
return collections.Counter(dis)

c = Opt().add(a.J).run(shots=100)

#カウント
counter(c)

#=>
Counter({'000000': 48, '100000': 49, '010000': 2, '110000': 1})


全体の分布

上記は4つの状態が出てきました。この分布を使って更新を行います。分布は温度によって影響されます。最近ではpause+quenchなどの機能もでていますが、今回は逆温度パラメータと呼ばれるものを使って分布を調整してみます。

def bars(cc):

index = []
value = []
for k,v in counter(cc).items():
index.append(k)
value.append(v)
plt.bar(index,value)
plt.show()

bars(c)

これで、

Figure_1-2.png

こんな感じです。逆温度パラメータはjijの前にパラメータβをつけます。これによって、分布を変えてみます。β=0.1としてみて、

beta = 0.1

c = Opt().add(beta*a.J).run(shots=100)
bars(c)

これにより、

2-2.png

より複雑な分布をとることができました。

Counter({'011000': 4,

'000000': 12,
'001000': 7,
'100000': 29,
'010000': 15,
'110000': 24,
'000001': 1,
'101000': 3,
'111000': 5})


学習について

学習に関しては、学習すべきパラメータが3種類あります。

1、visible layerのbias

2、hidden layerのbias

3、v*hの結合荷重

1は簡単でデータの入力値を使います。

2は一旦visible layerからsigmoidを通してhiddenlayerの値を計算します。

3がめんどいのですが、CD法を使う代わりにアニーリングを使って、v*hの期待値をだしてからの1/Nして平均を取ります。


学習式

学習式はバイアスと結合荷重で、

b_{i,t+1} = \alpha * b_{i,t} + \epsilon(v_{data} -v_{model})\\ c_{i,t+1} = \alpha * c_{i,t} + \epsilon(h_{data} -h_{model})\\ w_{ij,t+1} = \alpha * w_{ij,t} + \epsilon(vh_{data} -vh_{model})

αはモーメンタムでϵは学習率です。基本的にはvはそのまま計算し、hのdataはsigmoidを通してもとめ、wはNLLの計算というもので計算量が多いものです。wのdataはvからhを求めれば求まります。

また、最初にシグモイド関数を用意しておきます。

def sigm(xx):

return 1/(1+np.exp(-xx))


例題

1100というvisible layerを学習させてみます。まずは100回普通に今の回路を計算してみます。

a = rbm().setRBM(4,4)

c = Opt().add(a.J).run(shots=100)
counter(c)

#=>
Counter({'00000000': 100})

#=>
array([[0.22, 0. , 0. , 0. , 0.28, 0.58, 0.56, 0.02],
[0. , 0.33, 0. , 0. , 0.81, 0.87, 0.15, 0.86],
[0. , 0. , 0.21, 0. , 0.38, 0.03, 0.64, 0.33],
[0. , 0. , 0. , 0.15, 0.18, 0.7 , 0.71, 0.7 ],
[0. , 0. , 0. , 0. , 0.77, 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. , 0.51, 0. , 0. ],
[0. , 0. , 0. , 0. , 0. , 0. , 0.78, 0. ],
[0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.87]])

こうなりました。visible layerは最初の4量子ビットなので期待値は、

0000/100 = 0000

となります。hidden layerも同様です。visible layerのデータは1100なので、更新式は、それぞれの量子ビットに対して、

b_{0,t+1} = 0.9*0.22+0.1*(1-0)\\b_{1,t+1} = 0.9*0.33+0.1*(1-0)\\b_{2,t+1} = 0.9*0.21+0.1*(0-0)\\b_{3,t+1} = 0.9*0.15+0.1*(0-0)

一方hidden layerの更新をするためには事前にhidden layerのデータ値を出しておく必要があります。例えば、

h_4 = sigm(v_0*w_{04}+v_1*w_{14}+v_2*w_{24}+v_3*w_{34}+c_4)

その上で、

c_{4,t+1} = 0.9*0.77+0.1*(h_4-0)

このようにして更新を。あとは、wに関しても同様に、データからの期待値と上記のv∗hをしたものを使って地味にひとつずつ更新をしていきます。使いやすいように更新用のmatrixを作って、

J = \alpha*J + \epsilon*(M_{data}-M_{model})

と変数を導入して計算できるように、QUBOmatrixの形にしておきます。それぞれのMにはデータの期待値とモデルの期待値をQUBOmatrixの形にしておきます。getMdata(self,vdata)としてvdataをlistでいれたら計算できるように追加しました。

from blueqat.opt import Opt

import networkx as nx
import numpy as np
import collections

class rbm:
中略、、、

def getMdata(self,vdata):
hdata = [0 for i in range(self.h)]
for i in range(self.h):
hdata[i] += self.J[self.v+i][self.v+i]
for j in range(self.v):
hdata[i] += vdata[j]*self.J[j][self.v+i]
hdata[i] = sigm(hdata[i])
Mdata = np.diag(vdata+hdata)
for i in range(self.v):
for j in range(self.h):
Mdata[i][self.v+j] = vdata[i]*hdata[j]
return Mdata

あとはmodelの方のMatrixができればいいですが、アニーリングが入るためちょっとめんどいですね。

    def getMmodel(self):

c = Opt().add(self.J).run(shots=100)
ccc = counter(c)
print(ccc)
Mmodel = np.diag([0 for i in range(self.n)])
for i,j in ccc.items():
temp = [*i]
temp = list(map(lambda x:int(x)/j,temp))
Mmodel = Mmodel + np.diag(temp)
for k in range(self.v):
for l in range(self.h):
Mmodel[k][self.v+l] = temp[k]*temp[self.v+l]
return Mmodel

間違ってたら指摘してください。。。


早速使ってみる

更新はself.Jをどんどん更新します。

    def train(self,vdata,alpha=0.9,epsilon=0.1):

a.J = alpha*a.J-epsilon*(a.getMdata(vdata)-a.getMmodel())
return self

こちらを連続でやってみて値が固定されるかみてみます。


コード全体を確認

import matplotlib.pyplot as plt

from blueqat.opt import Opt
import networkx as nx
import numpy as np
import collections

class rbm:
def __init__(self):
self.J = []
self.v = 0
self.h = 0
self.n = 0

def setRBM(self,v=4,h=4):
n = v+h
self.v = v
self.h = h
self.n = n
self.J = np.diag([np.random.randint(99)/100 for i in range(n)])
for i in range(v):
for j in range(v,n):
self.J[i][j] = np.random.randint(99)/100
return self

def network(self):
G= nx.Graph()
edges = []
for i in range(self.v):
edges += [(i,j) for j in range(self.v,self.n)]
G.add_edges_from(edges)
pos = {}
for i in range(self.n):
if i<self.v:
pos[i] = [0.1,0.1*(self.v-i)]
else:
pos[i] = [0.2,0.1*(self.h+self.v-i)]
nx.draw_networkx(G,pos)
plt.show()

def getMdata(self,vdata):
hdata = [0 for i in range(self.h)]
for i in range(self.h):
hdata[i] += self.J[self.v+i][self.v+i]
for j in range(self.v):
hdata[i] += vdata[j]*self.J[j][self.v+i]
hdata[i] = sigm(hdata[i])
Mdata = np.diag(vdata+hdata)
for i in range(self.v):
for j in range(self.h):
Mdata[i][self.v+j] = vdata[i]*hdata[j]
return Mdata

def getMmodel(self):
c = Opt().add(self.J).run(shots=100)
ccc = counter(c)
print(ccc)
Mmodel = np.diag([0 for i in range(self.n)])
for i,j in ccc.items():
temp = [*i]
temp = list(map(lambda x:int(x)/j,temp))
Mmodel = Mmodel + np.diag(temp)
for k in range(self.v):
for l in range(self.h):
Mmodel[k][self.v+l] = temp[k]*temp[self.v+l]
return Mmodel

def train(self,vdata,alpha=0.9,epsilon=0.1):
a.J = alpha*a.J-epsilon*(a.getMdata(vdata)-a.getMmodel())
return self

def counter(narr):
dis = []
for i in range(len(narr)):
dis.append(''.join([str(x) for x in narr[i]]))
return collections.Counter(dis)

def bars(cc):
index = []
value = []
for k,v in counter(cc).items():
index.append(k)
value.append(v)
plt.bar(index,value)
plt.show()

#以下試してみます。
a = rbm().setRBM(4,4)
a.network()

a.J

#=>
array([[0.47, 0. , 0. , 0. , 0.25, 0. , 0.91, 0.23],
[0. , 0.1 , 0. , 0. , 0.21, 0.75, 0.09, 0.65],
[0. , 0. , 0.45, 0. , 0.43, 0.28, 0.11, 0.64],
[0. , 0. , 0. , 0.44, 0.57, 0.27, 0.06, 0.57],
[0. , 0. , 0. , 0. , 0.62, 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. , 0.45, 0. , 0. ],
[0. , 0. , 0. , 0. , 0. , 0. , 0.73, 0. ],
[0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.37]])

for i in range(10):
a.train([1,1,0,1])

#=>
Counter({'00000000': 99, '01000000': 1})
Counter({'00000000': 97, '01000000': 3})
Counter({'00000000': 67, '01000000': 33})
Counter({'01000000': 79, '01010000': 8, '01010010': 4, '10000100': 3, '11000000': 2, '00000000': 2, '10000000': 1, '10000101': 1})
Counter({'01010010': 100})
Counter({'01010010': 59, '01011010': 40, '01010011': 1})
Counter({'11011100': 98, '11011110': 2})
Counter({'11011111': 60, '11011101': 30, '11011110': 7, '11011100': 3})
Counter({'11011111': 100})
Counter({'11011111': 100})

この状態で、更新されたa.Jに局所磁場でq0=-10,q1=-10をかけてみて入力をして推論してみます。

c = Opt().add(a.J).add("-10*q0-10*q1").run(shots=100)

counter(c)

#=>
Counter({'11011111': 100})

再現できたような気がしますがどうでしょうか?


QAOAで解分布を取ってみる

上記を量子アニーリングからQAOAにするには、

c = Opt().add(self.J).run(shots=100)

の部分をQAOAにするだけです。試しに単体でやってみます。適当な例題を作ってみてやってみます。まずは普通のアニーリング

a = rbm().setRBM(4,4)

モデルを作ります。逆温度パラメータを設定して解分布を取ります。

c = Opt().add(0.1*a.J).run(shots=100)

Figure_1-6.png

次にQAOAをやってみます。先ほどのQUBOを量子ゲートのVQE+QAOAで計算します。

c = Opt().add(0.1*a.J).qaoa(shots=100)

結果は、

c.most_common(12)

#=>
(((0, 0, 0, 0, 0, 0, 0, 0), 0.4072213424296961),
((0, 0, 1, 0, 0, 0, 0, 0), 0.18526818450665403),
((0, 0, 0, 0, 0, 1, 0, 0), 0.11036387100498915),
((0, 0, 0, 0, 0, 0, 1, 0), 0.0563814515181246),
((0, 0, 1, 0, 0, 1, 0, 0), 0.037928302595366385),
((0, 1, 0, 0, 0, 0, 0, 0), 0.03717892125780595),
((0, 0, 0, 0, 0, 0, 0, 1), 0.01960392923651757),
((0, 0, 0, 0, 0, 1, 1, 0), 0.01683173270893821),
((0, 1, 1, 0, 0, 0, 0, 0), 0.016153879541206548),
((0, 0, 0, 1, 0, 0, 0, 0), 0.015227462186776007),
((0, 0, 1, 0, 0, 0, 1, 0), 0.014427956796671807),
((0, 0, 0, 0, 1, 0, 0, 0), 0.007943667895491113))

グラフにすると、

index = []

value = []
k=0
for i,j in c.most_common(12):
index.append(k)
value.append(j)
k += 1
plt.bar(index,value)
plt.show()

2-3.png

並びが変わりましたが、だいたい同じ結果になりました。

このようにして、量子ゲートモデルでも、確率分布をとって計算をすることができました。


QAOAを導入した量子ゲートRBMモデル

QAOAを用いた量子ゲートのRBMモデルの式は下記の部分を変更する必要があります。

    def getMmodel(self,beta):

c = Opt().add(beta*self.J).run(shots=100)
ccc = counter(c)
print(ccc)
Mmodel = np.diag([0 for i in range(self.n)])
for i,j in ccc.items():
temp = [*i]
temp = list(map(lambda x:int(x)/j,temp))
Mmodel = Mmodel + np.diag(temp)
for k in range(self.v):
for l in range(self.h):
Mmodel[k][self.v+l] = temp[k]*temp[self.v+l]
return Mmodel

上記の最適化部分を、

c = Opt().add(0.1*a.J).qaoa(shots=100)

に書き換えます。また、ループ部分に量子ビットの答えの格納を変更処理する必要があります。


RBMの推論

visiblelayerを入力と出力に対応させて、学習と推論を行います。推論の入力値の固定は局所磁場を利用して行います。

c = Opt().add(a.J).add("-10*q0-10*q1").qaoa(shots=100)

c.most_common(12)

#=>
(((0, 0, 0, 0, 0, 0, 0, 0), 0.4020990180474289),
((0, 0, 1, 0, 0, 0, 0, 0), 0.18650948659338823),
((0, 0, 0, 0, 0, 1, 0, 0), 0.10986637438485558),
((0, 0, 0, 0, 0, 0, 1, 0), 0.055083968265518916),
((0, 0, 1, 0, 0, 1, 0, 0), 0.03859321782313931),
((0, 1, 0, 0, 0, 0, 0, 0), 0.03773956476911518),
((0, 0, 0, 0, 0, 0, 0, 1), 0.021609758825796577),
((0, 0, 0, 0, 0, 1, 1, 0), 0.015584625036213188),
((0, 1, 1, 0, 0, 0, 0, 0), 0.01550494310602666),
((0, 0, 0, 1, 0, 0, 0, 0), 0.015501708822164914),
((0, 0, 1, 0, 0, 0, 1, 0), 0.014592269080913394),
((0, 0, 1, 0, 0, 0, 0, 1), 0.007134026224727637))