経路積分量子モンテカルロ法を用いた量子アニーリングのシミュレータ

0または1の状態しかとれないビットを用いた古典的な計算機において量子力学の原理に基づいた量子アニーリングをシミュレートできるのだろうか．
量子アニーリングの理論について学んだ（前回の記事）のちに，このような疑問が生まれた．調べてみると，どうやら手元のPCでシミュレータを作るには経路積分量子モンテカルロ法という手法がよく用いられるようだった．そのため今回は，経路積分量子モンテカルロ法について学び，Pythonでの既存実装に触れてみることで，量子アニーリングについての理解をより深めていくことを目的とする．

経路積分量子モンテカルロ法
シミュレーテッドアニーリングとの対比
計算の流れ
シミュレータを触ってみる
まとめ
参考

経路積分量子モンテカルロ法

今回，古典的な計算機を用いて量子アニーリングをシミュレートするために，経路積分量子モンテカルロ法という手法を用いる．量子系をシミュレートするために，まず量子ゆらぎを含む横磁場イジングモデルの分配関数 $Z$ を，対応する古典イジングモデルの分配関数で表す．そのとき，これから示す計算を追っていくと，横磁場イジングモデルの分配関数が1次元高い古典の分配関数に対応することがわかる．この結果を使って量子系の平衡状態における性質を古典系のマルコフ連鎖モンテカルロ法で調べる，というのが量子モンテカルロの基本的な考え方である．

以下，理論的な背景を少し詳細に見ていく．
スピン数が $N$ の場合，横磁場イジングモデルのハミルトニアンは以下のように表される．

$\displaystyle \begin{eqnarray} \hat{H} & = & \hat{H_0}+\hat{H_q} \\ \hat{H_0} & = & -\sum_{i,j=1}^NJ_{ij}\hat{\sigma}_i^z\hat{\sigma}_j^z-\sum_{i=1}^Nh_i\hat{\sigma}_i^z \\ \hat{H_q} & = & -\Gamma\sum_{i=1}^N\hat{\sigma}_i^x \\ \end{eqnarray}$

ここで， $\hat{\sigma_i^z}, \hat{\sigma_i^x}$ はパウリ行列の $z，x$ 成分であり，以下の行列で表される．

$\displaystyle \hat{\sigma_i^z} = \begin{pmatrix} 1 & 0 \\ 0 & -1 \\ \end{pmatrix} , \hat{\sigma_i^x} = \begin{pmatrix} 0 & 1 \\ 1 & 0 \\ \end{pmatrix}$

このハミルトニアンを用いて分配関数は以下のように表せる．

$\displaystyle \begin{eqnarray} Z & = & Tr[e^{-\beta\hat{H}}] \\ & = & Tr[e^{-\beta(\hat{H_0}+\hat{H_q})}] \\ & = & \sum_{\boldsymbol{\sigma_1}}\langle\boldsymbol{\sigma_1}\mid{e^{-\beta(\hat{H_0}+\hat{H_q})}}\mid\boldsymbol{\sigma_1}\rangle \\ \end{eqnarray}$

ここで， $\beta$ は温度の逆数（逆温度）である．
Suzuki・Trotter分解により，

$\displaystyle e^{-\beta(\hat{H_0}+\hat{H_q})} = \lim_{M\rightarrow\infty}\biggl(e^{-\frac{\beta}{M}\hat{H_0}}e^{-\frac{\beta}{M}\hat{H_q}}\biggl)^M$

と表せることを用いると，

$\displaystyle Z = \lim_{M\rightarrow\infty}\sum_{\boldsymbol{\sigma_1}}\langle\boldsymbol{\sigma_1}\mid\biggl(e^{-\frac{\beta}{M}\hat{H_0}}e^{-\frac{\beta}{M}\hat{H_q}}\biggl)^M\mid\boldsymbol{\sigma_1}\rangle$

ここでM個の積の間に以下の恒等演算子を挟み込むと，

$\displaystyle \boldsymbol{1} = \sum_{\boldsymbol{\sigma_k}}\mid\sigma_k\rangle\langle\sigma_k\mid (k = 1, 2, ..., M)$

分配関数は以下のように表される．

$\displaystyle \begin{eqnarray} Z & = & \lim_{M\rightarrow\infty}\sum_{\boldsymbol{\sigma_k},\boldsymbol{\sigma_k'}}\langle\boldsymbol{\sigma_1}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_1'}\rangle\langle\boldsymbol{\sigma_1'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_2}\rangle \\ & \times & \langle\boldsymbol{\sigma_2}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_2'}\rangle\langle\boldsymbol{\sigma_2'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_3}\rangle \\ & \times & ・・・ \\ & \times & \langle\boldsymbol{\sigma_M}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_M'}\rangle\langle\boldsymbol{\sigma_M'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_1}\rangle \end{eqnarray}$

ここで， $\mid\boldsymbol{\sigma_k}\rangle$ は $N$ 個のスピン系の直積空間を表す．

$\displaystyle \mid\boldsymbol{\sigma_k}\rangle=\mid\boldsymbol{\sigma^z_{1,k}}\rangle\otimes\mid\boldsymbol{\sigma^z_{2,k}}\rangle\otimes・・・\otimes\mid\boldsymbol{\sigma^z_{N,k}}\rangle$

分配関数の式から， $k = 1, 2, ..., M$ は元の問題の持つ変数が定義された空間とは別の独立な変数のインデックスであり，トロッタ方向または虚時間方向と呼ばれる．
下図はSuzuki・Trotter分解を行なった際にトロッタ方向に次元が追加されている様子を示した概略図である．

f:id:hirotsuru314:20191210223523p:plain:w500 — Suzuki-Trotter分解の概略図

これまでの結果から，分配関数を求めるには，以下の2つの値を求めればよいことがわかる．

$\displaystyle \begin{eqnarray} \langle\boldsymbol{\sigma_k}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_k'}\rangle \tag{1} \\ \langle\boldsymbol{\sigma_k'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_{k+1}}\rangle \tag{2} \end{eqnarray}$

$(1)$ は $\hat{H_0}$ が対角行列であることから，

$\displaystyle \begin{eqnarray} \hat{\sigma}_i^z\mid\sigma_k\rangle & = & \hat{\sigma}_{i, k}\mid\sigma_k\rangle \\ \langle\boldsymbol{\sigma_k}\mid e^{-\frac{\beta}{M}\hat{H_0}}\mid\boldsymbol{\sigma_k'}\rangle & = & exp\biggl[\frac{\beta}{M}\biggl(\sum_{i,j=1}^NJ_{ij}\hat{\sigma}_{i, k}\hat{\sigma}_{j, k}+\sum_{i=1}^Nh_i\hat{\sigma}_{i, k}\biggl)\biggl]\prod_{i=1}^N\delta(\hat{\sigma}_{i, k}, \hat{\sigma}_{i, k}') \end{eqnarray}$

となる．次に $(2)$ は以下の式が成立する．

$\displaystyle \begin{eqnarray} \langle\boldsymbol{\sigma_k'}\mid e^{-\frac{\beta}{M}\hat{H_q}}\mid\boldsymbol{\sigma_{k+1}}\rangle = \biggl[\frac{1}{2}\sinh\biggl(\frac{2\beta\gamma}{M}\biggl)\biggl]^{\frac{N}{2}}exp\biggl[\frac{1}{2}\log \coth\biggl(\frac{\beta\Gamma}{M}\sum_{i=1}^N\hat{\sigma}_{i, k}'\hat{\sigma}_{i, k+1}\biggl)\biggl] \\ \end{eqnarray}$

よって，分配関数 $Z$ は以下のように表すことができる．

$\displaystyle Z = \lim_{M\rightarrow\infty}\biggl[\frac{1}{2}\sinh\biggl(\frac{2\beta\Gamma}{M}\biggl)\biggl]^{\frac{NM}{2}}\sum_{\{\hat{\sigma}_{i, k}=\pm1\}}exp\biggl\{\sum_{k=1}^M\biggl[\sum_{i, j=1}^N\frac{\beta J_{ij}}{M}\hat{\sigma}_{i, k}\hat{\sigma}_{j, k}-\sum_{i=1}^N\frac{\beta h_{i}}{M}\hat{\sigma}_{i, k} \\ +\frac{1}{2}\ln \coth\biggl(\frac{\beta\Gamma}{M}\sum_{i=1}^N\hat{\sigma}_{i, k}\hat{\sigma}_{i, k+1}\biggl)\biggl] \\$

この結果から横磁場イジングモデルの分配関数は対応する磁場なしイジングモデルに1次元トロッタ方向に次元を追加したもので表されることがわかる．対応するハミルトニアンは以下の式で表される．

$\displaystyle H_{eff} = -\sum_{i, j=1}^N\sum_{k=1}^M\frac{J_{ij}}{M}\hat{\sigma}_{i, k}\hat{\sigma}_{j, k}-\sum_{i=1}^N\sum_{k=1}^M\frac{h_{i}}{M}\hat{\sigma}_{i, k}-\frac{1}{2\beta}\ln \coth\biggl(\frac{\beta\Gamma}{M}\biggl)\sum_{i=1}^N\sum_{k=1}^M\hat{\sigma}_{i, k}'\hat{\sigma}_{i, k+1} \\$

シミュレーテッドアニーリングとの対比

これまでの計算から，横磁場イジングモデルの分配関数は対応する古典イジングモデルの分配関数に対して虚時間方向に次元を一つ追加したもので表せることがわかった．このことから，量子アニーリングのシミュレータは，一見単にトロッタ数 $M$ で並列化したシミュレーテッドアニーリングと同じもののように感じるかもしれないが，上記で計算したハミルトニアンを読み解いてみると，そうでないことがわかる．
ハミルトニアン $H_{eff}$ の第1・2項を見ると，同じ層内（ $k$ が同一）での相互作用のみであるが，横磁場の効果を表す第三項では，層間の相互作用を持っている．（ $k$ と $k+1$ の相互作用）このことから時間発展とともに層同士が互いに独立ではなく，相互に作用し合うため，単にトロッタ数分の並列処理をしているわけではないようだ．
どのような相互作用でそれがどう基底状態探索に効いているのかなどについてはまた別途追っていきたい．

計算の流れ

計算の全体的な流れは以下のようになる．

パラメータの初期化
・温度 $T$
・横磁場の強さ $\Gamma$
スピンの初期化 $M$ 層のそれぞれ $N$ 個のスピンをランダムに初期化する．
ランダムに層と層中のスピンを選択し，スピンの向きを反転させる．
スピン反転前後のエネルギーの差分を $\Delta E$ とすると，確率 $min(1, e^\frac{\Delta E}{T})$ でスピンの向きの反転を受け入れる（受け入れない場合は棄却して変更なし）
任意の回数3，4のステップを繰り返す（モンテカルロステップ）
横磁場の強さを $\Gamma$ 小さくする．
任意の回数5，6のステップを繰り返す（アニーリングステップ）
各層のエネルギーを求め（横磁場の項は除く）一番エネルギーが低い層のスピン配位を解とする

シミュレータを触ってみる

GitHub上ではいくつかの既存実装が存在したが，今回は以下のものを触ってみることにした．

github.com

必要なライブラリ環境をそろえて，リポジトリのトップディレクトリで以下のように実行すると，ヘルプが見れる．ちなみにPython 2系で実装されている．

$ python solve.py qmc -h
usage: Ising Solvers qmc [-h] [-problem PROBLEM] [-steps STEPS] [-dump DUMP]
                         [-spins SPINS] [-P P] [-T T] [-G0 G0] [-Gf GF]
                         [-e0 E0] [-ef EF]

optional arguments:
  -h, --help        show this help message and exit
  -problem PROBLEM  a file containing problem data in three column format, see
                    README for details
  -steps STEPS      number of monte carlo steps (MCS) used, good values depend
                    on problem size
  -dump DUMP        a file to dump move acceptances to, if desired -- slows
                    down performance
  -spins SPINS      initial spin configuration, if desired
  -P P              number of trotter slices (replicas), 20 to 80 is typical
  -T T              ambient temperature in K, 0.015 is typical
  -G0 G0            initial transverse field strength
  -Gf GF            final tranverse field strength, typically near zero
  -e0 E0            initial coupling pre-factor
  -ef EF            final coupling pre-factor

多くの引数があるが，基本的には指定しなくてもデフォルト値が入るようになっていて，必要であれば実行者が上書きできるようにしてある．
マストで渡さないといけない引数は-problem であり，これは，以下のような形式のファイルのパスを渡す必要がある．（sample_dataディレクトリ配下にサンプルが置いてある）

$ head -5 sample_data/ising12.txt
random 12x12 spin ising model with solution energy of -18,972,276
1 2 -28420
2 3 17885
3 4 -146311
4 5 128784

2行以降は以下の形式で

i j J_ij

すなわち，スピン間の相互作用を示す $J_{ij}$ を設定するためのファイルであることがわかる．今回の最適化では実用的な問題が設定してあるのではなく，単にイジングモデルのスピン間にファイルで設定した相互作用を与え，そのイジングモデルの基底状態を探索するのみである．また， -dump 引数を指定しておくと，実行の過程がファイル出力される．-steps でアニーリングステップ数が指定できる．
以下のように実行すると，用いたパラメータ等の出力が返ってくる．

$ python solve.py qmc -problem sample_data/ising12.txt -dump output.txt -steps 1000

Ising Spin Glass
+----------------+-------------------------------------------------------------------+
| Field          | Value                                                             |
+----------------+-------------------------------------------------------------------+
| data file      | sample_data/ising12.txt                                           |
| description    | random 12x12 spin ising model with solution energy of -18,972,276 |
| initial config | 0xe696019ff2add3594c219fcbfd9aef9b0dd4                            |
| current config | 0x313fac59f68f16d3ec6a17d6aca964b2025d                            |
| initial energy | 1627456.0000000007                                                |
| current energy | -15899168.000000006                                               |
+----------------+-------------------------------------------------------------------+

Solver: Path-Integral Quantum Monte Carlo
+-----------+------------+
| Parameter | Value      |
+-----------+------------+
| G0        | 3          |
| Gf        | 1e-06      |
| P         | 40         |
| T         | 0.015      |
| dump      | output.txt |
| e0        | 1e-06      |
| ef        | 4          |
| steps     | 1000       |
+-----------+------------+
To replicate starting conditions, run with arguments:
    qmc -dump output.txt -G0 3 -ef 4 -P 40 -steps 1000 -T 0.015 -problem sample_data/ising12.txt -Gf 1e-06 -e0 1e-06 -spins 0xe696019ff2add3594c219fcbfd9aef9b0dd4

出力ファイルを見てみると以下のような形式で出力されている．

$ head output.txt
Path-Integral Quantum Monte Carlo
{'dump': 'output.txt', 'G0': 3, 'ef': 4, 'P': 40, 'steps': 1000, 'T': 0.015, 'Gf': 1e-06, 'e0': 1e-06}
<class 'ising.SpinGlass'>(data_file="sample_data/ising12.txt", spin_configuration="0xe696019ff2add3594c219fcbfd9aef9b0dd4")
3,16.27456,0xe696019ff2add3594c219fcbfd9aef9b0dd4,-1
3.0,15.64962,0xe696019ff2add3594c219fcbfd9aef9b0cd4,7
3.0,13.09112,0xe696019ff2add3594c219fcbfd9aef9b0ddc,28
3.0,13.02468,0xe69e019ff2add3594c219fcbfd9aef9b0dd4,12
3.0,16.31304,0xe696019ff2add3194c219fcbfd9aef9b0dd4,30
3.0,15.9128,0xe696019ff2add3594c219fcbfd9aefbb0dd4,0
3.0,14.19152,0xe697019ff2add3594c219fcbfd9aef9b0dd4,24

4行目以降の結果は，左から横磁場の強さ，全層の中から最も低いエネルギーの値，スピンの状態を16進数で表したもの，何層目のスピンを反転させたのか．
このライブラリの実装では，1アニーリングステップの中で，それぞれの層の中から1スピンずつ反転させている．すなわちモンテカルロステップ数はトロッタ数と同じである．そのモンテカルロステップごとに1行の出力がある．

ステップを重ねるにつれてイジングモデルが低エネルギー状態に向かっていることを確認するために，横軸にアニーリングステップ，縦軸はエネルギーでグラフを書くと以下のようになる．

f:id:hirotsuru314:20191216173856p:plain

アニーリングステップが進むにつれてエネルギーが低い状態に向かっていることから，基底状態を探索している様子が見て取れる．

まとめ

古典的な計算機で量子アニーリングをシミュレートする代表的な手法として，経路積分量子モンテカルロ法がある．経路積分を用いる過程でSuzuki・Trotter展開によって得られるハミルトニアンの分配関数の近似値から，マルコフ連鎖モンテカルロ法を用いることで，量子アニーリングの挙動をシミュレートすることができる．動作としては，問題とするイジングモデルをトロッタ数分だけコピーして，それらが互いに相互作用しながら基底状態を探索していく．今回用いたqmcのような既存実装を用いると，実用的な問題ではないものの手元のPCで簡単に量子アニーリングの挙動をシミュレートすることができる．

Fire Engine

消防士→ITエンジニア→研究者

経路積分量子モンテカルロ法を用いた量子アニーリングのシミュレータ

経路積分量子モンテカルロ法

シミュレーテッドアニーリングとの対比

計算の流れ

シミュレータを触ってみる

まとめ

参考