補遺: 量子アニーリングのシミュレーション手法

量子アニーリングを調べるとき、実機を動かすだけでは足りません。何が量子効果なのか、何が単なる熱緩和なのか、どの難しさが問題インスタンスに由来し、どの難しさが装置や埋め込みに由来するのかを切り分ける必要があります。そのために、古典計算機上でいくつかの「似せた力学」を走らせます。

ここで大切なのは、すべてのシミュレーションが同じものを再現しているわけではない、という点です。シミュレーテッドアニーリングは、最終的な古典エネルギー地形を熱で登り降りする方法です。量子モンテカルロによる量子アニーリングは、横磁場を含む量子系の平衡分布を、虚時間方向を持つ古典模型としてサンプルします。連続時間量子モンテカルロは、そのTrotter分割を極限まで細かくした表現です。スピンベクトルモンテカルロは、各量子ビットを連続的な古典ベクトルとして置き換え、半古典的なエネルギー地形の上をMetropolis法で動かします。

同じ「量子アニーリングを模擬する」と言っても、答えている問いが違います。実時間Schrodinger方程式を解いているのか、瞬間平衡分布をサンプルしているのか、熱活性化だけで説明できるかを試しているのでしょうか。この違いを曖昧にしたままベンチマークを見ると、量子スピードアップの議論はすぐに混乱します。

1. 基準となる量子アニーリングハミルトニアン

まず、横磁場Ising型の量子アニーリングを基準にします。目的ハミルトニアンを

$${\hat{H}}_{\mathrm{P}}=\sum _i{h}_i{\hat{\sigma }}_i^z+\sum _{i<j}{J}_{ij}{\hat{\sigma }}_i^z{\hat{\sigma }}_j^z$$

とします。符号の取り方は文献によって少し違いますが、ここでは $H_{\mathrm{P}}$ の小さい状態を探すという約束にします。横磁場ドライバーは

$${\hat{H}}_{\mathrm{D}}=-\sum _i{\hat{\sigma }}_i^x$$

です。典型的な時間依存ハミルトニアンは

$$\hat{H}(s)=A(s){\hat{H}}_{\mathrm{D}}+B(s){\hat{H}}_{\mathrm{P}},0\le s\le 1$$

と書けます。$A(s)$ は最初大きく最後に小さくなり、$B(s)$ は最初小さく最後に大きくなります。D-Wave型の実機で使われる書き方も基本的にはこの形です。

古典配置 $\mathit{\sigma }$ に対する目的エネルギーを

$$E_{\mathrm{P}}(\mathit{\sigma })=\sum _i{h}_i{\sigma }_i+\sum _{i<j}{J}_{ij}{\sigma }_i{\sigma }_j$$

と書くと、アニーリングの最後に欲しいのは $E_{\mathrm{P}}$ の低い配置です。評価では、最低エネルギー状態に到達する確率

$$p_{\mathrm{s}}=Pr[{\mathit{\sigma }}_{\mathrm{o}\mathrm{u}\mathrm{t}}\in \mathcal{G}]$$

を測ることが多いです。$\mathcal{G}$ は基底状態の集合です。何度も独立に走らせて少なくとも1回成功する確率を $p_{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}}$ にしたいなら、必要な反復回数は

$$R=\frac{\log (1-p_{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}})}{\log (1-p_{\mathrm{s}})}$$

です。1回の実行時間を $t_{\mathrm{r}\mathrm{u}\mathrm{n}}$ とすれば、time-to-solutionは

$$\mathrm{T}\mathrm{T}\mathrm{S}=t_{\mathrm{r}\mathrm{u}\mathrm{n}}\frac{\log (1-p_{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}})}{\log (1-p_{\mathrm{s}})}$$

となります。よく使われるのは $p_{\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{g}}=0.99$ です。ただし、この式を使うときには、何を $t_{\mathrm{r}\mathrm{u}\mathrm{n}}$ に含めるのかを明確にしなければなりません。実機ならアニーリング時間だけなのか、プログラム時間、読み出し時間、埋め込み計算まで含めるのでしょうか。古典シミュレーションならCPU時間なのか、Monte Carlo sweep数なのでしょうか。比較の結論は、この選び方にかなり敏感です。

2. シミュレーテッドアニーリング

シミュレーテッドアニーリングは、目的エネルギー $E_{\mathrm{P}}$ に温度を入れて動かします。時刻、あるいはMonte Carlo sweepの番号を $r$ とし、そのときの逆温度を ${\beta }_r$ とします。固定した ${\beta }_r$ での目標分布は

$${\pi }_{{\beta }_r}(\mathit{\sigma })=\frac{1}{Z_{{\beta }_r}}\exp [-{\beta }_r{E}_{\mathrm{P}}(\mathit{\sigma })]$$

です。${\beta }_r$ を小さい値から大きい値へ上げていけば、分布は高温の広い探索から低温の低エネルギー探索へ移ります。

1スピン反転を提案するMetropolis法なら、配置 $\mathit{\sigma }$ から ${\mathit{\sigma }}^{\prime }$ へ移る受理確率は

$$A(\mathit{\sigma }\to {\mathit{\sigma }}^{\prime })=min\{1,\exp [-{\beta }_r\mathrm{\Delta }E]\},\mathrm{\Delta }E=E_{\mathrm{P}}({\mathit{\sigma }}^{\prime })-E_{\mathrm{P}}(\mathit{\sigma })$$

です。エネルギーが下がるなら受け入れ、上がるならBoltzmann因子で受け入れます。これは「山を熱で越える」力学です。

SAは、量子アニーリングの比較対象として必ず出てきます。理由は単純です。もし量子アニーリングが本当に役に立つなら、少なくとも素朴な熱的探索とは違う振る舞いを示すはずだからです。ただし、SAには多くの改良版があります。温度スケジュール、提案分布、クラスタ更新、交換Monte Carlo、問題構造を使う局所探索を変えるだけで性能は大きく変わります。したがって「SAに勝った」という主張は、何のSAに勝ったのかまで言わないと弱いです。

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

量子モンテカルロによる量子アニーリングは、しばしばSQA、simulated quantum annealingと呼ばれます。ここでの発想は、横磁場を含む量子Ising模型の分配関数を、虚時間方向を持つ古典Ising模型へ写すことです。

一つの $s$ でハミルトニアンを固定し、簡単のため

$$\hat{H}=E_{\mathrm{P}}({\hat{\mathit{\sigma }}}^z)-\mathrm{\Gamma }\sum _i{\hat{\sigma }}_i^x$$

と書きます。$\mathrm{\Gamma }$ が横磁場の強さです。分配関数は

$$Z=\mathrm{T}\mathrm{r}\exp (-\beta \hat{H})$$

です。$E_{\mathrm{P}}({\hat{\mathit{\sigma }}}^z)$ と横磁場項は交換しないので、そのままでは古典的な確率分布になりません。そこで逆温度方向を $M$ 個のTrotterスライスに分けます。

$$Z_M=\sum _{\{{\sigma }_{i,k}\}}\exp [-\frac{\beta }{M}\sum _{k=1}^M{E}_{\mathrm{P}}({\mathit{\sigma }}_k)+K\sum _{i=1}^N\sum _{k=1}^M{\sigma }_{i,k}{\sigma }_{i,k+1}]$$

ここで ${\sigma }_{i,k}$ は、サイト $i$、Trotterスライス $k$ の古典スピンです。周期境界条件 ${\sigma }_{i,M+1}={\sigma }_{i,1}$ を入れます。横磁場は、虚時間方向の隣接スライスを結ぶ結合

$$K=\frac{1}{2}\log \coth \left(\frac{\beta \mathrm{\Gamma }}{M}\right)$$

へ変わります。

この式は、2日目のSQAで見た式そのものです。量子揺らぎは、古典スピンが虚時間方向にどれだけ切り替わるかとして表現されます。$\mathrm{\Gamma }$ が小さいと $K$ は大きくなり、同じサイトのTrotterコピーはそろいやすいです。$\mathrm{\Gamma }$ が大きいと虚時間方向に切れ目が入りやすくなり、量子揺らぎが強くなります。

SQAの実装では、$s$ を少しずつ変えながら、各 $s$ でこの古典模型に対するMonte Carlo更新を行います。1 sweepは、全サイト・全Trotterスライスに対して更新を1回ずつ試す単位として定義されることが多いです。

$$s_0\to s_1\to \cdots \to s_R,\mathrm{\Gamma }(s_0)\gg 1,\mathrm{\Gamma }(s_R)\simeq 0$$

この手順は、量子系の平衡分布を追うMonte Carloです。実時間Schrodinger方程式を解いているわけではありません。したがって、SQAが速いから実機の量子アニーリングも必ず速い、とは言えません。一方で、stoquasticな横磁場Ising模型におけるトンネル的な経路を、虚時間インスタントンとしてよく捉えることがあります。そのため、SQAは「量子アニーリング的な障壁通過が、古典計算機でどこまで模擬できるか」を調べる強い基準になります。

4. 連続時間量子モンテカルロ

Trotter数 $M$ を有限にしたSQAでは、虚時間刻み

$$\mathrm{\Delta }\tau =\frac{\beta }{M}$$

が残ります。$M$ が十分大きくなければ、離散化誤差が物理的でない利点や欠点を作ることがあります。そこで $\mathrm{\Delta }\tau \to 0$ の極限を直接扱うのが、連続時間量子モンテカルロです。

連続時間では、各スピンは虚時間 $\tau \in [0,\beta ]$ の上で区分的に一定のworldlineを持ちます。スピンが反転する点をkinkと呼びます。全worldlineの集合を $\mathcal{W}$、その中のkink数を $K(\mathcal{W})$ と書くと、分配関数は概念的に

$$Z=\sum _{\mathcal{W}}{\mathrm{\Gamma }}^{K(\mathcal{W})}\exp [-{\int }_0^{\beta }{E}_{\mathrm{P}}(\mathit{\sigma }(\tau ))d\tau ]$$

と書けます。周期境界条件のため、各スピンのworldlineは $\tau =0$ と $\tau =\beta$ でつながります。横磁場 $\mathrm{\Gamma }$ はkinkを作る重みとして現れ、目的エネルギーは虚時間に沿った作用として現れます。

連続時間QMCでは、格子点上の ${\sigma }_{i,k}$ を更新する代わりに、各スピンのworldlineを「どの虚時間区間で上向きか、どこで反転するか」という線分のリストとして持ちます。たとえば、スピン $i$ 以外のworldlineを固定すると、スピン $i$ が見る局所場は

$$b_i(\tau )=h_i+\sum _{j(\ne i)}{J}_{ij}{\sigma }_j(\tau )$$

です。このとき、スピン $i$ の一本のworldline ${\sigma }_i(\tau )$ の条件付き重みは

$$w_i[{\sigma }_i(\tau )\mid \{{\sigma }_j(\tau ){\}}_{j\ne i}]\propto {\mathrm{\Gamma }}^{K_i}\exp [-{\int }_0^{\beta }{b}_i(\tau ){\sigma }_i(\tau )d\tau ]$$

と見ればよいです。$K_i$ はそのworldline上のkink数です。横磁場 $\mathrm{\Gamma }$ は「反転点を作る重み」、縦方向のIsingエネルギーは「虚時間に沿って積分されるポテンシャル」として働きます。

もっとも素朴なMetropolis更新は、次のように書けます。まずスピン $i$ と二つの虚時刻 ${\tau }_a,{\tau }_b$ を選びます。区間 $({\tau }_a,{\tau }_b)$ の ${\sigma }_i(\tau )$ を反転させる、あるいは既存のkink位置を動かす候補worldline ${\mathcal{W}}^{\prime }$ を作ります。元のworldlineを $\mathcal{W}$ とすれば、受理確率は

$$\begin{array}{rl}{P}_{\mathrm{a}\mathrm{c}\mathrm{c}}(\mathcal{W}\to {\mathcal{W}}^{\prime })& =min\{1,R(\mathcal{W}\to {\mathcal{W}}^{\prime })\},\\ R(\mathcal{W}\to {\mathcal{W}}^{\prime })& ={\mathrm{\Gamma }}^{\mathrm{\Delta }K}\exp [-{\int }_0^{\beta }\{E_{\mathrm{P}}({\mathit{\sigma }}^{\prime }(\tau ))-E_{\mathrm{P}}(\mathit{\sigma }(\tau ))\}d\tau ]\frac{q({\mathcal{W}}^{\prime }\to \mathcal{W})}{q(\mathcal{W}\to {\mathcal{W}}^{\prime })}.\end{array}$$

です。$\mathrm{\Delta }K$ は候補で増減したkink数、$q$ は提案確率です。対称な提案なら最後の比は1になります。区間反転では端点にkinkが二つ挿入されることもあれば、既存のkinkが消えることもあります。実装では、端点が周期境界 $\tau =0,\beta$ をまたぐ場合も同じ重みで扱い、更新後に隣り合う同符号の線分をまとめておくと管理しやすいです。

1 sweepは、全スピンについてこのようなworldline更新を一通り試す単位として定義できます。局所更新だけではkinkが固まりやすいので、実際のベンチマークでは虚時間方向に長いクラスタを作る更新、worm型の更新、あるいは単一スピンの条件付きworldlineを熱浴的に引き直す更新が使われます。共通しているのは、有限個のTrotterスライスを明示的に持たず、kink位置を連続変数として動かす点です。

測定量もworldlineから直接計算します。縦磁化と横磁化は、たとえば

$$⟨{\sigma }_i^z⟩=⟨\frac{1}{\beta }{\int }_0^{\beta }{\sigma }_i(\tau )d\tau ⟩,⟨\sum _i{\sigma }_i^x⟩=⟨\frac{K(\mathcal{W})}{\beta \mathrm{\Gamma }}⟩$$

と評価できます。二つ目の式は、重みの ${\mathrm{\Gamma }}^{K(\mathcal{W})}$ を $\mathrm{\Gamma }$ で微分するとkink数が出ることから分かります。目的エネルギーは ${\beta }^{-1}{\int }_0^{\beta }{E}_{\mathrm{P}}(\mathit{\sigma }(\tau ))d\tau$ の平均です。

量子アニーリングを模擬するときは、$s=s_0,s_1,\dots ,s_R$ を用意し、各 $s_r$ で $\mathrm{\Gamma }=\mathrm{\Gamma }(s_r)$ と $E_{\mathrm{P}}$ の係数を更新します。前の $s_r$ で得たworldlineを次の $s_{r+1}$ の初期状態として使い、一定数のsweepを行います。これが連続時間版のSQAです。十分に平衡化してから測るなら平衡量の計算になり、平衡化を短くしてスケジュールに沿って進めるなら、Monte Carlo時間で定義された擬似的なアニーリングになります。

有限Trotter数に由来する誤差がないので、量子Ising模型の熱平衡性質を調べるには連続時間QMCがきれいです。

ただし、連続時間にしたから万能になるわけではありません。第一に、やはり虚時間の平衡サンプリングであって、実時間のユニタリ動力学そのものではありません。第二に、non-stoquastic ハミルトニアンでは符号問題が残ります。第三に、Monte Carlo sweep数と実機の物理時間は同じ単位ではありません。この区別は、2次元スピングラスのベンチマークで重要になりました。有限Trotter刻みのSQAが良く見せていた優位性が、連続時間極限では弱まることがあるからです。

5. スピンベクトルモンテカルロ

量子アニーリングの古典的な対照モデルとして、spin-vector Monte Carlo、略してSVMCがよく使われます。Shin--Smith--Smolin--Vaziraniの頭文字からSSSV modelとも呼ばれます。

SVMCでは、各量子ビットを量子状態として持つ代わりに、$xz$ 平面内の古典的な単位ベクトルで表します。角度を ${\theta }_i$ とし、

$$⟨{\hat{\sigma }}_i^z⟩=\cos {\theta }_i,⟨{\hat{\sigma }}_i^x⟩=\sin {\theta }_i$$

と置きます。すると、上の時間依存ハミルトニアンの期待値から半古典的なエネルギー関数

$$\begin{array}{rl}{E}_{\mathrm{S}\mathrm{V}\mathrm{M}\mathrm{C}}(s,\mathit{\theta })& =-A(s)\sum _i\sin {\theta }_i\\ & +B(s)(\sum _i{h}_i\cos {\theta }_i+\sum _{i<j}{J}_{ij}\cos {\theta }_i\cos {\theta }_j)\end{array}$$

が得られます。これは量子多体状態ではありません。各スピンを独立な古典ベクトルとして扱う平均場的な近似です。もとの量子ハミルトニアンにある重ね合わせ、干渉、エンタングルメント、多体トンネルは失われます。

実装では、まずアニーリング点 $s_0,s_1,\dots ,s_R$、各点でのsweep数、有効温度 $T_{\mathrm{e}\mathrm{f}\mathrm{f}}$、角度提案の幅を決めます。実機の古典対照として使うなら、$A(s),B(s)$ は実機のスケジュールに合わせ、温度は $k_{\mathrm{B}}{T}_{\mathrm{e}\mathrm{f}\mathrm{f}}$ としてエネルギー単位にそろえます。純粋なアルゴリズム比較なら、$\beta =1/T_{\mathrm{e}\mathrm{f}\mathrm{f}}$ を無次元パラメータとして固定してもよいです。

1回の局所更新では、サイト $i$ を選び、${\theta }_i$ を ${\theta }_i^{\prime }$ に変える候補を作ります。候補は $[0,\pi ]$ から一様に引いてもよいし、現在値の近くから ${\theta }_i^{\prime }={\theta }_i+\delta$ として引いてもよいです。後者では端点 $0,\pi$ で反射させると扱いやすいです。${\theta }_i$ だけを変えるとき、エネルギー差は全エネルギーを計算し直さなくても

$$\mathrm{\Delta }{E}_i=-A(s)\{\sin {\theta }_i^{\prime }-\sin {\theta }_i\}+B(s)(h_i+\sum _{j(\ne i)}{J}_{ij}\cos {\theta }_j)\{\cos {\theta }_i^{\prime }-\cos {\theta }_i\}$$

で求まります。括弧内は、古典的な $z$ 成分 $\cos {\theta }_j$ が作る局所場です。この差に対してMetropolis受理確率

$$P_{\mathrm{a}\mathrm{c}\mathrm{c}}=min\{1,\exp [-\beta \mathrm{\Delta }{E}_i]\}$$

を使います。これを全サイトに対して一通り試すと1 sweepです。サイトの順序を毎回ランダムにすれば、順序に由来する癖を弱められます。各 $s_r$ で決めたsweep数だけ更新し、最後の角度列を次の $s_{r+1}$ の初期状態にします。

初期状態は横磁場の基底状態に対応するので、全スピンを $x$ 方向、すなわち ${\theta }_i=\pi /2$ に置きます。完全に同じ角度から始めると対称性が残りすぎる場合があるので、数値実験では小さい乱数を加えることもあります。$s$ が進むと $B(s)$ のIsingエネルギーが強くなり、ベクトルは $z$ 方向へ倒れていきます。

最後の読み出しには二つの流儀があります。最適化のヒューリスティックとして使うなら

$${\sigma }_i^{\mathrm{o}\mathrm{u}\mathrm{t}}=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}(\cos {\theta }_i)$$

と決めればよいです。量子測定に近い読み出しを模擬したいなら、各スピンを独立に

$$Pr({\sigma }_i^{\mathrm{o}\mathrm{u}\mathrm{t}}=+1\mid {\theta }_i)=\frac{1+\cos {\theta }_i}{2},Pr({\sigma }_i^{\mathrm{o}\mathrm{u}\mathrm{t}}=-1\mid {\theta }_i)=\frac{1-\cos {\theta }_i}{2}$$

でサンプルします。アニーリング終端で $\cos {\theta }_i$ がほぼ $\pm 1$ まで倒れていれば、二つの読み出しはほとんど同じになります。途中で止める、または有限温度で角度が大きく揺れる場合には差が出ます。

SVMCが重要なのは、量子ではないのに、D-Wave実機の一部の統計をよく再現したことがあるからです。これは「実機は古典だ」と即断するためではなく、「成功確率分布が似ている」「特定の埋め込みで相関が高い」といった観測だけでは、量子性の証拠として十分でないことを示します。

SVMCは、半古典的なエネルギー地形の局所最小に沿って動きます。障壁を越えるには熱活性化が必要です。したがって、弱強クラスタ問題のように、古典ベクトルの道では高い障壁を越えなければならないが、量子多体状態では狭い障壁をトンネルできる、という状況を見分けるための対照群として有効です。

6. 何を再現し、何を捨てているか

4つの手法を並べると、次のように役割が分かれます。

手法	主に再現するもの	捨てているもの
SA	最終古典ハミルトニアン上の熱的探索	横磁場、虚時間経路、トンネル的障壁通過
SQA	stoquastic横磁場Ising模型の虚時間平衡分布	実時間位相、コヒーレント振動、物理時間との直接対応
連続時間QMC	Trotter誤差を除いた虚時間worldline分布	実時間動力学、non-stoquastic系の符号問題
SVMC	古典スピンベクトルの半古典的な熱活性化	エンタングルメント、多体トンネル、量子干渉

この表から分かるように、SQAとSVMCはどちらも「古典計算機上のモデル」ですが、かなり違うものです。SQAは虚時間方向を持ち、トンネル経路をworldlineの変形として表します。SVMCは連続角度の古典ベクトルを使い、障壁越えは熱活性化になります。量子アニーリングらしさを見たいなら、この二つを両方置くとよいです。SQAで説明できてSVMCで説明しにくいなら、少なくとも半古典的な熱活性化だけでは足りません。逆にSVMCでも説明できるなら、その観測だけを量子性の証拠にするのは危ういです。

もう一つの注意は、最適化アルゴリズムとしてのSQAと、物理実機のシミュレーションとしてのSQAを分けることです。有限Trotter数、有限sweep数、クラスタ更新、温度設定を工夫したSQAは、優秀な量子インスパイアード古典アルゴリズムになります。しかし、それが実機の開放系時間発展を忠実に再現しているとは限りません。このずれを利用して速く解くのは古典アルゴリズムとしては正当ですが、量子ハードウェアの物理を説明するモデルとしては慎重に扱う必要があります。

7. ベンチマークで何を比べるか

量子アニーリングのベンチマークでは、少なくとも次の量を分けて見る必要があります。

$$p_{\mathrm{s}}(N,t)=\text{サイズ }N\text{、実行時間 }t\text{ での成功確率}$$

$${ϵ}_{\mathrm{r}\mathrm{e}\mathrm{s}}=\frac{⟨E_{\mathrm{P}}⟩-E_{min}}{N}$$

$$\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{l}\mathrm{i}\mathrm{n}\mathrm{g}\sim \text{TTS}(N)\text{ の }N\text{ 依存性}$$

$p_{\mathrm{s}}$ は基底状態を当てる確率です。${ϵ}_{\mathrm{r}\mathrm{e}\mathrm{s}}$ は残留エネルギーで、最低エネルギーにどれだけ近いかを見ます。TTSのスケーリングは、サイズを大きくしたときに本当に差が広がるかを見る量です。

ここで、平均だけを見ると危険です。難しいインスタンスが少数混ざると、平均TTSは大きく動きます。逆に中央値だけを見ると、重い裾を持つ分布を見落とします。Rønnowらの量子スピードアップ論文で強調されたように、すべてのインスタンスをまとめた比較、同じインスタンス同士の比較、最適なアニーリング時間を各インスタンスで選んだ比較は、別々の主張になります。

代表的なベンチマークを、手法の読み方に沿って整理します。

研究	問題設定	主な読み方
Kirkpatrick--Gelatt--Vecchi 1983	配置最適化とスピングラス的エネルギー地形	SAを統計力学の最適化手法として定式化
Kadowaki--Nishimori 1998	小規模横磁場Ising模型の時間発展	熱ゆらぎではなく量子ゆらぎで探索するQAの原型
Martonak--Santoro--Tosatti 2002, Santoroら 2002	2次元ランダムIsing模型	PIMC-QAとSAの残留エネルギー比較、Landau-Zener列としての理解
Boixoら 2013, 2014	D-Wave One/Twoのランダム・設計Ising問題	実機統計がSAだけでは説明しにくいことを調べる
Rønnowら 2014	最大503量子ビットのランダムChimeraスピングラス	量子スピードアップの定義と、ベンチマーク設計の注意点
Shin--Smith--Smolin--Vazirani 2014	D-Wave Oneの入出力統計	SVMCが一部データと高相関を示し、量子性判定の基準を厳しくした
Heim--Rønnow--Isakov--Troyer 2015	2次元Isingスピングラス	離散時間SQAの見かけの優位性と連続時間極限の違いを示した
Denchevら 2016, Boixoら 2016	弱強クラスタと有限幅トンネル障壁	D-Wave、SQA、SVMC、SAを比較し、多体トンネルと熱活性化を切り分ける
Mandràら 2016	弱強クラスタの詳細比較	専用古典ヒューリスティックを入れると、スピードアップ主張が限定されることを示した
Crosson--Harrow 2016	spike問題	QAがSAに勝つ領域でもSQAが多項式時間で動く場合を証明
Albash--Lidar 2018	planted frustrated loops とtunneling gadget	D-Wave 2000Q、SA、SVMC、SQAを比較し、SAへのスケーリング優位とSQAの強さを示した

この表の読みどころは、結論が単純に一方向へ進んでいないことです。初期のPIMC-QAはSAより良く見えました。D-Wave実機の統計は、単純なSAだけでは説明しにくい面を示しました。一方でSVMCは、一部の実機統計を古典ベクトル模型でも説明できることを示しました。SQAは、量子トンネルに似た障壁通過を古典計算機上でかなり強く模擬できます。さらに専用古典ヒューリスティックを入れると、問題によっては量子優位の主張は弱まります。

したがって、ベンチマークで問うべきことは「量子か古典か」だけではありません。少なくとも次の問いに分ける必要があります。

1. SAのような単純な熱的探索よりよいか。
2. SVMCのような半古典的模型では説明できない特徴があるか。
3. SQAのような量子インスパイアード古典アルゴリズムよりよいか。
4. 問題構造を使う専用古典ソルバよりよいか。
5. 実機の物理時間、プログラム時間、読み出し時間、埋め込み時間をどう数えるか。
6. 成功確率、残留エネルギー、サンプリング分布、TTSのどれを目的にするか。

この分解をしておくと、結果が否定的に見える論文も役に立ちます。たとえば、SVMCで説明できるデータは「量子性がない」という結論ではなく、「この観測量だけでは量子性を示すには足りない」という情報です。SQAが非常に強いベンチマークは「実機が不要」という結論ではなく、「stoquasticなトンネル経路は古典計算機でも相当よく模擬できるので、実機はそれを越える物理、サイズ、温度、制御性を見せる必要がある」という情報です。

8. 手法を選ぶときの読み方

これらの手法は、2日目のSA/SQA、p-spin模型の平均場解析、適応的量子モンテカルロ法を読むときの物差しになります。

SAでは、状態は常に古典スピン配位です。p-spin模型なら、自由エネルギーの二つの谷の間を熱で越えます。障壁が高ければ遅いです。

SQAでは、状態は虚時間方向を持つ古典スピン配位です。p-spin模型で言えば、単一の磁化 $m$ だけでなく、虚時間方向に伸びたworldlineが障壁をすり抜ける経路を持ちます。stoquasticな横磁場なら、この重みは正で、Monte Carloで扱えます。

SVMCでは、状態は古典的な角度変数 ${\theta }_i$ です。横磁場によってベクトルが連続的に倒れるため、見た目は量子アニーリングに似ています。しかし、多体波動関数の重ね合わせではなく、積状態のエネルギー面を熱的に動いているだけです。だからこそ、SVMCで再現できる現象と再現できない現象を分けることが重要になります。

適応的量子モンテカルロ法は、この流れのさらに先にあります。通常のSQAはstoquastic横磁場に強いです。一方、反強磁性XX項のようなnon-stoquastic項を入れると符号問題が出ます。そこで、横磁化を補助変数として取り出し、有効横磁場を自己無撞着に更新します。つまり「QMCで扱える形へどこまで戻せるか」を考える方法です。

量子アニーリングのシミュレーション手法は、単なる代用品ではありません。実機の性能を測る物差しであり、どの物理を使えば難しい障壁を抜けられるかを調べる実験台でもあります。SA、SQA、連続時間QMC、SVMCを同じ土俵に置くと、量子アニーリングの主張はかなり鍛えられます。主張が鍛えられるほど、残った差分は研究する価値の高い差分になります。

参考文献: [1], [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12], [13]。

参考文献

1. S. Kirkpatrick, C. D. Gelatt, M. P. Vecchi, “Optimization by Simulated Annealing,” Science 220, 671–680 (1983). 出典

2. Tadashi Kadowaki, Hidetoshi Nishimori, “Quantum annealing in the transverse Ising model,” Physical Review E 58, 5355–5363 (1998). 出典

3. Roman Martoňák, Giuseppe E. Santoro, Erio Tosatti, “Quantum annealing by the path-integral Monte Carlo method: The two-dimensional random Ising model,” Physical Review B 66, 094203 (2002). 出典

4. Giuseppe E. Santoro, Roman Martoňák, Erio Tosatti, Roberto Car, “Theory of Quantum Annealing of an Ising Spin Glass,” Science 295, 2427–2430 (2002). 出典

5. Troels F. Rønnow, Zhihui Wang, Joshua Job, Sergio Boixo, Sergei V. Isakov, David Wecker, et al., “Defining and detecting quantum speedup,” Science 345, 420–424 (2014). 出典

6. Seung Woo Shin, Graeme Smith, John A. Smolin, Umesh Vazirani, “How "Quantum" is the D-Wave Machine?,” arXiv:1401.7087 (2014). 出典

7. Bettina Heim, Troels F. Rønnow, Sergei V. Isakov, Matthias Troyer, “Quantum versus classical annealing of Ising spin glasses,” Science 348, 215–217 (2015). 出典

8. Vasil S. Denchev, Sergio Boixo, Sergei V. Isakov, Nan Ding, Ryan Babbush, Vadim Smelyanskiy, et al., “What is the Computational Value of Finite-Range Tunneling?,” Physical Review X 6, 031015 (2016). 出典

9. Sergio Boixo, Vadim N. Smelyanskiy, Alireza Shabani, Sergei V. Isakov, Mark Dykman, Vasil S. Denchev, et al., “Computational multiqubit tunnelling in programmable quantum annealers,” Nature Communications 7, 10327 (2016). 出典

10. Sergei V. Isakov, Guglielmo Mazzola, Vadim N. Smelyanskiy, Zhang Jiang, Sergio Boixo, Hartmut Neven, et al., “Understanding Quantum Tunneling through Quantum Monte Carlo Simulations,” Physical Review Letters 117, 180402 (2016). 出典

11. Salvatore Mandrà, Zheng Zhu, Wenlong Wang, Alejandro Perdomo-Ortiz, Helmut G. Katzgraber, “Strengths and weaknesses of weak-strong cluster problems: A detailed overview of state-of-the-art classical heuristics versus quantum approaches,” Physical Review A 94, 022337 (2016). 出典

12. Elizabeth Crosson, Aram W. Harrow, “Simulated Quantum Annealing Can Be Exponentially Faster Than Classical Simulated Annealing,” 2016 IEEE 57th Annual Symposium on Foundations of Computer Science (FOCS), 714–723 (2016). 出典

13. Tameem Albash, Daniel A. Lidar, “Demonstration of a Scaling Advantage for a Quantum Annealer over Simulated Annealing,” Physical Review X 8, 031016 (2018). 出典