量子モンテカルロ法による量子アニーリングのシミュレーション

量子アニーリングとは
モンテカルロ法とは
【復習】シミュレテッドアニーリング
1. シミュレテッドアニーリングのアルゴリズム
2. 遷移確率
量子モンテカルロ法による量子アニーリング
量子アニーリングアルゴリズムのイメージ
シュレディンガー方程式とか，量子力学的なものはどこいった！？
【断熱条件】どれくらいゆっくり磁場を弱めればいいか？温度を下げればいいか？

量子アニーリングとは

量子アニーリングとは，組合せ問題最適化手法の一つです．

詳しくはこちらの記事にまとめました．

量子アニーリングのアルゴリズムはシュレディンガー方程式を解くため，量子アニーリングを量子アニーリングマシン（量子コンピュータみたいなもの）ではなく，古典コンピュータ（我々が今使っているコンピュータ）でシミュレーションするには，計算量の面から，工夫しなければなりません．

量子アニーリングを古典コンピュータで行う方法のひとつとして，量子モンテカルロ法があります．

量子アニーリングのアルゴリズムの何となくのイメージだけ知りたいひとは，後半の「量子アニーリングアルゴリズムのイメージ」の項目だけ読んでいただければOKです（目次から飛べます）．

モンテカルロ法とは

モンテカルロ法をざっと説明すると，乱数を用いてシミュレーションを何度も行い，その乱数を用いた確率で答えを出すといった計算方法です．

モンテカルロ法を用いた有名な計算例として，円の面積を求める問題があります．

下の図のように，四角形の中に乱数で発生させた点を打ち，（円の中に入った点の数）/（四角形の中の点の数）×（四角形の面積）で円の面積を近似的に求めることができます．

モンテカルロ法と言っても色々あるので，とりあえず「乱数を使ってシミュレーションするんだな」くらいに思っておいてください．

【復習】シミュレテッドアニーリング

量子アニーリングを説明する前に，よく似た最適化手法であるシミュレテッドアニーリング（SA）について述べておきます．

シミュレテッドアニーリングもモンテカルロ法の一種です（正確には，マルコフ連鎖モンテカルロ法を利用した最適化手法）．

シミュレテッドアニーリングのアルゴリズム

シミュレテッドアニーリングは，はじめのうちは色々な解を探索して，だんだんと良い解のみを探索していき，最適解を見つける方法です．

イメージとしては，以下のような図です．

シミュレテッドアニーリングのアルゴリズムは，以下の通りです．量子アニーリングも流れは同じなので，覚えておいてください．

シミュレテッドアニーリングのアルゴリズム

１．温度の初期値を決定

２．解を遷移させる（遷移確率については以下で説明）

３．一定回数２．を試行したら温度を下げる

４．温度が十分小さくなったら終了．そうでなければ，２．に戻る

つまり，解の候補を選びながら温度を下げていきます．

ここでカギとなるのが，２．で行う「どのように遷移させるか（解の候補を選ぶか）？」です．

新しい解の候補として何かを選んだとき，前の解（暫定解）での目的関数と比べて目的関数が小さくなればその新しい解を暫定解とすれば良さそうですが，それだと局所解（先ほどの図の小さな谷）で終わってしまいます．

我々が見つけたいのは，大きな谷のそこを示す解です．

そこで，遷移確率というものを定義して，新しい候補解での目的関数が暫定解のときの目的関数よりも大きくなったとしても（悪化したとしても），ある確率で遷移するようにします．

遷移確率

遷移確率として最も一般的に使われるのは，

$P_{SA} \propto \exp\left(-\frac{\Delta E}{T}\right)$

です．ここで， $\Delta E$ は目的関数の増加量（改悪度，（新しい解の目的関数値）-（前の解の目的関数値））， $T$ はその時の温度です．

つまり，良い解になる（目的関数値が小さくなる）ような遷移は起きやすく，また，温度が高いほど遷移確率が高くなります．

ただし，温度が高い時には，目的関数値が大きくなるようなときにも遷移する可能性があり，これによって局所的な最適解を避けることができます．

さて，上式のような遷移確率はどこからきているのでしょうか？いろいろな解候補を見ながら最適解を探すだけなら，このような遷移確率ではなくてもいいような気もします．

つまり，目的関数が減少するときと温度が高い時には遷移確率が大きく，温度が０になったときには遷移しないような関数を遷移確率としてもってこれば良さそうです．

実は，この式は「統計力学」をモデルとして使っており，一応理由づけされた遷移確率なのです．

この統計力学の考え方を用いて，定常分布（遷移しても同じ状態になるときの分布）がギブス・ボルツマン分布になるように遷移確率が設定されています（なぜギブス・ボルツマン分布なのかというと，この分布は物理に照らし合わせると最低エネルギーをとる基底状態だからです）．

ギブス・ボルツマン分布は，

$P_{GB} = \frac{1}{Z}\exp(-\beta H(\sigma))$

と書かれます．

$\beta$ は逆温度で， $\beta=\frac{1}{T}$ です． $H(\sigma)$ はエネルギー（量子論の場合はハミルトニアン演算子）です．

状態 $\sigma'$ から $\sigma$ への遷移確率を $M(\sigma|\sigma')$ とするとき，ギブス・ボルツマン分布が定常分布となるにはつり合い条件

$P_{GB} = \sum_{\sigma'}M(\sigma|\sigma')P_{GB}(\sigma')$

を満たす必要があり，これを満たす遷移確率 $M(\sigma|\sigma')$ は，例えば

$M(\sigma|\sigma') = \frac{e^{-\beta H(\sigma)}}{e^{-\beta H(\sigma)}+e^{-\beta H(\sigma')}}$

となります．これをすこし変形すれば，先ほどの遷移確率と一致するわけです．

この一連の理論は，マルコフ連鎖モンテカルロ法と呼ばれます（乱数は新しい解候補を探すときに使います）．

このような手続きは，一つ一つの状態（解）を運動方程式を解くことなく，統計的に扱っていることを意味しています．

つまり，「ちょっと粗いけど，マクロ的に見れば大体こんな感じになるよね」という立場の熱力学や統計力学を用い，エネルギーが最小化される物理事象をモデルに，最適化しようというわけです．

同様に，統計力学やモンテカルロ法の考えに沿うことで，シュレディンガー方程式を解かず，近似的に解の状態を追っていくことで，量子アニーリングを行う方法を以降で述べていきます．

量子モンテカルロ法による量子アニーリング

量子アニーリングを量子モンテカルロ法で行うときの考え方は

「シュレディンガー方程式をいちいち解くのではなく，量子統計的に解きましょう」というものです．

量子アニーリングを実行するにはその時々のハミルトニアンが分かればいいので，量子統計で出てくる分配関数からエネルギーを求めればOKです．

ただし，量子的な分配関数そのままでは古典的なシミュレーションができないので，古典的分配関数の形に直します．

量子モンテカルロ法の考え方

・シュレディンガー方程式を量子統計的に表現して計算する．

・量子アニーリングの量子的な分配関数を古典的な分配関数と同じ形に書き直す．

・古典的分配関数が分かれば，エネルギー（ハミルトニアン）が分かり，古典コンピュータで量子アニーリングが実行できる．

したがって，以降では

①量子アニーリング用イジング模型の分配関数を古典的な形で求める，

②その分配関数を見れば古典イジング模型としてのハミルトニアンは何かが分かる，

③そのハミルトニアンを使って量子アニーリングを実行する

の順で説明していきます．

そこそこ細かく書いてあるので，結論だけ知りたい人は，最後らへんだけ見て頂ければと思います．

量子力学を統計力学で扱うときは，分配関数 $Z$ は

$Z = \rm Tr \left[e^{-\beta \hat H}\right]$

です．

分配関数とは，分布の規格化定数の逆数であり，ギブス・ボルツマン分布の場合は

$P_{GB} = \frac{1}{Z}\exp(-\beta H(\sigma))$

の $Z$ が分配関数であり，古典力学の場合， $Z = \sum_{\sigma}\exp(-\beta H(\sigma))$ です．

さて，シミュレテッドアニーリングのときと同様にして，量子アニーリングの遷移確率を

$P_{QA} = \frac{1}{Z}<\sigma|e^{-\beta H}|\sigma>$

と書きます．あとはこれを計算できる形にしていくだけですが，この量子アニーリングにおいて，最小化すべきハミルトニアン及び遷移確率は分配関数からすぐに求められるので，分配関数を変形していきます．

つまり，

$Z= \rm Tr \left[e^{-\beta \hat H}\right]=\sum_{\sigma}<\sigma|e^{-\beta \hat H}|\sigma>$

を計算できる形に変形していきます．

ここで， $\hat H$ を復習しましょう．量子アニーリングにおいて，ハミルトニアン $\hat H$ はイジング模型のスピン相互作用による古典ハミルトニアン $\hat H_c$ と横磁場による量子ハミルトニアン $\hat H_q$ の和だったので，

$\hat H = \hat H_c + \hat H_q$

です．ハミルトニアンの内約は，＜↑｜猫でもわかる量子アニーリングの理論｜↓＞を見てください．

ここでは， $\hat H=\hat H_c + \hat H_q = \hat H_c-\Gamma\sum_{i}^{N}\hat\sigma_i^x$ として話を進めていきます．

$\hat H_c=-\sum_{i,j}J_{ij}\hat \sigma^z_i \hat \sigma^z_j$ です．

$<\sigma|e^{-\beta \hat H}|\sigma>=<\sigma|e^{-\beta (\hat H_c + \hat H_q)}|\sigma>$ を考えたとき，もしハミルトニアンが古典のみで $<\sigma|e^{-\beta \hat H_c}|\sigma>$ だけだったら， $\hat H_c$ は対角行列なので $<\sigma|e^{-\beta \hat H_c}|\sigma>=e^{-\beta E_c(\sigma)}$ と簡単に計算できるのですが， $\hat H_q$ は対角行列ではなく，かつ $\hat H_c$ と $\hat H_q$ は非可換であるため， $e^{-\beta \hat H} = e^{-\beta \hat H_c}e^{-\beta \hat H_q}$ と分割することはできず，単に $<\sigma|e^{-\beta \hat H}|\sigma>=<\sigma|e^{-\beta (\hat H_c + \hat H_q)}|\sigma>=e^{-\beta E(\sigma)}$ のように計算することはできません．

この問題を解決するには，『鈴木・トロッタ分解』を使用します．

鈴木・トロッタ分解による経路積分表示

鈴木・トロッタ分解を用いると，以下のように $e^{\hat H} = e^{-\beta (\hat H_c+\hat H_q)}$ を分解することができます．

$e^{-\beta (\hat H_c+\hat H_q)}=\lim_{M \to \infty}\left(e^{-\frac{\beta}{M}\hat H_c}e^{-\frac{\beta}{M}\hat H_q}\right)^M$

鈴木・トロッタ分解の証明は省略しますが，指数関数を $M^{-1}$ で展開したときに $M$ の2次以降の項が $M \rightarrow \infty$ で十分小さく無視できることを用いています．

これを用いて，分配関数を以下のように変形していきます．

$\begin{split} Z &= \sum_{\sigma}<\sigma|e^{-\beta \hat H}|\sigma> \\&= \lim_{M \to \infty}\sum_{\sigma}<\sigma|\left(e^{-\frac{\beta}{M}\hat H_c}e^{-\frac{\beta}{M}\hat H_q}\right)^M|\sigma> \end{split}$

ここで，とある状態 $\sigma_l$ に関して，単位行列 $I$ を用いて

$I = \sum_{\sigma_l}|\sigma_l><\sigma_l|$

$(l = 1,...,M)$

が成り立つことを用いると（※ $|\sigma_l>$ はスピン $l$ の状態ではなく，N個のスピンのある状態を表す），

$\begin{split} Z&= \lim_{M \to \infty}\sum_{\sigma}<\sigma|e^{-\frac{\beta}{M}\hat H_c}e^{-\frac{\beta}{M}\hat H_q}\cdot I \cdot e^{-\frac{\beta}{M}\hat H_c}e^{-\frac{\beta}{M}\hat H_q}\cdot I \cdot ... e^{-\frac{\beta}{M}\hat H_c}e^{-\frac{\beta}{M}\hat H_q} |\sigma> \\&=\lim_{M \to \infty}\sum_{\sigma_1,...,\sigma_M}<\sigma_1|e^{-\frac{\beta}{M}\hat H_c}e^{-\frac{\beta}{M}\hat H_q}|\sigma_M><\sigma_M| e^{-\frac{\beta}{M}\hat H_c}e^{-\frac{\beta}{M}\hat H_q}|\sigma_{M-1}> ... \\&...<\sigma_2|e^{-\frac{\beta}{M}\hat H_c}e^{-\frac{\beta}{M}\hat H_q} |\sigma_1> \\&=\lim_{M \to \infty}\sum_{\sigma_1,...,\sigma_M}<\sigma_1|e^{-\frac{\beta}{M}\hat H_q}|\sigma_M>e^{-\frac{\beta}{M}\hat H_c(\sigma_M)}<\sigma_M| e^{-\frac{\beta}{M}\hat H_q}|\sigma_{M-1}>e^{-\frac{\beta}{M}\hat H_c(\sigma_{M-1})} ... \\&...<\sigma_2|e^{-\frac{\beta}{M}\hat H_q} |\sigma_1>e^{-\frac{\beta}{M}\hat H_c(\sigma_1)}\end{split}$

となります．

ここで2行目から3行目では $\hat \sigma^z$ が対角化された表示での和が $\sigma_l$ についての和であることを使いました．より詳細に書けば，

$<\sigma_l|e^{-\frac{\beta}{M}\hat H_c}|\sigma_l>=<\sigma_l|e^{\frac{\beta}{M}J_{ij}\hat \sigma^z_i \hat \sigma^z_j}}|\sigma_l>=e^{\frac{\beta J}{M}\sigma_{i,l} \sigma_{j,l}}}$

となっています． $\sigma_{i,l}$ や $\sigma_{j,l}$ は±１をとる古典イジング変数です．

ここで，イジング変数に $l$ という次元が加わっていることがわかります．

つまり，古典イジン模型に書き直すと，次元が1つ増えるのです．これをトロッタ次元と呼びます．

この $l$ の意味ついては，後の「トロッタ次元「 $l$ 」はどんな意味をもつのか（量子力学的説明）」や「1次元量子統計学的イジング模型は2次元古典イジン模型と同じ」のところで説明します．

さて，分配関数の変形の続きをしていきましょう．

$(\hat \sigma_j^x)^2=1$ であることを使えば，

$<\sigma_{l+1}|e^{-\frac{\beta}{M}\hat H_x}|\sigma_{l}>=<\sigma_{l+1}|(\cosh h+\sinh h \hat \sigma_j^x)|\sigma_{l}>=A(h)e^{K^x(h)\sigma_{j,l}\sigma_{j,l+1}}$

という形に書くことができます．

ここで， $h=\frac{\beta \Gamma}{M}$ を用いて

$K^x(h)=-\frac{1}{2}\log \tanh{h}$

$A(h)^2=\frac{1}{2}\sinh(2h)$

です．

今までのことをまとめると，量子統計的横磁場イジング模型の分配関数は，次のような古典イジング模型の分配関数で表現できることになります．

$Z=\rm Tr \left[e^{-\beta \hat H}\right]=\lim_{M \to \infty}A(h)^{MN}\sum_{\sigma}e^{-\beta H_c(\sigma)}$

そして，この時のハミルトニアンは

$H_c(\sigma)=-\sum_{l=1}^{M}\sum_{j=1}^{N-1}\left( \frac{\beta J}{M}\sigma_{j,l} \sigma_{j+1,l}}+\frac{K^x}{\beta}\left(\frac{\beta\Gamma}{M}\right)\sigma_{j,l} \sigma_{j,l+1}} \right)$

であり，目標のハミルトニアンが得られました．このハミルトニアンについては後で少し考察します．

さて量子アニーリングを実行するには，横磁場を弱めていきながら状態を求めていきますが，先ほど得られたハミルトニアンから決められた転移確率をもとに一定の確率でスピン上下を反転させていけばOKです．

つまり，量子アニーリングを実行するアルゴリズムは

量子アニーリングのアルゴリズム

１．磁場の初期値を決定．温度パラメータに比べて十分大きくする．

２．遷移確率 $P_q = min(1,e^{-\beta \delta H})$ に従って解を転移させる

３．一定回数２．を試行したら磁場を弱める

です． $\delta H$ は先ほど得られたハミルトニアンの変化量です．※具体的な形は，後日改めて書きます．

これを見て分からないところがある人は，＜↑｜猫でもわかる量子アニーリングの理論｜↓＞を読んでください．

1次元量子統計学的イジング模型は2次元古典イジン模型と同じ

1次元の横磁場イジング模型の分配関数を古典イジング模型の分配関数で表現すると，対応するハミルトニアンは以下のようになることが分かりました．

横磁場イジング模型を古典イジング模型に書き直したときのハミルトニアン

$H_c(\sigma)=-\sum_{l=1}^{M}\sum_{j=1}^{N-1}\left( \frac{\beta J}{M}\sigma_{j,l} \sigma_{j+1,l}}+\frac{K^x}{\beta}\left(\frac{\Gamma}{M}\right)\sigma_{j,l} \sigma_{j,l+1}} \right)$

元の横磁場イジング模型ではスピン番号 $i$ という1つの次元（＝スピンを横一直線に並べた方向）だけを持っていましたが，古典イジング模型で表現したら， $l$ という次元が増えています．

このハミルトニアンは，2次元の古典イジング模型です．つまり，古典系に直すと以下の図のように変換されていることになります．

実際に解いているのは横磁場ありの1次元イジング模型ですが，数学的には2次元イジング模型となっているわけですね．

トロッタ次元「l」はどんな意味をもつのか（量子力学的説明）

さて，古典イジング模型表記にすると次元が１つ増えた，というのは分かりましたが，「 $l=1, ... ,M$ とは何なのか？」という疑問が残ります．

もちろん，実際にN×M個のスピンに増えたわけではありません．

これは，量子力学の経路積分と対比すれば，虚時間方向です．経路積分では色々な経路で作用を位相にもつ因子を虚時間方向に足し上げますが，それと対応しています．

つまり， $l$ は色々な経路を指します．

量子アニーリングアルゴリズムのイメージ

シミュレテッドアニーリングは，以下のように一つの解の可能性から始めて，アルゴリズムによって試行を進め（iteration）ます．

一方，量子アニーリングでは $l$ 次元方向（虚時間方向，色々な経路）に $M$ パターン用意して，互いの経路が干渉しながら試行を進めていきます．

経路同士の干渉具合は，各径路の解になる確率で結び付けられ，各径路の解が実現する確率の和は１になります．

これは，＜↑｜猫でもわかる量子アニーリングの理論｜↓＞で説明した量子アニーリングと非常に似ていることが分かりますね（下図）．

詳細は省略しますが，経路の重ね合わせ度は以下の $f$ のような形を取ります（名前は筆者が適当に付けたものです）．

重ね合わせ度が小さい時ときには，たくさんの経路が実現しようとします．一方，重ね合わせ度が大きい時には，実現経路は一つになろうとします（下図）．

これを用いると，以下のような関係があることが分かります．

磁場 $\Gamma$ が小さい時には $f$ が大きくなりますから，経路はほとんど1つになります．経路が一つということは，シミュレテッドアニーリング（SA）と実質的に等しくなります．

一方で，磁場が大きい時にはたくさんの経路の解が混合されます．つまり，磁場が大きい時には量子力学っぽくなります．

この経路は，量子ゲート方式量子コンピュータの各ビットの経路と似ていますね（これが対応しているかは分かりませんが，実際，量子アニーリングと量子ゲートは実質的に等しいことが証明されているらしいです）．

シュレディンガー方程式とか，量子力学的なものはどこいった！？

あまり量子統計力学を学んだことがないひとは，

「いやいや，シュレディンガー方程式とか，量子力学とかどこいった！？シュレディンガー方程式出てきてない！ただの分配関数やん！シミュレテッドアニーリングとどう違うの！？」

と思うかもしれません．

しかし，経路積分表示は量子力学の別の記述方法で，シュレディンガー方程式と同じ意味を持ちます．

経路積分表示とシュレディンガー方程式が等価である証明は，以下をご覧ください。

https://www.phys.chuo-u.ac.jp/labs/nakano/kaisekirikigaku/sec9(2010).pdf

つまり，量子モンテカルロ法は経路積分表示で量子力学を扱っているだけで，本質は保存されています．

それに，単に分配関数を記述しているのではなく，古典コンピュータで計算できる形で書いているので，先ほど説明した通り次元が1つ増えた形に変わっており，そのままというわけではありません．

【断熱条件】どれくらいゆっくり磁場を弱めればいいか？温度を下げればいいか？

量子アニーリングにしろ，シミュレテッドアニーリングにしろ，適切な解を得るには制御パラメータ（磁場や温度）を「ゆっくり」小さくしていかなければなりません．

そこで，どれくらい「ゆっくり」ならいいの？という疑問が残ります．

これについては，以下の記事に書きましたのでぜひ参考にしてみてください．

fit より:
6月 11, 2021 7:55 pm
量子アニーリングを量子モンテカルロ法で実装を試みていたので大変参考になりました．
その上で一つ質問をさせていただきたいのですが，温度やガンマを具体的にどのくらいの値に設定もしくは変動させるのが良いのでしょうか？
返信
- 理系リアルタイムより:
  6月 12, 2021 1:29 am
  fit様
  ご質問ありがとうございます．
  温度やガンマは，シミュレテッドアニーリングと同様，人の感覚で決定することが多いのが現状だと思います．
  また，問題設定によって値は変わってくると考えられます．
  これらのパラメータを機械学習などの手法で最適化する研究も行われていたりしますが，個人的にはまだ難しい印象を持っています．
  また，本記事の内容は量子アニーリングのモンテカルロシミュレーションの最も基本部分しか話しておらず，現在では色々な手法が提案されているので，パラメータの設定はさらに複雑化していると思われます．
  ちなみに，以下の記事には収束条件を記述しています．
  https://atsblog.org/qa_shusoku/
  記事をご覧いただき，ありがとうございました．
  返信