Simulated Quantum Annealingを導出する

こんにちは、株式会社Jijの山城です。

この記事はJij Inc. Advent Calendar 2023の7日目の記事です。
前回はJijの広報による "ChatGPTを使って広報活動レポートを作成し業務効率化" でした。

Jij(ジェイアイジェイ)は数理最適化、量子技術を用いたソフトウェアの開発を行っています。
Jijの社名の由来となっているのは統計物理のイジングモデルに関連して今回はシミュレーテッド量子アニーリングに関する記事を書いていきます。

はじめに

想定読者

物理学科の学部生くらいで、基礎的な量子力学、統計力学を学んでいることを想定します。簡単なブラケット演算や分配関数などは既知とします。

量子アニーリング

量子アニーリングを模したアルゴリズムとしてシミュレーテッド量子アニーリングというアルゴリズムがあります。量子アニーリングを通常のコンピュータでシミュレーションする方法の一つです。量子アニーリングは量子力学つまり、シュレディンガー方程式に従って実行されるのですが、シュレディンガー方程式を愚直にシミュレートするのは通常のコンピュータでは困難です。なのでシュレディンガー方程式を使わず平衡統計力学の手法を使ってシミュレートします。

しかし平衡統計力学とシュレディンガー方程式によるダイナミクスはそもそも異なる物理です。しかしなぜ平衡統計力学を持ち出すのでしょうか？大雑把には量子アニーリングにおいて断熱時間発展を考えて時間発展中もハミルトニアンの基底状態付近にいるとすると、低温の平衡状態（基底状態）で近似してもよさげな結果になるだろうという考えから来ています。

この記事では横磁場のみの量子アニーリングに対するシミュレーテッド量子アニーリング (Simulated Quantum Annealing; SQA) を紹介します。特にアルゴリズムを導出するのに必要なトロッター展開を用いた量子モンテカルロ法を中心とした記事です。

SQAを実際に実行してみたい方はOpenJijをぜひ利用してみてください。

SQAの導出

Hamiltonian of the transverse field quantum annealing

横磁場量子アニーリングのハミルトニアンは以下で記述されます。

$$H^(t)=−A(t)∑iσ^ix+B(t)(∑iσ^iz+∑i,jJijσ^izσ^jz)=−A(t)∑iσ^ix+B(t)Hp(σ^z)\hat H (t) = -A(t) \sum_i \hat \sigma_i^x + B(t) \left(\sum_i \hat \sigma_i^z + \sum_{i,j}J_{ij}\hat \sigma^z_i \hat \sigma^z_j\right) = -A(t) \sum_i \hat \sigma_i^x + B(t)H_p(\hat \sigma^z)$$

where

$$∑iσ^iz+∑i,jJijσ^izσ^jz=Hp(σ^z)\sum_i \hat \sigma_i^z + \sum_{i,j}J_{ij}\hat \sigma^z_i \hat \sigma^z_j = H_p (\hat \sigma^z)$$

このハミルトニアン下での分配関数を考えましょう。

Partition function

平衡統計力学でのこのハミルトニアン下での状態を考えるわけなので分配関数を考えましょう。

$$Z(β)=Trexp⁡(−βH^(t))=Trexp⁡(βA(t)∑iσ^ix−βB(t)Hp(σ^z))Z(\beta) = \mathrm{Tr}\exp\left( -\beta \hat H(t) \right) = \mathrm{Tr} \exp\left(\beta A(t) \sum_i \hat \sigma_i^x - \beta B(t)H_p(\hat \sigma^z)\right)$$

実は以下の計算の手続きでこの$Tr\mathrm{Tr}$ を処理して古典的なスピン配位の和に変換することで、量子系の分配関数を古典系の分配関数へ書き換えることができます。

そうすることで古典系の分配関数になってしまえばあとはそれをそのままモンテカルロ法すればよいわけです。シミュレーテッドアニーリングがモンテカルロ法をベースであることを思い出すと、この系に対してのモンテカルロ法をそのままアニーリングのようにパラメータを変化させながら実行すれば量子アニーリングっぽいアルゴリズムが作れることが想像できると思います。では早速この分配関数に対しての計算を進めていきましょう。

Suzuki-Trotter decomposition

そのままでは指数の中身の2つのハミルトニアンが非可換なので計算を進めることができません。ここで鈴木トロッター展開と呼ばれる方法を使って指数関数を2つに分解します。

$$Z=Trexp⁡(βA(t)∑iσ^ix−βB(t)Hp(σ^z))=Tr[exp⁡(βmA(t)∑iσ^ix)exp⁡(−βmB(t)Hp(σ^z))]m+O((βm)2)Z =\mathrm{Tr} \exp\left(\beta A(t) \sum_i \hat \sigma_i^x -\beta B(t)H_p(\hat \sigma^z)\right)\\ = \mathrm{Tr} \left[\exp\left(\frac{\beta}{m}A(t) \sum_i \hat \sigma_i^x\right)\exp\left(-\frac{\beta}{m}B(t)H_p(\hat \sigma^z)\right)\right]^m + \mathcal O\left(\left(\frac{\beta}{m}\right)^2\right)$$

あるスピン配位を $σ⃗\vec\sigma$ と書くことにします。トレースはこのスピン配位の全パターンでの状態和なので、

$$Z≃∑σ⃗⟨σ⃗∣[exp⁡(βmA(t)∑iσ^ix)exp⁡(−βmB(t)Hp(σ^z))]m∣σ⃗⟩Z \simeq \sum_{\vec \sigma} \bra{\vec \sigma} \left[\exp\left(\frac{\beta}{m}A(t) \sum_i \hat \sigma_i^x\right)\exp\left(-\frac{\beta}{m}B(t)H_p(\hat \sigma^z)\right)\right]^m \ket{\vec \sigma}$$

です。中の演算は同じ形の行列の指数関数が$mm$個ありますが、その積の間に完全系の式

$$1^=1^k=∑σ⃗∣σ⃗k⟩⟨σ⃗k∣\hat 1 = \hat 1_k = \sum_{\vec \sigma} | \vec \sigma_k\rang \lang \vec \sigma_k|$$

を挟みます。ここで添字$kk$は$1^\hat 1$をたくさん入れていくので区別がつくように入れました。分配関数に$1^\hat 1$を挟んでいくと

$$Z≃∑σ⃗⟨σ⃗∣1^1[exp⁡(βmA(t)∑iσ^ix)exp⁡(−βmB(t)Hp(σ^z))]1^2[exp⁡(βmA(t)∑iσ^ix)exp⁡(−βmB(t)Hp(σ^z))]1^3⋯ ∣σ⃗⟩ Z \simeq \sum_{\vec \sigma} \left \lang\vec \sigma\left| \hat 1_1 \left[\exp\left(\frac{\beta}{m}A(t) \sum_i \hat \sigma_i^x\right) \exp\left(-\frac{\beta}{m}B(t)H_p(\hat \sigma^z)\right)\right] \hat 1_2 \left[\exp\left(\frac{\beta}{m}A(t) \sum_i \hat \sigma_i^x\right)\exp\left(-\frac{\beta}{m}B(t)H_p(\hat \sigma^z)\right)\right]\hat 1_3\cdots \right| \vec \sigma\right\rang$$

$$=∑σ⃗⟨σ⃗∣∑σ⃗1∣σ⃗1⟩⟨σ⃗1∣[exp⁡(βmA(t)∑iσ^ix)exp⁡(−βmB(t)Hp(σ^z))]∑σ⃗2∣σ⃗2⟩⟨σ⃗2∣⋯ ∣σ⃗⟩=∏k=1m∑σ⃗k⟨σ⃗k∣[exp⁡(βmA(t)∑iσ^ix)exp⁡(−βmB(t)Hp(σ^z))]∣σ⃗k+1⟩ = \sum_{\vec \sigma} \left \lang \vec \sigma\left| \sum_{\vec \sigma_1}| \vec \sigma_1\rang \lang \vec \sigma_1| \left[\exp\left(\frac{\beta}{m}A(t) \sum_i \hat \sigma_i^x\right) \exp\left(-\frac{\beta}{m}B(t)H_p(\hat \sigma^z)\right)\right] \sum_{\vec \sigma_2}| \vec \sigma_2\rang \lang \vec \sigma_2|\cdots \right| \vec \sigma\right\rang\\ =\prod_{k=1}^m \sum_{\vec \sigma_k}\left \lang \vec \sigma_k\left| \left[\exp\left(\frac{\beta}{m}A(t) \sum_i \hat \sigma_i^x\right) \exp\left(-\frac{\beta}{m}B(t)H_p(\hat \sigma^z)\right)\right] \right| \vec \sigma_{k+1}\right\rang$$

となります。ここで $m+1=1m+1 = 1$ とします。また

$$exp⁡(−βmB(t)Hp(σ^z))∣σ⃗k⟩=exp⁡(−βmB(t)Hp(σ⃗k))∣σ⃗k⟩ \exp\left(-\frac{\beta}{m}B(t)H_p(\hat \sigma^z)\right) |\vec \sigma_k\rang = \exp\left(-\frac{\beta}{m}B(t)H_p(\vec \sigma_k)\right) |\vec \sigma_k\rang$$

なので$HpH_p$に関するexpは行列からただの数に変わります。よってそのまま外にくくりだします。次に

$$⟨σ⃗k∣exp⁡(βmA(t)∑iσ^ix)∣σ⃗k+1⟩ \lang\vec \sigma_k|\exp\left(\frac{\beta}{m}A(t) \sum_i \hat \sigma_i^x\right) |\vec \sigma_{k+1}\rang$$

を考えます。あとはこの横磁場項の一体問題を解けばよいです。スピン配位を1スピンごとのテンソル積 $∣σ⃗k⟩=⊗i∣(σi)k⟩\ket{\vec \sigma_k} = \otimes_i \ket{(\sigma_i)_k}$ で記述すると

$$⟨σ⃗k∣exp⁡(βmA(t)∑iσ^ix)∣σ⃗k+1⟩=∏i⟨(σi)k∣exp⁡(βmA(t)σ^ix)∣(σi)k+1⟩ \lang\vec \sigma_k|\exp\left(\frac{\beta}{m}A(t) \sum_i \hat \sigma_i^x\right) | \vec \sigma_{k+1}\rang = \prod_{i} \lang(\sigma_i)_k|\exp\left(\frac{\beta}{m}A(t) \hat \sigma_i^x\right) |(\sigma_i)_{k+1}\rang$$

となります。ではこの一体問題を考えましょう。行列の指数関数はテイラー展開で定義されるのでそのままexpを展開してやって、$(σ^ix)2=1(\hat \sigma_i^x)^2 = 1$であることを使うと、

$$exp⁡(βmA(t)σ^ix)=1^cosh⁡(βA(t)m)+σ^ixsinh⁡(βA(t)m) \exp\left(\frac{\beta}{m}A(t) \hat \sigma_i^x\right) = \hat 1 \cosh\left(\frac{\beta A(t)}{m}\right) + \hat \sigma_i^x \sinh\left(\frac{\beta A(t)}{m}\right)$$

となります。よって

$$⟨(σi)k∣exp⁡(βmA(t)σ^ix)∣(σi)k+1⟩=⟨(σi)k∣(1^cosh⁡(βA(t)m)+σ^ixsinh⁡(βA(t)m))∣(σi)k+1⟩=cosh⁡(βA(t)m)⟨(σi)k∣(σi)k+1⟩+sinh⁡(βA(t)m)⟨(σi)k∣−(σi)k+1⟩ \lang(\sigma_i)_k|\exp\left(\frac{\beta}{m}A(t) \hat \sigma_i^x\right) |(\sigma_i)_{k+1}\rang = \lang(\sigma_i)_k|\left(\hat 1 \cosh\left(\frac{\beta A(t)}{m}\right) + \hat \sigma_i^x \sinh\left(\frac{\beta A(t)}{m}\right)\right) |(\sigma_i)_{k+1}\rang \\ = \cosh\left(\frac{\beta A(t)}{m}\right) \lang(\sigma_i)_k|(\sigma_i)_{k+1}\rang + \sinh\left(\frac{\beta A(t)}{m}\right) \lang(\sigma_i)_k|-(\sigma_i)_{k+1}\rang$$

となります。ここで$σ^ix\hat \sigma_i^x$はスピンを反転させるので$σ^ix∣σi⟩=∣−σi⟩\hat \sigma_i^x|\sigma_i\rang = |-\sigma_i\rang$を使いました。

ここまでで分配関数は

$$Z≃∏k∑σ⃗kexp⁡(−βmB(t)Hp(σ⃗k))(cosh⁡(βA(t)m)⟨(σi)k∣(σi)k+1⟩+sinh⁡(βA(t)m)⟨(σi)k∣−(σi)k+1⟩) Z \simeq \prod_k\sum_{\vec \sigma_k}\exp\left(-\frac{\beta}{m}B(t)H_p(\vec \sigma_k)\right) \left(\cosh\left(\frac{\beta A(t)}{m}\right) \lang(\sigma_i)_k|(\sigma_i)_{k+1}\rang + \sinh\left(\frac{\beta A(t)}{m}\right) \lang(\sigma_i)_k|-(\sigma_i)_{k+1}\rang\right)$$

となりました。$σ^x\hat \sigma^x$ からの寄与も指数の肩に載せて整理します。

ここで、

$$cosh⁡(βA(t)m)=cosh⁡(βA(t)m)sinh⁡(βA(t)m)(tanh⁡(βA(t)m))−1=12sinh⁡(2βA(t)m)(tanh⁡(βA(t)m))−1sinh⁡(βA(t)m)=12sinh⁡(2βA(t)m)tanh⁡(βA(t)m) \cosh\left(\frac{\beta A(t)}{m}\right) = \sqrt{\cosh\left(\frac{\beta A(t)}{m}\right) \sinh\left(\frac{\beta A(t)}{m}\right)\left(\tanh\left(\frac{\beta A(t)}{m}\right)\right)^{-1}} =\sqrt{\frac{1}{2}\sinh\left(2\frac{\beta A(t)}{m}\right)\left(\tanh\left(\frac{\beta A(t)}{m}\right)\right)^{-1}}\\ \sinh\left(\frac{\beta A(t)}{m}\right) = \sqrt{\frac{1}{2}\sinh\left(2\frac{\beta A(t)}{m}\right)\tanh\left(\frac{\beta A(t)}{m}\right)}$$

であることを使います。この時 $cosh⁡\cosh$ と $sinh⁡\sinh$ の違いは $tanh⁡\tanh$ のべきだけなので以下のように $σ^x\hat \sigma^x$ からの寄与を計算することができます。整理すると

$$cosh⁡(βA(t)m)⟨(σi)k∣(σi)k+1⟩+sinh⁡(βA(t)m)⟨(σi)k∣−(σi)k+1⟩=12sinh⁡(2βA(t)m)tanh⁡(βA(t)m)−(σi)k(σi)k+1=12sinh⁡(2βA(t)m)exp⁡(−12ln⁡tanh⁡(βA(t)m)(σi)k(σi)k+1) \cosh\left(\frac{\beta A(t)}{m}\right)\lang(\sigma_i)_k|(\sigma_i)_{k+1}\rang + \sinh\left(\frac{\beta A(t)}{m}\right) \lang(\sigma_i)_k|-(\sigma_i)_{k+1}\rang\\ = \sqrt{\frac{1}{2}\sinh\left(2\frac{\beta A(t)}{m}\right)\tanh\left(\frac{\beta A(t)}{m}\right)^{-(\sigma_i)_k (\sigma_i)_{k+1}}}\\ =\sqrt{\frac{1}{2}\sinh\left(2\frac{\beta A(t)}{m}\right)}\exp\left(-\frac{1}{2}\ln\tanh\left(\frac{\beta A(t)}{m}\right){(\sigma_i)_k (\sigma_i)_{k+1}}\right)$$

となります。よって分配関数は

$$Z≃∏k∑σ⃗kexp⁡(−βm(B(t)Hp(σ⃗k)+m2βln⁡tanh⁡(βA(t)m)(σi)k(σi)k+1))=∑σ⃗1,⋯ ,σ⃗mexp⁡(−βm(B(t)∑kHp(σ⃗k)+m2βln⁡tanh⁡(βA(t)m)∑k(σi)k(σi)k+1)) Z \simeq \prod_k\sum_{\vec \sigma_k} \exp\left(-\frac{\beta}{m}\left( B(t)H_p(\vec \sigma_k) +\frac{m}{2\beta}\ln\tanh\left(\frac{\beta A(t)}{m}\right){(\sigma_i)_k (\sigma_i)_{k+1}} \right)\right)\\ =\sum_{\vec \sigma_1,\cdots,\vec \sigma_m} \exp\left(-\frac{\beta}{m}\left( B(t)\sum_k H_p(\vec \sigma_k) +\frac{m}{2\beta}\ln\tanh\left(\frac{\beta A(t)}{m}\right)\sum_k{(\sigma_i)_k (\sigma_i)_{k+1}} \right)\right)$$

となります。スピン配位にかかわらない係数は省略しました。
ここで量子項があるハミルトニアンの分配関数から古典スピン系に関する分配関数へと書き換えることができました。
この古典スピン系は元の系からスピン変数が$mm$倍に増えて、添え字が追加され、$k,k+1k, k+1$ 間に相互作用が追加された系となっています。
イメージとしては量子の重ね合わせ状態を複数のレプリカでシミュレートしているというイメージを持っても良いかもしれません。

$βA(t)m≥0\beta A(t) m \geq 0$なので、$ln⁡tanh⁡(βA(t)/m)≤0\ln \tanh\left(\beta A(t) / m\right) \leq 0$ となります。つまり量子項によって導入された項はレプリカ間の強磁性相互作用（スピン同士が同じ向きになるとエネルギーが低くなる）となっています。

今回はこのSQAの導出までとし、今度さらにこの古典系のイメージについての説明や実際にOpenJijでSQAを実行する方法について書きたいと思います。SQAでは$A(t)A(t)$は時間に対してだんだん大きくなるように設計されるため、アニーリングの後半ではレプリカ間の相互作用が強くなり、最終的には全てのスピンが同じ向きに揃うようになります。これは最初は量子の重ね合わせを独立なレプリカで表しておき、だんだん同じ向きに揃えられていくことで一つの状態に収束していくという描像をイメージしてもらうといいかもしれません。

最後に

Jijでは上記にあるような量子アニーリングに関する技術はもちろん数理最適化を使ったソリューション開発を行っています。
ぜひ以下のリンクからJijの求人情報をご覧ください！

＼Rustエンジニア・数理最適化エンジニア募集中！／
株式会社Jijでは、数学や物理学のバックグラウンドを活かし、量子計算と数理最適化のフロンティアで活躍するRustエンジニア、数理最適化エンジニアを募集しています！
詳細は下記のリンクからご覧ください。皆さんのご応募をお待ちしております！
Rustエンジニア: https://open.talentio.com/r/1/c/j-ij.com/pages/51062
数理最適化エンジニア: https://open.talentio.com/r/1/c/j-ij.com/pages/75132