$$
\def\bra#1{\mathinner{\left\langle{#1}\right|}}
\def\ket#1{\mathinner{\left|{#1}\right\rangle}}
\def\braket#1#2{\mathinner{\left\langle{#1}\middle|#2\right\rangle}}
$$
概要
アニーリング方式の量子コンピュータで用いている手法を古典計算機で疑似的に再現することができます。ここでは量子コンピュータのアルゴリズムをそのまま古典計算に近似した量子モンテカルロ法と呼ばれるものの導出を行います。
量子アニーリングでは横磁場をかけて局所解からの脱出を試みますが、熱運動を利用して局所解から脱出するシミュレーテッドアニーリングという手法もあります。こちらは別の記事で紹介します。
参考
- 量子アニーリングを用いたクラスタ分析 田中宗，栗原賢一，宮下精二(http://www.shutanaka.com/papers_files/ShuTanaka_DEXSMI_10.pdf)
- Can quantum Monte Carlo simulate quantum annealing? (2017) Evgeny Andriyash, Mohammad H. Amin (https://arxiv.org/abs/1703.09277)
- 量子アニーリングで組合せ最適化
手順
ハイゼンベルクモデルのハミルトニアン
アニーリング方式の量子コンピュータではハイゼンベルクモデルという、粒子に備わっているスピンの物理現象のモデルを使用します。結晶の表面の2次元平面にたくさん粒子が格子状に並んでいるとします（厳密には1次元の列を2次元平面に並べているだけであり、したがってこれは1次元ハイゼンベルクモデルとなります）。それぞれの粒子はスピン相互作用がはたらくので、（相互作用の係数が正のとき）スピンが同じ向きにそろう力が働きます。
用いるスピンは$\mathrm{SU}(2)$群の1/2表現になっています。$\mathrm{SU}(2)$群はパウリ行列$\hat \sigma^x,\hat \sigma^y,\hat \sigma^z$で生成されます。観測軸をz軸に取ることにすると、スピンの固有状態は$\hat \sigma^z$の固有状態（固有ベクトル）になっています。
今回はスピンのアップを$+1$、ダウンを$-1$と表記します。すると1粒子状態に作用する演算子$\hat \sigma^z$の作用は次のように書けます。
\hat \sigma^z\ket{1} = \ket{1}, \ \hat \sigma^z\ket{-1} = -\ket{-1}.
複数の粒子がある場合は状態は1粒子状態のテンソル積になります。この記事では格子に並んだ粒子$N$個に番号をつけます。
例えば上図の状態は$N=15$で状態は次のように書けます。
\ket{\boldsymbol \sigma} = \ket{1}\ket{1}\ket{-1}\ket{1}\ket{-1}\ket{1}\ket{1}\ket{-1}\ket{-1}\ket{1}\ket{-1}\ket{1}\ket{-1}\ket{-1}\ket{1}.
$i$番目のスピンに作用する演算子には添え字$i$をつけて$\hat \sigma^z_i$と表記します。スピン間の相互作用をモデル化したハミルトニアンは次のように書けます。
\begin{align}
\hat H = -\sum_{i,j=1}^N J_{ij}\hat\sigma_i^z \hat\sigma_j^z - \sum_{i=1}^N h_i \hat\sigma_i^z -\Gamma \sum_{i=1}^N \hat\sigma_i^x.
\end{align}
第一項は二つのスピンの向きによって決まる相互作用を表しています。第二項は格子全体に働く力でを表していて、係数が正ならば格子のスピン全体を上に向かせようとする力が働きます。第三項は横磁場の項でz軸方向のスピンの向きを強制的にひっくり返す力を表しています。量子アニーリングの肝は第三項です。前二項がスピンを揃えようと働く一方で、第三項がときどき向きをばらばらにします。$\Gamma$を徐々に小さくすることで一番エネルギーが小さい真のスピンの向き（ポテンシャルのglobal minimum）にたどり着くことができます。
量子モンテカルロ法の導出
方針ですが、上記のハミルトニアンの分配関数を求めます。計算していくとわかりますが、1次元量子ハイゼンベルクモデルの分配関数が2次元古典ハイゼンベルクモデルの分配関数になることがわかります。そこで、2次元古典ハイゼンベルクモデルの分配関数から古典ハイゼンベルクモデルのハミルトニアンを読み取ります。この古典ハイゼンベルクモデルを古典コンピュータ上に実装すれば量子モンテカルロ法となります。
まず、ハミルトニアンを対角化可能部分とそうでない部分に分けます。
\begin{align}
\hat H &= \hat H_0 + \hat H_1\\
\hat H_0 & = -\sum_{i,j=1}^N J_{ij}\hat\sigma_i^z \hat\sigma_j^z - \sum_{i=1}^N h_i \hat\sigma_i^z \\
\hat H_1 & =-\Gamma \sum_{i=1}^N \hat\sigma_i^x.
\end{align}
ここで、$\Gamma$は時間に依存しないものとします。分配関数は次のように書けます。
\begin{align}
Z &= \sum_{\boldsymbol{\sigma_1}}\bra{\boldsymbol{\sigma_1}}e^{-\beta \hat H}\ket{\boldsymbol{\sigma_1}}\\
&= \sum_{\boldsymbol{\sigma_1}}\bra{\boldsymbol{\sigma_1}}e^{-\beta \hat H_0-\beta \hat H_1}\ket{\boldsymbol{\sigma_1}}\\
&= \sum_{\boldsymbol{\sigma_1}}\bra{\boldsymbol{\sigma_1}}\left(e^{-\frac{\beta \hat H_0}{m}-\frac{\beta \beta \hat H_1}{m}}\right)^m\ket{\boldsymbol{\sigma_1}} + \mathcal{O}\left(\frac{\beta^2}{m}\right).
\end{align}
最後の等号は鈴木・トロッター分解を用いました。
\begin{align}
e^{x(\hat A+ \hat B)} = \left(e^{\frac{x \hat A}{L}}e^{\frac{x \hat B}{L}}\right)^L + \mathcal{O}\left(\frac{x^2}{L}\right).
\end{align}
$L\rightarrow \infty$で誤差なく両辺一致しますが、古典コンピュータ上に無限の変数は実装できないので、実装上は適当に大きな数を選びます。
それぞれの$e^{\cdots}$の間に完全系
\sum_{\boldsymbol{\sigma_2}}\ket{\boldsymbol{\sigma_2}}\bra{\boldsymbol{\sigma_2}}=1,\sum_{\boldsymbol{\sigma_3}}\ket{\boldsymbol{\sigma_3}}\bra{\boldsymbol{\sigma_3}}=1,\cdots ,\sum_{\boldsymbol{\sigma_{2m}}}\ket{\boldsymbol{\sigma_{2m}}}\bra{\boldsymbol{\sigma_{2m}}}=1,
を挿入します。
\begin{align}
Z = \sum_{\boldsymbol{\sigma_1},\boldsymbol{\sigma_2},\cdots,\boldsymbol{\sigma_{2m}}}
\bra{\boldsymbol{\sigma_1}}e^{-\frac{\beta \hat H_0}{m}}\ket{\boldsymbol{\sigma_2}}
\bra{\boldsymbol{\sigma_2}}e^{-\frac{\beta \hat H_1}{m}}\ket{\boldsymbol{\sigma_3}}
\bra{\boldsymbol{\sigma_3}}e^{-\frac{\beta \hat H_0}{m}}\ket{\boldsymbol{\sigma_4}}
\bra{\boldsymbol{\sigma_4}}e^{-\frac{\beta \hat H_1}{m}}\ket{\boldsymbol{\sigma_5}}
\cdots
\bra{\boldsymbol{\sigma_{2m-1}}}e^{-\frac{\beta \hat H_0}{m}}\ket{\boldsymbol{\sigma_{2m}}}
\bra{\boldsymbol{\sigma_{2m}}}e^{-\frac{\beta \hat H_1}{m}}\ket{\boldsymbol{\sigma_1}}.
\end{align}
$\hat H_0$の部分について見てみます。
\begin{align}
\bra{\boldsymbol{\sigma_{k}}} e^{-\frac{\beta \hat H_0}{m}} \ket{\boldsymbol{\sigma_{l}}}
=\bra{\boldsymbol{\sigma_{k}}} \exp\left[
\frac{\beta}{m}\left( \sum_{i,j=1}^N J_{ij}\hat\sigma_i^z \hat\sigma_j^z + \sum_{i=1}^N h_i \hat\sigma_i^z \right)
\right] \ket{\boldsymbol{\sigma_{l}}}.
\end{align}
$k$番目の完全系の状態には添え字$k$をつけます。$N$個の1粒子状態それぞれにも添え字$k$をつけます。
\begin{align}
\ket{\boldsymbol{\sigma_{k}}} = \ket{\sigma_{1,k}}\ket{\sigma_{2,k}}\cdots \ket{\sigma_{i,k}}\cdots \ket{\sigma_{N,k}}.
\end{align}
$k$番目の完全系の$i$個目のスピンに対するスピン演算子の固有値を$\hat \sigma_i^z\ket{\sigma_{i,k}}=\sigma_{i,k}\ket{\sigma_{i,k}}$と書くことにします。$\sigma_{i,k}$は$\sigma_i^z$の固有値である$1$または$-1$をとる変数です。
\begin{align}
\hat \sigma_i^z \ket{\boldsymbol{\sigma_{k}}}
&= \ket{\sigma_{1,k}}\ket{\sigma_{2,k}}\cdots \left( \hat \sigma_i^z\ket{\sigma_{i,k}}\right)\cdots \ket{\sigma_{N,k}}\\
&= \ket{\sigma_{1,k}}\ket{\sigma_{2,k}}\cdots \left(\sigma_{i,k}\ket{\sigma_{i,k}}\right)\cdots \ket{\sigma_{N,k}}\\
&= \sigma_{i,k}\ket{\boldsymbol{\sigma_{k}}}.
\end{align}
よって次のように書けます。
\begin{align}
\bra{\boldsymbol{\sigma_{k}}} e^{-\frac{\beta \hat H_0}{m}} \ket{\boldsymbol{\sigma_{l}}}
= \exp\left[
\frac{\beta}{m}\left( \sum_{i,j=1}^N J_{ij}\sigma_{i,k} \sigma_{j,k} + \sum_{i=1}^N h_i \sigma_{i,k} \right)
\right] \delta(\boldsymbol{\sigma_{k}},\boldsymbol{\sigma_{l}}).
\end{align}
ここで$\delta$はクロネッカーのデルタです。
次に、$\hat H_1$の部分について見てみます。
\begin{align}
\bra{\boldsymbol{\sigma_{k}}} e^{-\frac{\beta \hat H_1}{m}} \ket{\boldsymbol{\sigma_{l}}}
&= \bra{\boldsymbol{\sigma_{k}}} \exp\left[
\frac{\beta \Gamma}{m}\sum_{i=1}^N \hat\sigma_i^x
\right] \ket{\boldsymbol{\sigma_{l}}} \\
&= \bra{\boldsymbol{\sigma_{k}}} \left(\prod_{i=1}^N \exp\left[
\frac{\beta \Gamma}{m}\hat\sigma_i^x
\right] \right)\ket{\boldsymbol{\sigma_{l}}} \\
&= \bra{\boldsymbol{\sigma_{k}}} \left(\prod_{i=1}^N \left(
\cosh\frac{\beta \Gamma}{m} + \hat \sigma^x_i \sinh\frac{\beta \Gamma}{m}
\right) \right)\ket{\boldsymbol{\sigma_{l}}} \\
&=\prod_{i=1}^N \bra{\sigma_{i,k}} \left(
\cosh\frac{\beta \Gamma}{m} + \hat \sigma^x_i \sinh\frac{\beta \Gamma}{m}
\right)\ket{\sigma_{i,l}}.
\end{align}
$\hat \sigma^x_i$は$i$番目のスピンをひっくり返す演算子なので次のように場合分けできます。
\begin{align}
\bra{\sigma_{i,k}} \left(
\cosh\frac{\beta \Gamma}{m} + \hat \sigma^x_i \sinh\frac{\beta \Gamma}{m}
\right)\ket{\sigma_{i,l}}
=\left\{
\begin{array}{ll}
\cosh\frac{\beta \Gamma}{m} & (\sigma_{i,k}=\sigma_{i,l})\\
\sinh\frac{\beta \Gamma}{m} & (\sigma_{i,k}=-\sigma_{i,l})
\end{array}
\right..
\end{align}
双曲線関数は次のように書き直すことができます。
\begin{align}
\cosh x &= \left( \frac{1}{2}\sinh 2x \right)^{\frac{1}{2}}e^{-\frac{1}{2}\ln(\tanh x)} \\
\sinh x &= \left( \frac{1}{2}\sinh 2x \right)^{\frac{1}{2}}e^{\frac{1}{2}\ln(\tanh x)} .
\end{align}
$\sigma_{i,k}$は$1$または$-1$をとることと場合分けを鑑みると次のように式変形できます。
\begin{align}
\bra{\boldsymbol{\sigma_{k}}} e^{-\frac{\beta \hat H_1}{m}} \ket{\boldsymbol{\sigma_{l}}}
&=\prod_{i=1}^N \left( \frac{1}{2}\sinh 2\frac{\beta \Gamma}{m} \right)^{\frac{1}{2}}e^{-\frac{1}{2}\ln(\tanh \frac{\beta \Gamma}{m})\sigma_{i,k}\sigma_{i,l}}\\
&=\left( \frac{1}{2}\sinh 2\frac{\beta \Gamma}{m} \right)^{\frac{N}{2}}e^{-\frac{1}{2}\ln(\tanh \frac{\beta \Gamma}{m})\sum_{i=1}^N \sigma_{i,k}\sigma_{i,l}}.
\end{align}
$\hat H_0$と$\hat H_1$の計算を合わせると分配関数は次のようになります。
\begin{align}
Z &=\left( \frac{1}{2}\sinh 2\frac{\beta \Gamma}{m} \right)^{\frac{mN}{2}} \sum_{\boldsymbol{\sigma_1},\boldsymbol{\sigma_2},\cdots,\boldsymbol{\sigma_{2m}}}
\prod_{k=1}^m\exp\left[
\frac{\beta}{m}\left( \sum_{i,j=1}^N J_{ij}\sigma_{i,2k-1} \sigma_{j,2k-1} + \sum_{i=1}^N h_i \sigma_{i,2k-1} \right)
-\frac{1}{2}\ln\left(\tanh \frac{\beta \Gamma}{m}\right)\sum_{i=1}^N \sigma_{i,2k}\sigma_{i,2k+1}\right] \delta(\boldsymbol{\sigma_{2k-1}},\boldsymbol{\sigma_{2k}})\\
&=\left( \frac{1}{2}\sinh 2\frac{\beta \Gamma}{m} \right)^{\frac{mN}{2}} \sum_{\boldsymbol{\sigma_1},\boldsymbol{\sigma_3},\cdots,\boldsymbol{\sigma_{2m-1}}}
\prod_{k=1}^m\exp\left[
\frac{\beta}{m}\left( \sum_{i,j=1}^N J_{ij}\sigma_{i,2k-1} \sigma_{j,2k-1} + \sum_{i=1}^N h_i \sigma_{i,2k-1} \right)
-\frac{1}{2}\ln\left(\tanh \frac{\beta \Gamma}{m}\right)\sum_{i=1}^N \sigma_{i,2k-1}\sigma_{i,2k+1}\right].
\end{align}
ここで$\boldsymbol{\sigma_{2m+1}}:=\boldsymbol{\sigma_{1}}$としました。総和はクロネッカーのデルタによって半分になりました。
分かりやすさのために添え字の付け替えを行います。完全系の添え字に1を足して2で割ったものを新たな添え字とします。
\begin{align}
Z &=\left( \frac{1}{2}\sinh 2\frac{\beta \Gamma}{m} \right)^{\frac{mN}{2}} \sum_{\boldsymbol{\sigma_1},\boldsymbol{\sigma_2},\cdots,\boldsymbol{\sigma_{m}}}
\prod_{k=1}^m\exp\left[
\frac{\beta}{m}\left( \sum_{i,j=1}^N J_{ij}\sigma_{i,k} \sigma_{j,k} + \sum_{i=1}^N h_i \sigma_{i,k} \right)
-\frac{1}{2}\ln\left(\tanh \frac{\beta \Gamma}{m}\right)\sum_{i=1}^N \sigma_{i,k}\sigma_{i,k+1}\right]\\
&=\left( \frac{1}{2}\sinh 2\frac{\beta \Gamma}{m} \right)^{\frac{mN}{2}} \sum_{\boldsymbol{\sigma_1},\boldsymbol{\sigma_2},\cdots,\boldsymbol{\sigma_{m}}}
\exp\left[
\frac{\beta}{m}\left( \sum_{i,j=1}^N\sum_{k=1}^m J_{ij}\sigma_{i,k} \sigma_{j,k} + \sum_{i=1}^N \sum_{k=1}^m h_i \sigma_{i,k} \right)
-\frac{1}{2}\ln\left(\tanh \frac{\beta \Gamma}{m}\right)\sum_{i=1}^N \sum_{k=1}^m \sigma_{i,k}\sigma_{i,k+1}\right].
\end{align}
$\boldsymbol{\sigma_{2m+1}}=\boldsymbol{\sigma_{1}}$の関係は$\boldsymbol{\sigma_{m+1}}=\boldsymbol{\sigma_{1}}$と改められます。
分配関数の大きさそのものには意味がないので、結局分配関数は次のように書けます。
\begin{align}
Z
=\sum_{\boldsymbol{\sigma_1},\boldsymbol{\sigma_2},\cdots,\boldsymbol{\sigma_{m}}}
\exp\left[ -
\frac{\beta}{m}\left( -\sum_{i,j=1}^N\sum_{k=1}^m J_{ij}\sigma_{i,k} \sigma_{j,k} -\sum_{i=1}^N \sum_{k=1}^m h_i \sigma_{i,k}
+\frac{m}{2\beta}\ln\left(\tanh \frac{\beta \Gamma}{m}\right)\sum_{i=1}^N \sum_{k=1}^m \sigma_{i,k}\sigma_{i,k+1}\right)\right].
\end{align}
これは2次元ハイゼンベルクモデルの特殊な相互作用の場合の分配関数になっています（$\sigma$に二つの添え字がついているので、これは本当に2次元のモデルです）。逆温度を$\beta'=\beta/m$で再定義し、分かりやすさのために$\tanh$を$\coth$で書き直します。ハミルトニアンは次のように書くことができます。
\begin{align}
H
&=
-\sum_{i,j=1}^N\sum_{k=1}^m J_{ij}\sigma_{i,k} \sigma_{j,k} -\sum_{i=1}^N \sum_{k=1}^m h_i \sigma_{i,k}
-\frac{1}{2\beta'}\ln\left(\coth (\beta'\Gamma)\right)\sum_{i=1}^N \sum_{k=1}^m \sigma_{i,k}\sigma_{i,k+1}\\
&=-\sum_{i,j=1}^N\sum_{k=1}^m J_{ij}\sigma_{i,k} \sigma_{j,k} -\sum_{i=1}^N \sum_{k=1}^m h_i \sigma_{i,k}
-\frac{T'}{2}\ln\left(\coth \left(\frac{\Gamma}{T'}\right)\right)\sum_{i=1}^N \sum_{k=1}^m \sigma_{i,k}\sigma_{i,k+1}
\end{align}
$T'$は有効温度です。第三項は量子ハイゼンベルクモデルの横磁場の効果を表しています。
量子モンテカルロ法ではこのハミルトニアン最小の値をとるようなスピンの配位を求めることになります。
量子モンテカルロ法（アニーリング）の手順
- 温度$T'$と横磁場の強さ$\Gamma$を決めます。
- $m$個の層のそれぞれの$N$個のスピンの配位をランダムに決めます（初期化）。
- ハミルトニアンの値を計算し、エネルギーを求め、$E$とします。ランダムに層とその層の中のスピンを選びスピンの向きを変えます。ハミルトニアンを計算し系のエネルギーを求め、$E'$とします。確率$\mathrm{min}(1,e^{\frac{E-E'}{T'}})$でスピンの向きの変更を受容し、状態変更をコミットします。棄却した場合はスピンの向きの変更をなかったことにします。
- ある決められた回数（モンテカルロステップ数）だけ3.を繰り返しエネルギーを下げていきます。
- 横磁場の強さ$\Gamma$を小さくします。
- 4.と5.をある決められた回数（アニーリングステップ数）行います。
- 最後に第三項を除いた各層のエネルギーを求め、一番エネルギーの低い層の配位を解として採用します。
横磁場の寄与
ハミルトニアンの第三項の係数の振る舞いについて見てみましょう。
温度固定で$\Gamma$依存性を見てみます。
g(x) =\ln(\coth(x)),
のグラフは下記画像になります。$\Gamma$が小さくなるほど無限大に近づく関数になっています。
微分係数の絶対値が大きいので、$\Gamma = 1.0\times (0.9)^x$とすれば、なだらかになります。
微分係数はあるステップ以上でほぼ一定となる。
以上の係数の振る舞いは実装のときに重要になります。あまりに急激に変化すると局所解につかまってしまいます。一方でゆっくり変化するとなかなか解に収束せず、計算に時間がかかってしまいます。
おまけに$\Gamma$固定での温度依存性も記載しておきます。
f(x) = \frac{x}{2}\ln(\coth(\frac{1}{x})),
のグラフは下記画像になります。温度が下がるほどゼロに近づく関数になっています。
横磁場の効果の解釈
\begin{align}
H
&=-\sum_{i,j=1}^N\sum_{k=1}^m J_{ij}\sigma_{i,k} \sigma_{j,k} -\sum_{i=1}^N \sum_{k=1}^m h_i \sigma_{i,k}
-\frac{T'}{2}\ln\left(\coth \left(\frac{\Gamma}{T'}\right)\right)\sum_{i=1}^N \sum_{k=1}^m \sigma_{i,k}\sigma_{i,k+1}
\end{align}
ハミルトニアンの第三項は横磁場の効果と解釈できますが、もう少し詳しく見てみます。元々のハミルトニアンは下記のように書けました。
\begin{align}
\hat H = -\sum_{i,j=1}^N J_{ij}\hat\sigma_i^z \hat\sigma_j^z - \sum_{i=1}^N h_i \hat\sigma_i^z -\Gamma \sum_{i=1}^N \hat\sigma_i^x.
\end{align}
これは2次元平面に格子状にスピンが並んでいる様子を表していました。添え字$i,j$がスピンの位置を表しています。一方で、求めた古典ハイゼンベルクモデルは添え字$i,j$に加え、添え字$k$を持っています。これが時間方向の次元を表しています。
ハミルトニアンの第一項と第二項にも添え字$k$がありますが、同じ層の相互作用になっていますので、ただのハイゼンベルクモデルを複数用意しただけになっています。各層間の相互作用はハミルトニアンの第三項、すなわち横磁場の効果を表している項によって行われます。これは横磁場によってスピンの配位が時間発展していくことを表しています。
分配関数を求めるときの鈴木・トロッター分解、
\begin{align}
Z = \sum_{\boldsymbol{\sigma_1},\boldsymbol{\sigma_2},\cdots,\boldsymbol{\sigma_{m}}}
\bra{\boldsymbol{\sigma_1}}e^{-\frac{\beta \hat H_0}{m}}e^{-\frac{\beta \hat H_1}{m}}\ket{\boldsymbol{\sigma_2}}
\bra{\boldsymbol{\sigma_2}}e^{-\frac{\beta \hat H_0}{m}}e^{-\frac{\beta \hat H_1}{m}}\ket{\boldsymbol{\sigma_3}}
\cdots
\bra{\boldsymbol{\sigma_{m}}}e^{-\frac{\beta \hat H_0}{m}}e^{-\frac{\beta \hat H_1}{m}}\ket{\boldsymbol{\sigma_1}},
\end{align}
は量子状態を発展させる時間を$m$個に分割したことに相当します。また、平衡状態を表現するために時間方向は周期的になります。周期は逆温度に一致します。$\sum_{\boldsymbol{\sigma_1},\boldsymbol{\sigma_2},\cdots,\boldsymbol{\sigma_{m}}}$が各層の重ね合わせを表しています。つまり分配関数は次のように図示できます。
重ね合わせのうち、最もエネルギーが低くなる状態の実現確率が高くなります。また、ハミルトニアンの第三項（横磁場の効果）の係数は$\Gamma$を小さくするほど大きくなるので、各層の状態は同じ状態に近づいていきます。つまり各層でエネルギーの低い状態を探索して、各層の状態を揃えることでエネルギーを低くしていくのです。$m\rightarrow \infty$とすることで探索パターンが増えるので、最低エネルギー状態にたどり着きやすくなります。また、このとき量子モデルの分配関数に一致します。
古典コンピューターでは重ね合わせを作り出すことができないので、モンテカルロで状態探索を行います。
実装
ではこの量子モンテカルロ法がどのように役立つか、実装はどうすればいいのかというと、、、これについては他の方々が詳しく解説していますので、ご紹介にとどめます。
- 巡回セールスマン問題への応用：量子アニーリングで組合せ最適化
- クラスタ分析：「量子アニーリングを用いた クラスタ分析」(2015) 田中宗 (http://www.butsuri.it-chiba.ac.jp/~yasutake/matter/tanaka.pdf)
- 因数分解：D-Waveで素因数分解をした
他には包括的な紹介記事もあります。
おわりに
今回は量子モンテカルロ法で用いるハミルトニアンの導出、モンテカルロ法の手順を紹介しました。量子モンテカルロ法はハイゼンベルクモデルを解くのに利用されます。そして最適化問題のいくつかはハイゼンベルクモデルと等価になるので、最適化問題を解くのに量子モンテカルロ法を活用することができます。
量子モンテカルロ法は問題のサイズ$N$が大きくなったとき力を発揮します。(NASA&Googleの量子コンピュータは「一億倍速い」の論文（量子アニーリング）)