@yuichirominato 2018.12.05更新 271views

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

D-Wave RBM イジング パラメータ最適化 量子アニーリング

はじめに

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

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

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

D-Waveとは?

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

企業情報:
本社所在地:カナダブリティッシュコロンビア州バーナビー
設立:1999年
事業内容:量子コンピュータのハードウェアとソフトウェアを商用製品として提供
従業員数:160名(Ph.D 55名)

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

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

詳細はカナダのホームページを参照してくださいませ。
https://www.dwavesys.com/

日本語サイトもあります。
http://dwavejapan.com/

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

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

興味ある方は下記を参照していただきまして、、、

D-Waveで深層学習の基礎となるRBMのボルツマン学習を実行してみた
http://blog.mdrft.com/post/258

簡単に確認すると、、、下記のような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)$

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

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

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

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

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

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

何を最適化するのか

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

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

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

キメラグラフ参考資料:https://qiita.com/YuichiroMinato/items/e6952fec1a9965156873

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

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

例題

簡単な例題で上記の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$の場合の数でエネルギーを求めて見ると、

解析的にかんがえて見ると、ここではすべての場合の数が出せているので、まず規格化定数$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がでます。

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

まずは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()    

結果は良好。なぜならパラメータ調整した後だからです。
コード中にbetaというパラメータを0.03にしていますが、これを知らない状態で最適化したいです。
ちなみにbeta = 1としてやってみると、、、

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


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

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

あわせて読みたい

SERVICE

汎用量子コンピュータ向けアプリケーション開発SDK

詳しく見る Githubで入手する(無料)

汎用量子コンピュータ向け高速シミュレータ

詳しく見る

量子コンピュータ向けクラウドサービス(準備中)

詳しく見る

イジングマシン向けアプリケーション開発SDK

詳しく見る Githubで入手する(無料)

COMMUNITY

量子コンピュータのことを学びたい、仕事としたいなどの情報交換の場を設け、コミュニティの形成を進めています。オフラインの勉強会と、オンラインのチャットコミュニティの2種類あります。オフラインのConnpassは1400名、オンラインのSlackは880名を超える参加があります。どちらも無料ですのでお気軽にご参加ください。

CONNPASS SLACK

CONTACT

弊社あての連絡はinfo@mdrft.comより用件を明記の上、メールで連絡を頂けますと幸いです。

ブログトップに戻る

Recent Posts