はじめに
量子コンピュータを利用した機械学習のうちで、量子断熱計算を利用した量子断熱機械学習を見てみたいと思います。
量子断熱時間発展計算
量子断熱計算は、状態ベクトルを断続的に変化させることで基底状態をキープしたまま時間発展を行うことができる理論です。
初期状態のハミルトニアンを
としたときに
時間発展計算は、
となります。課題は基底状態と第一励起状態が最接近する部分ですが、
ボルツマン分布
上記時間発展において、最終的に出てくる分布がボルツマン分布に従うと仮定して、ある状態の出現する確率は、エネルギーに対応して、
Zは規格化定数。となっています。量子コンピュータにおけるイジングモデルのエネルギーと上記のエネルギーを対応させて、機械学習を行います。
Yoshua Bengio先生のいうとおり
また、Yoshua Bengio先生の論文にもある通り、ソフトウェアシミュレーションが活用できます。実機やシミュレーションを通じて量子ゲートや量子アニーリングマシンを活用するための量子断熱機械学習を見てみたいと思います。
RBMとは
RBMは制限ボルツマンマシン、制限付きボルツマンマシンなどと呼ばれ、ボルツマンマシンと呼ばれるタイプのNNに層を制限して2層の可視層と隠れ層に制限したものです。
参考:
Restricted Boltzmann machine
量子コンピュータとボルツマンマシンの本格的な検証は、モントリオール大学のこの論文が始まりです。ディープラーニングでの有名な先生がD-Waveでの深層学習の実装可能性として下記の論文を発表し、D-Waveマシンへの実装が検討されました。
On the Challenges of Physical Implementations of RBMs
Vincent Dumoulin and Ian J. Goodfellow and Aaron Courville and Yoshua Bengio
この際にはD-Waveを模したソフトウェアシミュレータでNLLとよばれる学習の際に必要なRBMの勾配の推定にサンプリングが活用できるという方針は模索されましたが、実機は使われませんでした。
かなりいい論文で、ボルツマンマシンの基本が書かれています。この論文とGeoffrey Hinton先生の下記の論文はRBMを学ぶ上で大変役立ちました。
A Practical Guide to Training Restricted Boltzmann Machines
Geoffrey Hinton
量子ボルツマンサンプリングの提唱
次に、実際に上記の問題をD-Wave実機で行なった結果として、
Application of Quantum Annealing to Training of Deep Neural Networks
Steven H. Adachi, Maxwell P. Henderson
があります。これは量子アニーラをボルツマンサンプリングマシンとして、上記のNLLの推定に役立てようという具体的な方針が示され、実機で学習が実際にされています。
ボルツマンサンプリング、RBM、DBMを活用した強化学習
さらに、上記のサンプリング機能を使って、RBMを複数接続したDBMを学習し、簡単な強化学習を解くという論文も出ています。
Reinforcement Learning Using Quantum Boltzmann Machines
Daniel Crawford, Anna Levit, Navid Ghadermarzy, Jaspreet S. Oberoi, Pooya Ronagh
詳細なアルゴリズムも掲載されていますので、実装に関してはとても有用かと思います。
モデルと結合について
実は最初に前提としてモデルを構成する結合数に関しての考察がされているので、それは遵守したいと思います。
On the Challenges of Physical Implementations of RBMs
Vincent Dumoulin and Ian J. Goodfellow and Aaron Courville and Yoshua Bengio
この論文中で、D-Waveの結合数とRBMの結合に関して、RBMの結合を簡略化して学習ができないかどうかという試みがされました。実際にはD-Waveはキメラグラフと呼ばれる結合を持っており、1量子ビットからの接続数が少ないので、RBMの結合を簡略化してD-Waveに合わせて学習ができないかということが検討されました。
しかし、上記の論文での結論は結合数を減らすこと(モデルを簡略化すること)では学習がうまくいかないという結論でしたので、以降はRBMの結合数をきちんと守った形で結合を遵守するという方針があります。
下記にRBMとchimera restricted RBMの模式図を描いてみました。
8量子ビットある場合、上記の完全二部グラフは結合数が、N*Nありますが、下のキメラグラフは元々の結合が少ないので、二部グラフがつくりにくいという事情があります。
ただ、下のChimera Restricted RBMはうまくはいかないので、上記のRBMモデルを忠実に再現するということを、D-Waveのキメラグラフを使いながら実現していく必要があります。
早速仕組みを少しずつみていきたいと思います。
モデル概略
まず、モデルは可視層と呼ばれる層と隠れ層と呼ばれる2層からなるモデルです。結合に方向性がなく無向グラフです。
引用:https://arxiv.org/pdf/1510.06356.pdf
結合分布がギブス分布、ボルツマン分布に従います。
結合分布
確率論において、同時分布(どうじぶんぷ)または結合分布(けつごうぶんぷ, joint distribution)とは、確率変数が複数個ある場合に、複数の確率変数がとる値の組に対して、その発生の度合いを確率を用いて記述するもので、確率分布の一種である。
引用:https://ja.wikipedia.org/wiki/同時分布
ボルツマン分布
ボルツマン分布(ボルツマンぶんぷ、英語: Boltzmann distribution)は、一つのエネルギー準位にある粒子の数(占有数)の分布を与える理論式の一つである。ギブス分布とも呼ばれる。気体分子の速度の分布を与えるマクスウェル分布をより一般化したものに相当する。
量子統計力学においては、占有数の分布がフェルミ分布に従うフェルミ粒子と、ボース分布に従うボース粒子の二種類の粒子に大別できる。ボルツマン分布はこの二種類の粒子の違いが現れないような条件におけるフェルミ分布とボーズ分布の近似形(古典近似)である。ボルツマン分布に従う粒子は古典的粒子とも呼ばれる。
引用:https://ja.wikipedia.org/wiki/同時分布
RBMの結合分布
ここで、この結合分布はエネルギー関数で規定され、可視層のノード数nと隠れ層のノード数mとで下記のように示されます。
また、正規化のための規格化定数、分配関数は下記の通りになります。
また、完全二部グラフより、条件付き確率分布はそれぞれ、vとhについてシグモイド関数を用いて下記の通りになります。
学習について
次に上記の確率分布からなるNNのモデルの学習方法を確認したいと思います。複層のDBN(DBM)でも、RBMの形にして学習をします。これらの学習はlogPを最大にするように、トレーニングデータとモデルの誤差計算をします。結合係数やバイアスの勾配計算はlogPを用いて下記のようになります。
ここで、上記の結合係数の勾配計算で、下記モデルの期待値の計算はあまり効率的な計算がありません。
実際には、この値を直接求めるのが大変なので、近似計算でCD法など、Gibbsサンプリングを応用した手法で、順番に隠れ層と可視層の計算を行なって値を取るという方法が行われます。
参考記事:https://qiita.com/t_Signull/items/f776aecb4909b7c5c116
なので、CD法などでは計算量を節約する意味で、かなり近似的な計算が行われ、これを精度改善しようとすると計算量と時間がかかります。
上記は結合係数の更新ですが、バイアスの更新でも同様です。
量子ボルツマンサンプリングを活用したパラメータ更新
ここで、このCD法の代わりに量子断熱計算を使おうというのが実機を使ったボルツマンサンプリング学習の基本です。ボルツマンサンプリングを使った際の結合係数やバイアスの更新式は、下記の通りです。αはモーメンタムでϵは学習率です。
RBMを学習した後は古典計算機でバックプロパゲーションで仕上げの学習をするようです。
量子断熱計算を活用したサンプリング手法
量子断熱過程を計算し、分布を求めるサンプリングマシンとして活用するというのが、サンプリング学習の基本です。この励起状態のばらつきがボルツマン分布として仮定できると、下記式に近似できこれと初期に出てきた分布の式を対応させることで、サンプリング手法を更新式に導入できそうです。
Hfは最終的に求めるハミルトニアン(エネルギーコスト関数)で、βeffはサンプリングの分布を調整する変数です。
これを使用することで、一番計算量のかかる部分を下記のように近似します。
これを<vihj>modelに適用します。他の可視層や隠れ層のモデル期待値も同様です。
QUBOの作成と、逆温度パラメータ
イジングモデルを01のバイナリ値に拡張した、QUBOmatrixというものを作ります。QUBOmatrixは問題をバイナリ値でかんがえた時に作る行列で、これを少し変形して実機にかけます。
結果として、対角項にバイアスbとcがはいり、非対角項に結合係数wが入ります。バイアスは局所磁場とか磁場項とよび、結合係数は相互作用項と呼んだりします。
QUBOからイジングへ
次に、逆温度のパラメータはここで見ると、行列の外にだせるので全体のスケールを調整するパラメータになってます。QUBOからイジングに直す時にはまだかんがえなくていいので、イジングmatrixに直します。
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)
<Figure size 432x288 with 1 Axes>
ぼちぼちいい感じではないでしょうか。
次に、量子ビットを指定し、バイアスと結合荷重を初期化してみます。使用するのは量子ビットを表すリストの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()
<Figure size 432x288 with 1 Axes>
まぁ、こんなもんでしょう。
アニーリング計算
この状態でアニーリングをかけてみます。アニーリングは任意の状態からスタートして、バイアスと結合荷重で結果が決まります。また、RBMの場合には近い値での局所解がたくさんありそうです。
from blueqat.wq import Opt
a = rbm().setRBM(4,4)
a.J
array([[0.35, 0. , 0. , 0. , 0.84, 0.05, 0.89, 0.82],
[0. , 0.29, 0. , 0. , 0.63, 0.27, 0.49, 0.89],
[0. , 0. , 0.65, 0. , 0.78, 0.17, 0.71, 0.93],
[0. , 0. , 0. , 0.8 , 0.71, 0.17, 0.32, 0.61],
[0. , 0. , 0. , 0. , 0.55, 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. , 0.81, 0. , 0. ],
[0. , 0. , 0. , 0. , 0. , 0. , 0.49, 0. ],
[0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.33]])
これをアニーリングにかけてみます。
#複数回を使うカウンター
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({'00000000': 100})
全体の分布
上記は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 size 432x288 with 1 Axes>
こんな感じです。逆温度パラメータはjijの前にパラメータβをつけます。これによって、分布を変えてみます。β=0.1としてみて、
beta = 0.1
c = Opt().add(beta*a.J).run(shots=100)
bars(c)
<Figure size 432x288 with 1 Axes>
これにより、より複雑な分布をとることができました。
学習について
学習に関しては、学習すべきパラメータが3種類あります。
1、visible layerのbias
2、hidden layerのbias
3、v*hの結合荷重
1は簡単でデータの入力値を使います。
2は一旦visible layerからsigmoidを通してhiddenlayerの値を計算します。
3がめんどいのですが、CD法を使う代わりにアニーリングを使って、v*hの期待値をだしてからの1/Nして平均を取ります。
学習式
学習式はバイアスと結合荷重で、
αはモーメンタムでϵは学習率です。基本的には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': 82, '00001000': 18})
0000/100 = 0000
となります。hidden layerも同様です。visible layerのデータは1100なので、更新式は、それぞれの量子ビットに対して、
一方hidden layerの更新をするためには事前にhidden layerのデータ値を出しておく必要があります。例えば、
その上で、
このようにして更新を。あとは、wに関しても同様に、データからの期待値と上記のv∗hをしたものを使って地味にひとつずつ更新をしていきます。使いやすいように更新用のmatrixを作って、
と変数を導入して計算できるように、QUBOmatrixの形にしておきます。それぞれのMにはデータの期待値とモデルの期待値をQUBOmatrixの形にしておきます。getMdata(self,vdata)としてvdataをlistでいれたら計算できるように追加しました。
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.wq 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
<Figure size 432x288 with 1 Axes>
array([[0.96, 0. , 0. , 0. , 0.5 , 0.29, 0.27, 0.77],
[0. , 0.03, 0. , 0. , 0.07, 0.25, 0.7 , 0.83],
[0. , 0. , 0.58, 0. , 0.85, 0.97, 0.15, 0.88],
[0. , 0. , 0. , 0.91, 0.35, 0.75, 0.18, 0.15],
[0. , 0. , 0. , 0. , 0.43, 0. , 0. , 0. ],
[0. , 0. , 0. , 0. , 0. , 0.83, 0. , 0. ],
[0. , 0. , 0. , 0. , 0. , 0. , 0.69, 0. ],
[0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.55]])
for i in range(10):
a.train([1,1,0,1])
Counter({'00000000': 75, '01000000': 25})
Counter({'01000000': 96, '00000000': 4})
Counter({'01000000': 100})
Counter({'01001000': 71, '01000000': 29})
Counter({'01001000': 100})
Counter({'01001000': 80, '01001100': 16, '01011001': 2, '01011000': 1, '11001100': 1})
Counter({'01011011': 46, '00010011': 21, '01010011': 19, '01000000': 6, '01011001': 3, '01001000': 2, '11010011': 1, '10010011': 1, '01010001': 1})
Counter({'00000000': 81, '01000000': 17, '01000100': 2})
Counter({'11011110': 93, '11011111': 7})
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で解分布を取ってみる
適当な例題を作ってみてやってみます。
a = rbm().setRBM(4,4)
モデルを作ります。逆温度パラメータを設定して解分布を取ります。
c = Opt().add(0.1*a.J).qaoa(shots=100)
結果は、
c.most_common(12)
(((0, 0, 0, 0, 0, 0, 0, 0), 0.904533974890349),
((1, 0, 0, 0, 0, 0, 0, 0), 0.03532018667399832),
((0, 0, 0, 0, 0, 1, 0, 0), 0.027620747673709752),
((0, 0, 0, 1, 0, 0, 0, 0), 0.00890092450276286),
((0, 0, 0, 0, 0, 0, 1, 0), 0.004226074666371333),
((0, 1, 0, 0, 0, 0, 0, 0), 0.0032522717855628675),
((0, 0, 0, 0, 0, 0, 0, 1), 0.0025049020734404194),
((1, 0, 0, 0, 0, 1, 0, 0), 0.0017011392124368463),
((1, 0, 0, 1, 0, 0, 0, 0), 0.0014345796256494785),
((1, 0, 0, 0, 0, 0, 0, 1), 0.000982528470307196),
((1, 1, 0, 0, 0, 0, 0, 0), 0.0006476938976537597),
((0, 0, 0, 0, 0, 1, 1, 0), 0.000616404779804135))
グラフにすると、
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()
<Figure size 432x288 with 1 Axes>
並びが変わりましたが、だいたい同じ結果になりました。
このようにして、量子ゲートモデルでも、確率分布をとって計算をすることができました。
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.6328815743438415),
((1, 0, 0, 0, 0, 0, 0, 0), 0.06465485757472096),
((1, 0, 0, 0, 0, 0, 0, 1), 0.02970296805595352),
((0, 0, 0, 0, 0, 1, 0, 0), 0.02216373052807483),
((1, 0, 0, 0, 1, 0, 0, 0), 0.02024237592807581),
((0, 1, 0, 0, 0, 0, 0, 0), 0.014245501009297284),
((0, 1, 0, 0, 1, 0, 0, 0), 0.013999664790584593),
((0, 0, 1, 0, 0, 0, 0, 1), 0.013501852101370603),
((0, 1, 0, 0, 0, 0, 0, 1), 0.012691519297638355),
((0, 0, 1, 0, 0, 0, 1, 0), 0.012290472014287204),
((0, 0, 0, 1, 0, 0, 1, 0), 0.012075742401000194),
((0, 0, 0, 0, 0, 0, 0, 1), 0.01041598515808674))