D-Waveの量子ボルツマンマシンの逆温度パラメータ最適化でPFNのOptunaつかってみた


はじめに

量子コンピュータや量子アニーラを触っていると組合せ最適化問題や機械学習分野において多数のパラメータ調整にであいます。前回は私たち業務に携わるものとして喫緊の課題であった量子ゲートモデルでのVQEやQAOAでのOptuna利用を切実に検討してみました。

量子コンピュータゲートモデルの量子古典ハイブリッド計算のVariational Quantum Eigensolverの古典パラメータ最適化にPFNのOptuna使ってみた。

https://qiita.com/YuichiroMinato/items/afee437cad318699593e

今回はそれとはうってかわって量子アニーラでの切実なイジングの制約条件パラメータの最適化、、、ではなく全く需要のなさそうな量子ボルツマンマシンのRBMの確率分布を求める際の逆温度パラメータの最適化にチャレンジしてみたいと思います。検討はカナダのD-WaveマシンとMDRのWildqatSDKを使います。


D-Waveとは?

カナダのベンチャー企業で、量子アニーリングと言われる量子効果を活用したアルゴリズムをハードウェア実装したマシンです。

企業情報:

本社所在地:カナダブリティッシュコロンビア州バーナビー

設立:1999年

事業内容:量子コンピュータのハードウェアとソフトウェアを商用製品として提供

従業員数:160名(Ph.D 55名)

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

https://www.dwavesys.com/

日本語サイトもあります。

http://dwavejapan.com/


今回D-Waveで何をやるのか?

ずばり量子ボルツマンマシンでのRBM(制限付きボルツマンマシン)の学習用の確率分布を正しく求めたいという動機です。

簡単に確認すると、、、下記のような2層の制限付きボルツマンマシンのモデルがあった際、

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

こちらのモデルは、結合分布がボルツマン分布に従い、目的関数としてエネルギー関数を持ちます。

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

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のモデルの学習は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)


ここでD-Wave

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)}

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

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

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}]


何を最適化するのか

求めたい確率分布はβできまります。本来は逆温度パラメータなのでアニーリングの温度を調整したいところですが、それはできないので擬似的に調整しています。本来は学習がきちんと行われているように確認をしますが、今回は小さな問題で実際の厳密解から得られる確率分布に対して適切なβを設定するように最適化するトイモデルを確認したいと思います。


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

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

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

1.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の場合の数でエネルギーを求めて見ると、

1.png

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

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がでます。

2.png

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


まずはWildqatで

Wildqatで自動的にこれを求めたい!実装してみます。

まずは自動化しないでサンプリング。

from wildqat import * 

import numpy as np
import matplotlib.pyplot as plt

def binarycount(binarr):
if(binarr[0]==0 and binarr[1]==0):
return 0
if(result[0]==0 and result[1]==1):
return 1
if(result[0]==1 and result[1]==0):
return 2
if(result[0]==1 and result[1]==1):
return 3

a = opt()
beta = 0.03

a.qubo = beta * np.array([[-1,2],[0,-1]])
count = []

for i in range(20):
result = a.sa()
count.append(binarycount(result))

fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.hist(count,bins=8)
fig.show()

2-3.png

結果は良好。なぜならパラメータ調整した後だからです。

コード中にbetaというパラメータを0.03にしていますが、これを知らない状態で最適化したいです。

ちなみにbeta = 1としてやってみると、、、

3.png

このようにばらつきは出ません。うまく係数を最適化する必要があります。

ということで最適化しました。

from wildqat import * 

import numpy as np
import optuna
from scipy import stats

def binarycount(binarr):
if(binarr[0]==0 and binarr[1]==0):
return 0
if(binarr[0]==0 and binarr[1]==1):
return 1
if(binarr[0]==1 and binarr[1]==0):
return 2
if(binarr[0]==1 and binarr[1]==1):
return 3

def objective(trial):
a = opt()
a.R = 0.5
beta = trial.suggest_uniform('x',0.0001,1.0)
a.qubo = beta * np.array([[-1,2],[0,-1]])
count = [0,0,0,0]

for i in range(100):
result = a.sa()
count[binarycount(result)] += 1

val1 = [13,37,37,13]
print(beta,count)
return entropy(val1,count)

study = optuna.create_study()
study.optimize(objective, n_trials=100)

やっぱりいい感じでパラメータを最適化してくれました。結構面倒だなと思っていたので大変助かります。

(略)

[I 2018-12-05 23:17:17,454] Finished a trial resulted in value: inf. Current best value is 0.0027771501937655374 with parameters: {'x': 0.046462775556521}.
0.043734352042027934 [7, 38, 38, 17]
[I 2018-12-05 23:17:36,767] Finished a trial resulted in value: 0.025866275994701357. Current best value is 0.0027771501937655374 with parameters: {'x': 0.046462775556521}.
0.002543606466022884 [18, 24, 30, 28]
[I 2018-12-05 23:17:54,784] Finished a trial resulted in value: 0.09570822500371039. Current best value is 0.0027771501937655374 with parameters: {'x': 0.046462775556521}.
0.000988252005248394 [23, 23, 29, 25]
[I 2018-12-05 23:18:13,409] Finished a trial resulted in value: 0.10686566600428929. Current best value is 0.0027771501937655374 with parameters: {'x': 0.046462775556521}.
0.3628006425173081 [0, 45, 55, 0]
[I 2018-12-05 23:18:31,934] Finished a trial resulted in value: inf. Current best value is 0.0027771501937655374 with parameters: {'x': 0.046462775556521}.
0.004903371809014956 [23, 27, 21, 29]
[I 2018-12-05 23:18:51,355] Finished a trial resulted in value: 0.14767043994842627. Current best value is 0.0027771501937655374 with parameters: {'x': 0.046462775556521}.
0.21284808355208942 [0, 55, 45, 0]
[I 2018-12-05 23:19:11,059] Finished a trial resulted in value: inf. Current best value is 0.0027771501937655374 with parameters: {'x': 0.046462775556521}.
0.03954824784065314 [12, 34, 42, 12]
[I 2018-12-05 23:19:30,075] Finished a trial resulted in value: 0.00519920647901957. Current best value is 0.0027771501937655374 with parameters: {'x': 0.046462775556521}.
0.9611240639846177 [0, 48, 52, 0]
[I 2018-12-05 23:19:49,451] Finished a trial resulted in value: inf. Current best value is 0.0027771501937655374 with parameters: {'x': 0.046462775556521}.
0.011988166865586326 [22, 24, 39, 15]
[I 2018-12-05 23:20:05,338] Finished a trial resulted in value: 0.05368631692014822. Current best value is 0.0027771501937655374 with parameters: {'x': 0.046462775556521}.


つぎにD-Waveで

D-Waveではきちんと設定すればサンプリング機能がありますが、しょうもなく毎回呼び出しちゃいましたがやってみたところ、

from wildqat import * 

import numpy as np
import optuna
from scipy import stats

def binarycount(binarr):
if(binarr[0]==-1 and binarr[1]==-1):
return 0
if(binarr[0]==-1 and binarr[1]==1):
return 1
if(binarr[0]==1 and binarr[1]==-1):
return 2
if(binarr[0]==1 and binarr[1]==1):
return 3

def objective(trial):
a = opt()
a.dwavetoken = "token here"
beta = trial.suggest_uniform('x',0.001,1.0)
a.qubo = beta * np.array([[-1,2],[0,-1]])
count = [0,0,0,0]

for i in range(10):
result = a.dw()
count[binarycount(result)] += 1

val1 = [1.3,3.7,3.7,1.3]
print(beta,count)
return stats.entropy(val1,count)

study = optuna.create_study()
study.optimize(objective, n_trials=10)

きちんとパラメータが最適化されました。お任せでいいので楽ですね。。。

0.6042506628927251 [1, 3, 5, 1]
[I 2018-12-05 23:46:19,987] Finished a trial resulted in value: 0.03440242089486224. Current best value is 0.03440242089486224 with parameters: {'x': 0.6042506628927251}.
0.673443366429525 [4, 1, 2, 3]
[I 2018-12-05 23:46:41,130] Finished a trial resulted in value: 0.45687867402306015. Current best value is 0.03440242089486224 with parameters: {'x': 0.6042506628927251}.

少なくともパラメータレンジのあたりはつくので、計算量で済ませられるのはとても助かります。

以上です。