物理とか

Index

量子モンテカルロ法の負符号問題

[prev:量子モンテカルロ] [index]

量子モンテカルロ法と負符号問題

量子モンテカルロ法について調べだすと、誰もが必ず出くわすキーワードが

「負符号問題」

だと思う。ほとんどどの資料にも現れる単語なので重要そうなのは間違いないんだが、「負符号問題」とか英語で "sign problem" と検索してもふわっとした説明ばかり出てくる。特に日本語でしっかりと説明した資料をみつけることはできなかった。

そこで自分が調べた限りのことを、このページにまとめようと思う。

量子モンテカルロ法の大雑把なまとめ

負符号問題が現れる量子モンテカルロ法の種類は、前のページで説明した経路積分量子モンテカルロの仲間たち (グリーン関数量子モンテカルロ・拡散モンテカルロ) だ。まず経路積分量子モンテカルロ法について復習しよう。

基底状態を求めたいハミルトニアンを \(H\)、その基底状態を \(\ket{\psi_g}\) とする。経路積分量子モンテカルロは、適当な初期状態 \(\ket{\psi_0}\) に対して、虚時間発展演算子 \(e^{-\beta H}\) (\(\beta\gt 0\)) を作用させると、十分大きな \(\beta\) で、 \begin{align} \ket{\psi_g} \approx \frac{e^{-\beta H}\ket{\psi_0}}{\|e^{-\beta H}\ket{\psi_0}\|} = \frac{e^{-\beta H}\ket{\psi_0}}{\sqrt{\bra{\psi_0}e^{-2\beta H}\ket{\psi_0}}} \end{align} となることを利用する手法だった。\(H\) は指数関数的に大きな行列であり、\(e^{-\beta H}\) という行列をそのまま計算することは難しいので、なんらかの近似をしなければならない。前のページ では \(H=\sum_i H_i\) と小さなパーツに分解して、小さな \(H_i\) については \(e^{-\beta H_i}\) が計算できるという仮定のもとに、Trotter 展開を使って近似するアプローチを取った。今回もおなじやり方を使ってもいいのだが、説明を簡単にするために Trotter 展開ではなくて Taylor 展開 \begin{align} e^{-\Delta\beta H} = I - \Delta\beta H + O(\Delta\beta^2) \end{align} (\(\Delta\beta\)は十分小さい実数) を使う。すると任意の \(\beta\) に対する虚時間発展は、十分大きな分割数 \(n\) を使って、 \begin{align} e^{-\beta H} &= \left(e^{-\beta H/n}\right)^n\\ &\approx \left(I-\frac{\beta H}{n}\right)^n \end{align} と近似できる。したがって \begin{align} \ket{\psi_g} &\approx \left(I-\frac{\beta H}{n}\right)^n \ket{\psi_0}\\ &\equiv \ket{\tilde{\psi}_g} \end{align} とできる。これを前のページと同じように、基底 \(\ket{j}\) で展開すれば、 \begin{align} \ket{\tilde{\psi}_g} &= \sum_{j_1,j_2,\cdots,j_{n+1}} \ket{j_{n+1}}\bra{j_{n+1}}\left(I-\frac{\beta H}{n}\right)\ket{j_{n}}\bra{j_{n}}\cdots \ket{j_2}\bra{j_2}\left(I-\frac{\beta H}{n}\right)\ket{j_1}\braket{j_1}{\psi_0} \end{align} となる。つまり \begin{align} G_{ij} = \bra{i}\left(I-\frac{\beta H}{n}\right)\ket{j} \end{align} という行列を、初期状態 \(\ket{\psi_0}\) に掛け続ければ、基底状態 \(\ket{\psi_g}\) が近似的に得られることがわかる。つまり \begin{align} \ket{\tilde{\psi}_g} &= \sum_{j_1,j_2,\cdots,j_{n+1}} G_{j_{n+1}j_n}\cdots G_{j_2j_1}\braket{j_1}{\psi_0}\ket{j_{n+1}} \end{align} である。今回は簡単のため \(G_{ij}\) は実数であると仮定する。

さて、量子状態は上の式によって得られるが、物理量を得るためにはもう少し計算が必要である。物理量 \(A\) の基底状態に関する期待値 \(\bra{\psi_g}A\ket{\psi_g}\) を計算するには、\( \ket{\psi_g} \approx \frac{e^{-\beta H}\ket{\psi_0}}{\|e^{-\beta H}\ket{\psi_0}\|}\) だったので、 \begin{align} \bra{\psi_g}A\ket{\psi_g} \approx \frac{\bra{\psi_0}e^{-\beta H}Ae^{-\beta H}\ket{\psi_0}}{\bra{\psi_0}e^{-2\beta H}\ket{\psi_0}} \end{align} を計算しなくてはならない。いま、\(e^{-\beta H}\ket{\psi_0}\) は上の \(\tilde{\ket{\psi_g}}\) で近似することにしたので、 \begin{align} \bra{\psi_g}A\ket{\psi_g} \approx \frac{\bra{\tilde{\psi}_g}A\ket{\tilde{\psi}_g}}{\braket{\tilde{\psi}_g}{\tilde{\psi}_g}} \end{align} である。

負符号問題の定式化に向けて

さて、量子モンテカルロでは、\(\ket{\tilde{\psi}_g}\) を確率過程によって作り出すというアプローチなので、当然のことながら、上式の分子 \(\bra{\tilde{\psi}_g}A\ket{\tilde{\psi}_g}\) と分母 \(\braket{\tilde{\psi}_g}{\tilde{\psi}_g}\) もサンプリングによって計算することになる。

このとき現れるのが負符号問題で、以下のような問題である。

負符号問題:


\(w_{ij}\) に負の値が現れるとき、すなわち \(G_{ij} = \bra{i}\left(I-\beta H/n\right)\ket{j}\) が負の値が現れるときに、これらの量、特に \(\braket{\tilde{\psi}_g}{\tilde{\psi}_g}\) をサンプリングによって計算しようとすると、その分散が非常に大きくなってしまう現象のことを、負符号問題と呼ぶ。
以下ではこれがどういうことか、具体的な式によって説明していく。\(\braket{\tilde{\psi}_g}{\tilde{\psi}_g}\) の計算を見ていこう。\(\braket{\tilde{\psi}_g}{\tilde{\psi}_g}\) を計算するためには、 \begin{align} \braket{\tilde{\psi}_g}{\tilde{\psi}_g} = \sum_{j_1,j_2,\cdots,j_{2n+2}} G_{j_{2n+2}j_{2n+1}}\cdots G_{j_3j_2}G_{j_2j_1}\braket{\psi_0}{j_{2n+2}}\braket{j_1}{\psi_0} \end{align} 話を簡単にするために、初期状態を最大混合状態 \(I\) であることにして、\(\braket{\psi_0}{j_{2n+2}}\braket{j_1}{\psi_0} \to \delta_{j_1j_{2n+2}}\) としよう。 すると \begin{align} \braket{\tilde{\psi}_g}{\tilde{\psi}_g} &= \sum_{j_1,j_2,\cdots,j_{2n+2}} G_{j_{2n+2}j_{2n+1}}\cdots G_{j_3j_2}G_{j_2j_1}\delta_{j_1j_{2n+2}} \\ &= \sum_{j_1,j_2,\cdots,j_{2n+1}} G_{j_1j_{2n+1}}\cdots G_{j_3j_2}G_{j_2j_1} \end{align} となる。さらに式の見た目を怖くなくするために \begin{align} \b{j} &= (j_1, j_2, \cdots, j_{2n+1}) \\ G_{\b{j}} &:= G_{j_{1}j_{2n+1}}\cdots G_{j_3j_2}G_{j_2j_1} \end{align} とおこう。すると \begin{align} \braket{\tilde{\psi}_g}{\tilde{\psi}_g} &= \sum_{\b{j}} G_{\b{j}} \end{align} である。 量子モンテカルロでは、上の和 \(\sum_{\b{j}} G_{\b{j}}\) をサンプリングによって近似する。

最適な推定

\(Z:=\sum_{\b{j}} G_{\b{j}}\) をサンプリングによって推定するとき、最適なサンプリング方法はなんだろうか？これが分かれば、\(Z\) をモンテカルロで推定するときに到達可能な原理的限界がわかる。

そこでまず、\(\sum_{\b{j}} G_{\b{j}}\) サンプリングによって推定する、ということの意味を明確にしておこう。サンプリングによって \(\sum_{\b{j}} G_{\b{j}}\) を推定するということは、適当な確率分布 \(p_{\b{j}}\) について確率変数 \(\b{j}\) をサンプリングし、 \begin{align} \mathbb{E}_{\b{j}\sim p_{\b{j}}}[\hat{Z}_{\b{j}}] = \sum_{\b{j}} G_{\b{j}} \end{align} となるような \(\b{j}\) に依存する確率変数 \(\hat{Z}_{\b{j}}\) の平均をとるということだ。ここで \(\hat{Z}_{\b{j}}\) は \begin{align} \sum_{\b{j}} p_{\b{j}}\hat{Z}_{\b{j}} = \sum_{\b{j}} G_{\b{j}} \end{align} を満たさないといけないので、 \begin{align} \hat{Z}_{\b{j}} = \frac{G_{\b{j}}}{p_{\b{j}}} \end{align} と取ればよいだろう。

具体的なアルゴリズムでいうと、\(p_{\b{j}}\) にしたがって \(S\) 個のサンプル \(\{\b{j}_s\}_{s=1}^S\) を取ってきて、これらをもとに、 \begin{align} \bar{Z} = \frac{1}{S} \sum_{s} \hat{Z}_{\b{j}_s} \end{align} を計算すれば、\(S\) が十分大きいとき \(\bar{Z}\) は \(Z\) の良い近似となるはずである。このとき、「最適なサンプリング」とは「近似の精度」が最も高くなる確率分布 \(p_\b{j}\) である。さらに「近似の精度」は \(\bar{Z}\) の分散に対応すると考えられるが、\(\bar{Z}\) の分散には以下の確率論でよく知られた事実が使える。
独立同分布な確率変数 \(\{X_s\}_{s=1}^S\) を取ってきて、 \begin{align} \bar{X} = \frac{1}{S}\sum_{s=1}^S X_s \end{align} を作ると、分散は \(\mathrm{Var}[\bar{X}] = \mathrm{Var}[X]/S\) となる。
つまり、\(\hat{Z}_{\b{j}}\) の分散が最も小さくなるとき、\(\bar{Z}\) の分散も最も小さくなる。したがって、「最適なサンプリング」とは \(\hat{Z}_{\b{j}}\) の分散を最も小さくするものになる。

実は次のことが言える。
\(\hat{Z}_{\b{j}}\) の分散を最も小さくする \(p_{\b{j}}\) は \begin{align} p_{\b{j}} = \frac{|G_{\b{j}}|}{\sum_{\b{j}}|G_{\b{j}}|}. \end{align} またこのとき \(\hat{Z}_{\b{j}} = \mathrm{sign}(G_{\b{j}})\sum_{\b{j}}|G_{\b{j}}|\). ただし \(\mathrm{sign}(G_{\b{j}})\) は \(G_{\b{j}}\) の符号を表し、\(\mathrm{sign}(G_{\b{j}}) =G_{\b{j}}/|G_{\b{j}}|\).
このことは \(\hat{Z}_{\b{j}}\) の分散を直接計算することで示せる。 \begin{align} \mathrm{Var}[\hat{Z}_{\b{j}}] &= \mathbb{E}[\hat{Z}_{\b{j}}^2] - \mathbb{E}[\hat{Z}_{\b{j}}]^2 \\ &= \sum_{\b{j}} \frac{G_{\b{j}}^2}{p_{\b{j}}} - \left(\sum_{\b{j}} G_{\b{j}}\right)^2 \end{align} となるので、\(\sum_{\b{j}} \frac{G_{\b{j}}^2}{p_{\b{j}}}\) が最小となれば分散が最小となる。この最小化には "シュワルツの不等式の応用公式|高校数学の美しい物語" のページの最初の公式がそのまま使える。結果だけ使わせていただくと、 \begin{align} \sum_{\b{j}} \frac{G_{\b{j}}^2}{p_{\b{j}}} \geq \left(\sum_{\b{j}} |G_{\b{j}}|\right)^2 \end{align} で、最小値を達成する等号成立条件は \(p_{\b{j}} \propto |G_{\b{j}}| \) である。これで上の最適条件が示された。

ここの部分の参考文献: arXiv:1503.07525

負符号問題の定式化

最適な \(p_{\b{j}}\) を使ったとき、 \begin{align} \mathrm{Var}[\hat{Z}_{\b{j}}] &= \left(\sum_{\b{j}} |G_{\b{j}}|\right)^2 - \left(\sum_{\b{j}} G_{\b{j}}\right)^2 \end{align} となる。したがって \(Z=\sum_{\b{j}}\) の値に対する相対誤差は \begin{align} \frac{\mathrm{Var}[\hat{Z}_{\b{j}}]}{Z^2} &= \frac{\left(\sum_{\b{j}} |G_{\b{j}}|\right)^2}{\left(\sum_{\b{j}} G_{\b{j}}\right)^2} - 1 \\ &= \left(\frac{\sum_{\b{j}} |G_{\b{j}}|}{\sum_{\b{j}} G_{\b{j}}}\right)^2 - 1 \end{align} となる。

上の式に現れた量 \begin{align} \frac{\sum_{\b{j}} |G_{\b{j}}|}{\sum_{\b{j}} G_{\b{j}}} \end{align} の逆数 \begin{align} \langle\mathrm{sign}_{G}\rangle := \frac{\sum_{\b{j}} G_{\b{j}}}{\sum_{\b{j}} |G_{\b{j}}|} \end{align} は、\(G_{\b{j}}\) の「平均的な符号」 と考えられることがわかると思う。\(G_{\b{j}}\) が全て正のとき、\(\langle\mathrm{sign}_{G}\rangle=+1\) となり、全て負のとき、\(\langle\mathrm{sign}_{G}\rangle = -1\) になる。また、\(G_{\b{j}}\) に正の成分と負の成分が同じくらい出てきて打ち消しあってしまうときには、\(\langle\mathrm{sign}_{G}\rangle\approx 0\) になり、分散が非常に大きくなる。ここまでの議論をまとめると、負符号問題とは以下のことであるといえる。
最適な確率分布 \(p_{\b{j}} = \frac{|G_{\b{j}}|}{\sum_{\b{j}}|G_{\b{j}}|}\) からの \(S\) 個のサンプルを使って、\(Z=\sum_{\b{j}} G_{\b{j}}\) を推定しようとしたとき、その推定値 \(\bar{Z}\) の分散は \begin{align} \mathrm{Var}[\bar{Z}] &= \frac{\mathrm{Var}[\hat{Z}]}{S} \\ &= \frac{1}{S} \left(\frac{1}{\langle\mathrm{sign}_{G}\rangle^2} -1\right) \end{align} になる。したがって、平均的な符号 \(\langle\mathrm{sign}_{G}\rangle\) が 0 に近くなると、\(Z\) の推定を確率過程を使って行うという方針が、原理的に破綻する。これを量子モンテカルロ法における

負符号問題

と呼ぶ。

負符号問題と問題の難しさ

量子モンテカルロ法の資料を色々読んでいると、だいたいのもので、

「量子モンテカルロ法には負符号問題と呼ばれる問題があって、これをできる限り避けるための手法が色々開発されている」

という趣旨の文章が書いてある。こういう文章を僕は最初、「負符号問題さえ解決すれば、量子モンテカルロ法はうまく行って、色々な量子系の計算ができるのだ」と読んでしまっていたが、上の議論を踏まえると、どうもそういうことではないようだと思った。

負符号問題が無い (=\(G_{\b{j}}\) が全て0以上) からといって、計算がうまくいくわけではない。このことは次のようにしてわかる。\(G_{\b{j}}\) の定義に戻って計算すると、 \begin{align} \sum_{\b{j}}|G_{\b{j}}| &= \sum_{\b{j}} G_{\b{j}}~ (G_{\b{j}}\geq 0 \text{の仮定から}) \\ &= \Tr\left[\left(I-\frac{\beta H}{n}\right)^{2n}\right] \\ &\approx \Tr\left[e^{-2\beta H}\right] \end{align} となる。このとき負符号問題が生じないならば、\(G_{\b{j}}\) が全て 0 以上という条件から、ハミルトニアンには \(I-\frac{\beta H}{n}\) の成分が全て 0 以上であるという条件が課される。この条件は、対角なハミルトニアンを持ってきて、\(\beta/n\) を十分小さく取れば満たせる。例えば対角なハミルトニアンとして、一般のイジングモデル \begin{align} H = \sum_{i} B_i Z_i + \sum_{i,j} J_{ij} Z_i Z_j \end{align} (\(Z_i\)は\(i\)番目のスピンに作用するパウリ演算子のz成分) を取ってくることもできる。ちょっと調べてみると、この \(H\) に関して一般に \(\Tr\left[e^{-2\beta H}\right]\) を計算することは非常に困難であることがわかる。例えば arXiv:1409.5627 を見ると、\(\Tr\left[e^{-2\beta H}\right]\) の近似計算は #P-hard であるらしい。つまりこの計算は負符号問題を持たないが、一方で非常に困難である。

ある計算が #P-hard であるとは、その計算ができると、#P という計算量クラスに入っている問題が全て解けるということ。#P は解くのが難しいことで有名な NP をもっと難しくしたような計算量クラスなので、一般に \(H\) が対角であったとしても \(\Tr\left[e^{-2\beta H}\right]\) を計算するのはとんでもなく難しい。

ちょっと待って、最適な分布 \(p_{\b{j}} = \frac{|G_{\b{j}}|}{\sum_{\b{j}}|G_{\b{j}}|}\) を使えば、負符号問題が無いとき分散が 0 になるはずなのになぜこんなことが起きるの？と思う方もいるかも知れない。確かにこのとき \(\langle\mathrm{sign}_{G}\rangle=1\) なので分散は 0 になるはずである。しかしよく考えてみると、\(p_{\b{j}}\) の分母 \(\sum_{\b{j}}|G_{\b{j}}|\) は、負符号問題が無いとき \(\sum_{\b{j}}|G_{\b{j}}|=\sum_{\b{j}}G_{\b{j}}\) となり、これは今求めようとしていた \(\Tr\left[e^{-2\beta H}\right]\) そのものである。つまり実際問題、最適な分布からサンプリングすることは不可能なのだ。実際には、前のページのようにマルコフ連鎖を使うなど、他の手法でサンプリングするほかなく、これによって分散が (難しい問題であればあるほど) 最適値より大きくなってしまう。

ということで、確かに負符号問題は量子系の計算の難しさの一端を示すものではあるものの、「負符号問題を取り除けば問題が全て解ける」ということでは無いことがわかった。

ここの話は arXiv:0408370 [cond-mat] を参考にした。

余談