量子アニーリング - BLOG by Hideyuki Miyahara

　組み合わせ最適化問題を数値的に解く手法として、シミュレーテッド・アニーリング（SA）と呼ばれる手法が古くから知られており、広く使われている。（量子アニーリングとの比較で古典アニーリング（CA）と呼ばれることもある。）

　一方で量子アニーリングと呼ばれる手法がイジング模型という統計力学の模型を解く手法として、１９９８年にKadowaki、Nishimori[1]によって提案された。この手法では、古典的模型である”イジング模型”を量子的模型である”横磁場イジング模型”に変更することでエネルギー最適化問題を解く。SAとQAの大きな違いの一つに評価関数の局所解から局所解への遷移確率が異なるということがあげられる[2]。これについては後で簡単に触れるが、特に評価関数が複雑な場合はQAは有効な手法になり得る。さらに、QAの適用範囲は広い。例えば、有限次元のグラフとして定義される最適化問題の評価関数はイジング模型のハミルトニアンのようにみなせるので、QAは適用可能である。

　上記の理由から多くの適用先があると考え、QAに関する論文[1, 2, 3, 4, 5, 6, 7, 8, 11, 12]を10本程度読んだ。（ただ何を読んだか忘れてしまったので、随時更新する。）論文以外にネットですぐ見つかるものとして解説記事[9]を読んだ。

　QAの概要に入る前に簡単にSAの概要を述べる。SAで解くイジング模型は以下のように定義される。

$H=-{\sum}_{ij} J_{ij} {\sigma}^z_i {\sigma}^z_j$

一般には $H$ はハミルトニアンと呼ばれるエネルギー演算子である。ただ古典的模型では演算子である必要はないので、物理量としてのエネルギーである。最適化の文脈では評価関数である。定義を見てみると、右辺の初めに ${\sum}_{ij}$ とあるが、これは隣接したサイトに関しての合計である。 ${\sigma}^z_i$ はサイト $i$ におけるz軸方向のスピンである。本来スピンはパウリ行列で表現されるが、古典的なイジング模型を考える場合、 ${\sigma}^z_i=\pm 1$ で表現できる。

　SAは上記のイジング模型におけるエネルギーを下げるための乱数アルゴリズムである。初期値として各サイトのスピンをランダムに設定する。1回のステップで各サイトのスピンを反転させ、その前後でエネルギーを計算し、ある確率でスピンの反転を受容あるいは拒絶する。スピン反転させる前のエネルギーを $E_1$ とし、スピン反転後のエネルギーを $E_2$ とすると、そのステップでのスピン反転を受容する確率を

$p=\frac{1}{1+\exp(-\frac{E_2-E_1}{k_b T})}$

と定義する。ここで $k_b$ はボルツマン定数、 $T$ は温度を表す。初め温度 $T$ を高温に設定し、スピン反転を繰り返すたびに温度を下げていく。あらかじめ温度やエネルギーのばらつきに閾値を設けて、計算を終了する。これはざっくりと説明すると以下の2点を行っている。

（１）初めは温度を高温に設定し、エネルギーが上昇する（評価関数が悪くなる）スピン反転を受容する確率を高くすることで、局所解に陥りにくくする。

（２）徐々に低温にすることで局所解を探すという戦略に切り替える。

　このアルゴリズムでは局所解から局所解へ飛び移る確率が

$p \approx \exp(-\frac{\Delta }{k_b T})$

と計算される[2, 12]。ここで $\Delta$ とはポテンシャル障壁の高さである。

　本題のQAの概要である。QAは古典的イジング模型（z軸スピンのみ）にx軸方向（y軸方向でもよい）に磁場をかけ、アニーリングの過程でその磁場を徐々に弱くし、最終的には0にする。ハミルトニアンが以下のように変更し、 $\Gamma$ を最終的に0にする。

$H=-{\sum}_{ij} J_{ij} {\sigma}^z_i {\sigma}^z_j+\Gamma {\sum}_i {\sigma}^x_i$

${\sigma}^{x}_i$ を含む項が付け加わっていることが分かる。 ${\sigma}^{x}_i$ はサイト $i$ におけるx軸方向のスピンであり、 $\Gamma$ は地場の強さである。この時点ではじめの古典的模型は量子的模型になっている。先に触れたが、物理の用語でいうところの”横磁場イジング模型”である。スピン ${\sigma}^{z(x)}_i$ は以下のようにパウリ行列で表現される。

${\sigma}^z_i = \begin{pmatrix} 1 \ \ \ \ \ \ 0 \\ 0 \ \ -1 \end{pmatrix}, \ \ {\sigma}^x_i = \begin{pmatrix} 0 \ \ 1 \\ 1 \ \ 0 \end{pmatrix}$

このときイジング模型のときのように ${\sigma}^{z(x)}_i=\pm 1$ とは表現できない。 ${\sigma}^{z(x)}_i$ が交換しないことが量子的模型として本質的である。

　しかしこのような模型をどのように解けばいいのかという問題が生じる。それに対する解が以下である。量子的模型は鈴木トロッター展開を行うことで古典的模型にマップすることができる。重要なことはトロッター展開によって最初に定義された模型は1次元高い模型になることである。今回の場合、3次元の横磁場イジング模型（量子的模型）が4次元のイジング模型（古典的模型）にマップされる。実際にトロッター展開を行うと、新しい次元（トロッター次元、トロッター方向などと言う）はもともとの模型を並べて、同じスピン間に

$J=-\frac{PT}{2}\log \tanh \frac{\Gamma}{PT}$

というカップリングを導入することで定義される。ここで $P$ はトロッター方向のレプリカの数、 $T$ はシステムの温度である。計算は参考文献[10]に詳細が示されている。わかりやすい論文なので参照していただきたい。結果として先の”横磁場イジング模型”は

$H=-{\sum}_{ij} J'_{ij} {\sigma}^z_i {\sigma}^z_j$

のように表現することができる。 $J'$ は上記の新しい次元を考慮して $J$ を取り直している。後はこれをSAと同じ要領で解くだけである。

　実際には、以下のように $\tau$ 、 $t$ を導入してパラメータを取り直す。

$H=-\frac{t}{\tau}{\sum}_{ij} J_{ij} {\sigma}^z_i {\sigma}^z_j+(1-\frac{t}{\tau})\Gamma {\sum}_i {\sigma}^x_i$

$\tau$ がアニーリングパラメータで、 $0$ からスタートして、 $t$ で終了する。単純に温度は初めから低温に固定する。そのため断熱的であると表現される。当然であるが、温度もアニーリングパラメータとして問題ない。この場合はアニーリングパラメータが複数になるので、アルゴリズムに修正が必要である。

　このアルゴリズムでは局所解から局所解へ飛び移る確率が

$p \approx \exp(-\frac{\sqrt{\Delta}w }{k_b T})$

と計算される[2, 12]。ここで $\Delta$ とはポテンシャル障壁の高さであり、 $w$ はポテンシャル障壁の壁の厚さである。この手法は相対的にポテンシャル障壁の高さが高く、その壁の厚さが薄いときに局所解に止まりにくいということであり、評価関数が複雑であるときに局所解に止まりにくいということを意味している。これはまさにトンネル効果であり、モデルを量子化したことによる恩恵である。この点は非常に重要であり、QAの優位性を述べている。

　ここで一点SAとQAに共通する重要なことを付け加える。SAとQAにおいて最適化されるのは実際にはハミルトニアンあるいはエネルギーではない。自由エネルギーである。自由エネルギーは以下のように定義される。

$F=H-TS$

ここで $S$ はエントロピーである。これは両者に通じて温度を下げなければいけないもう一つの理由を示している。温度が高いとエントロピーが大きい状態が計算結果として得られてしまうが、温度を下げればエネルギーと自由エネルギーは近づくからである。

　次に数値計算についてである。厳密にはトロッター方向に無限にレプリカを並べなければいけないが、数値計算上は不可能である。実際、数値計算上はある程度の細かさで十分である。

　論文を読むとレプリカの数は20あるいは30程度とるのが良いとある。また、 $PT$ はコンスタントになるように計算することが多いようである（これは必然性はない）。私はレプリカの数を10程度としてSAで実装した。それなりの結果が得られたようである。一つのアイデアとして遺伝的アルゴリズム（GA）など他の手法を組み合わせるのは面白いのではないかと考える。

　当然のことであるが、レプリカの数だけ一回のモンテカルロステップでQAのほうが時間がかかる。それを考慮してモンテカルロステップを減らしてもQAの方がSAよりいい結果が出ることがあり、実時間で短くできるはずということである。ただ、論文ではモンテカルロステップについて述べているものが多く、実時間で比較しないで意味があるのだろうかという疑問を感じた。

　SAとQAの優劣についてである。問題によっては良くなる場合もあるし、悪くなる場合もあると論文には書かれている。ただ残念ながら私はどのようなケースで良くなり、悪くなるかということに関して知見はない。

　このアルゴリズムで重要な概念として、鈴木トロッター展開、経路積分展開、分配関数があり、量子力学と統計力学の知識が多く盛り込まれている。これについては物理の参考書あるいは論文を参照していただきたい。

　QAと多スタート局所探索法など他の最適化手法との比較は面白いと感じている。古典的手法と量子的手法の比較ができると考えるからである。実際、鈴木トロッター展開をした後は古典的手法にマップされているので、他の手法と簡単に比較できる。

　最後に簡単なコードを書いてみた。下記にサンプルコードをあげた。稚拙なコードであるがご了承願いたい。コードは書き直す予定である。

SAのコード：https://github.com/hmiyahara512/simulated-annealing_for_3d-ising-model

QAのコード：https://github.com/hmiyahara512/quantum-annealing_for_3d-ising-model

[1] T. Kadowaki and H. Nishimori, "Quantum annealing in the transverse Ising model"

Phys. Rev. E 58, 5355 (1998)

[2] A. Das, B.K. Chakrabarti, and R.B. Stinchcombe, "Quantum annealing in a kinetically constrained system", Phys. Rev. E 72 art. 026701 (2005)

[3] G. E. Santoro, R. Martonak, E. Tosatti, and R. Car, "Theory of Quantum Annealing of an Ising Spin Glass" DOI: 10.1126/science.1068774

[4] Edward Farhi, Jeffrey Goldstone, Sam Gutmann, "Quantum Adiabatic Evolution Algorithms versus Simulated Annealing" http://arxiv.org/abs/quant-ph/0201031

[5] Edward Farhi, Jeffrey Goldstone, Sam Gutmann, Joshua Lapan, Andrew Lundgren, Daniel Preda, "A Quantum Adiabatic Evolution Algorithm Applied to Random Instances of an NP-Complete Problem" DOI: 10.1126/science.1057726

[6] Edward Farhi, Jeffrey Goldstone, Sam Gutmann, "A Numerical Study of the Performance of a Quantum Adiabatic Evolution Algorithm for Satisfiability" http://arxiv.org/abs/quant-ph/0007071

[7] Roman Martonak (1), Giuseppe E. Santoro (2,3), Erio Tosatti (2,3) ((1) Dept. of Chemistry and Applied Biosciences, ETH Zuerich, Lugano, CH, (2) SISSA and INFM-Democritos, Trieste, Italy, (3) ICTP, Trieste, Italy) "Quantum annealing of the Traveling Salesman Problem" DOI: 10.1103/PhysRevE.70.057701

[8] S. Morita and H. Nishimori, "Mathematical foundation of quantum annealing", J.Math. Phys. 49, 125210 (2008)

[9] http://www-adsys.sys.i.kyoto-u.ac.jp/mohzeki/QA.pdf

[10] M. Suzuki, "Relationship between d-Dimensional Quantal Spin Systems and (d+1)-Dimensional Ising System : Equivalence, Critical Exponents and Systematic Approximants of the Partition Function and Spin Correlations" http://ci.nii.ac.jp/lognavi?name=nels&lang=en&type=pdf&id=ART0001548127

[11] Arnab Das and Bikas K. Chakrabarti
"Quantum annealing and analog quantum computation" http://journals.aps.org/rmp/abstract/10.1103/RevModPhys.80.1061

[12] V.N. Smelyanskiy, E.G. Rieffel, S.I. Knysh, C.P. Williams, M.W. Johnson, M.C. Thom, and K.L.P.W. G. Macready (2012), “A Near-Term Quantum Computing Approach for Hard Computational Problems in Space Exploration”

＊解釈あるいはコードに間違いがあれば、ご一報をお願いします。

＊8月12日に内容と出典の誤りを修正しました。随時修正します。一発目から長いエントリを書いたことを後悔してます。

＊モンテカルロステップやエネルギーの計算に不適切な部分があります。時間があれば修正するかもしれません。（許容範囲かなと思います。）