マルチカノニカル法

マルチカノニカル法（英: multicanonical algorithm）とは、統計学および物理学において、マルコフ連鎖モンテカルロ法におけるサンプリング法の一つで、メトロポリス・ヘイスティングス法により複数の局所安定点を持つ被積分関数を積分するときに用いられる。サンプリングを状態密度の逆数に従って行う^[1]。状態密度はあらかじめワン・ランダウ法などの他の方法により計算しておく。^[2]イジングモデルやスピングラスなどのスピン系において重要なアルゴリズムである^[1]^[3]^[4]。マルチカノニカルアンサンブル（英: multicanonical ensemble）、マルチカノニカルサンプリング（英: multicanonical sampling）、フラットヒストグラム法（英: flat histogram）とも。

背景

たとえばスピン系などの自由度の大きい系をシミュレートする場合、モンテカルロ積分（英語版）が必要とされる。この積分法では、重点サンプリング法（英語版）、とくにメトロポリス法が非常に重要である^[3]。しかし、メトロポリス法では $β$ を逆温度として $\exp(-\beta E)$ に従ってサンプリングを行なうため、エネルギー障壁 $\Delta E$ を乗り越えることはその大きさに対して指数関数的に困難となる^[1]。ポッツモデル（英語版）のような複数の局所安定点を持つ系では、このメトロポリス法はある特定の局所安定点に囚われてしまい、取り扱いが困難となる^[3]。このため、他のサンプリング法によるアプローチが必要となる。

概要

マルチカノニカル法では、メトロポリス・ヘイスティングス法を用いるが、メトロポリス法におけるサンプリング分布 $\exp(-\beta E)$ の代わりに、状態密度の逆数を用いる^[1]。このように選ぶことで、平均的な意味において、エネルギー毎に等しい状態数をサンプリングできる。言い換えれば、エネルギーに対して「フラットなヒストグラム」を得ることができる。このことから、エネルギー障壁を乗り越えることが困難ではなくなる。メトロポリス法と比べて優れている点としてもう一つ、サンプリングが系の温度によらないことから、全ての温度における熱力学的変数を推計することができる（「マルチカノニカル」とはここから来ている）。このことは一次相転移を研究する上において大きな利点となる^[1]。

マルチカノニカル法の実行上において最大の問題は、状態密度を予め知る必要があることである^[2]^[3]。マルチカノニカルサンプリングにおいて重要な貢献の一つとして、ワン・ランダウ法が挙げられる。このアルゴリズムは、収束するまで状態密度を計算しつづけることで、漸近的にマルチカノニカル集合を得ることができる^[2]。

マルチカノニカル法の応用範囲は物理系のみに留まらない。コスト関数 $F$ が定義できる抽象的な系にも適用することができる。 $F$ に対する状態密度を用いることにより、多次元空間上の積分や最小値探索に用いることができる^[5]。

動機

ある系とその位相空間 $\Omega$ （系のコンフィギュレーション ${\boldsymbol {r}}\in \Omega$ により張られる空間）とその位相空間上定義される「コスト」関数 $F$ 、そのスペクトル $\Gamma$ : $F(\Omega )=\Gamma =[\Gamma _{\min },\Gamma _{\max }]$ を考える。

ある量の位相空間上の平均値 $\langle Q\rangle$ を計算するには、以下のような積分を実行する必要がある。

\langle Q\rangle ={\frac {1}{V}}\int _{\Omega }Q({\boldsymbol {r}})p({\boldsymbol {r}})\,\mathrm {d} {\boldsymbol {r}}

ここで $p({\boldsymbol {r}})$ は体積あたりの重み、

p({\boldsymbol {r}})={\frac {1}{V}}\int _{\Omega }p({\boldsymbol {r}})\,\mathrm {d} {\boldsymbol {r}}

である。 $F$ に対する状態密度は以下のように与えられる。

\rho (f)={\frac {1}{V}}\int _{\Omega }\delta (f-F({\boldsymbol {r}}))\,\mathrm {d} {\boldsymbol {r}}

これは、もし $Q$ と $p$ の両方が特定の状態には依存せず、以下のようにその状態における $F$ の値 $F({\boldsymbol {r}})=F_{\boldsymbol {r}}$ にのみ依存するならば、

Q({\boldsymbol {r}})=Q(F_{\boldsymbol {r}}),p({\boldsymbol {r}})=p(F_{\boldsymbol {r}}),

$\langle Q\rangle$ の定義式をデルタ関数と f についての積分を追加することにより以下のように書き換えることができる。

{\begin{aligned}\langle Q\rangle &=\int _{\Gamma _{\min }}^{\Gamma _{\max }}\int _{\Omega }Q(F_{\boldsymbol {r}})p(F_{\boldsymbol {r}})\delta (f-F_{\boldsymbol {r}})\,\mathrm {d} {\boldsymbol {r}}\,df\\&=\int _{\Gamma _{\min }}^{\Gamma _{\max }}Q(f)p(f)\int _{\Omega }\delta (f-F_{\boldsymbol {r}})\,\mathrm {d} {\boldsymbol {r}}\,df\\&=\int _{\Gamma _{\min }}^{\Gamma _{\max }}Q(f)p(f)\rho (f)\,\mathrm {d} f\\\end{aligned}}

つまり、状態密度を知っていれば、それは $\Omega$ 上の状態数の $\Gamma$ への射影であるから、 $F$ の値についての一次元積分により、位相空間上の多次元積分そすることなく平均値が計算できることになる。

自由度の高い系については状態数が大きくなるため、解析的な表式を得ることは難しく、 $\langle Q\rangle$ の計算も困難となる。この問題は多次元積分であるから、モンテカルロ積分がよく用いられる。最も単純な表式では、一様分布に従う $N$ 個のランダムな状態 ${\boldsymbol {r}}_{i}\in \Omega$ を生成し、推定量

{\overline {Q}}_{N}={\frac {1}{N}}\sum _{i=0}^{N}Q({\boldsymbol {r}}_{i})p({\boldsymbol {r}}_{i})

を用いて $\langle Q\rangle$ を計算することができる。 ${\overline {Q}}_{N}$ は強い大数の法則により、ほぼ確実に $\langle Q\rangle$ に収束する。

\lim _{N\rightarrow \infty }{\overline {Q}}_{N}=\langle Q\rangle .

収束に関するよくある問題のひとつとして、 $Q$ の分散が非常に大きい場合、妥当な結果を得るまでの計算コストが大きくなってしまうという問題がある。

収束を速くするための方法として、メトロポリス・ヘイスティングス法がある。モンテカルロ法は一般的に重点サンプリングを用いて推定量の収束を加速する。つまり、 ${\overline {Q}}_{N}$ を適当な分布 $P({\boldsymbol {r}})$ （上述の $p$ とは異なることに注意）を用いて以下のように定義しなおす。

{\overline {Q}}_{N}={\frac {1}{X}}\sum _{i=0}^{N}Q({\boldsymbol {r}}_{i})P^{-1}({\boldsymbol {r}}_{i})p({\boldsymbol {r}}_{i})

ここで、 $X=\sum _{i=0}^{N}P^{-1}({\boldsymbol {r}}_{i})$ である。

当然、 $P$ が一様分布のとき、推定量は一様サンプリングを用いた推定量と一致する。

重要な例として、コスト関数をエネルギーとみて P をある温度に対するボルツマン因子と一致させる選び方がある。

P({\boldsymbol {r}})=p_{\text{Boltzmann}}({\boldsymbol {r}})={\frac {e^{-\beta F({\boldsymbol {r}})}}{\int _{\Omega }\,\mathrm {d} {\boldsymbol {r}}e^{-\beta F({\boldsymbol {r}})}}}

このように選ぶと、コスト関数が低いような状態がよりサンプルされやすくなる。

歴史的には、もともとのアイデア（英語版）^[6]ではメトロポリス・ヘイスティングス法は熱浴に接続された系における平均値を計算するために用いていたため、状態の重みはボルツマン因子と一致していた。このような系では、上述した状態の選び方により物理系の研究がずいぶん進展した^[3]。

しかし、サンプリング分布を重み関数と一致させる必要は必ずしもない。重み関数と一致しないサンプリング分布を用いる理由の一つは、メトロポリス・ヘイスティングス法はコスト関数が複数の極小値を持つ場合に収束しないことである^[1]。これは、このアルゴリズムが局所的ステップによるランダムウォークを用いているためである。つまり、ランダムウォークは通常、エネルギーの差分がオーダーひとつ分であるようなステップを実行する。このことから、このアルゴリズムにおいてある極小値近傍の領域から抜け出すための計算コストはコスト関数の極小値の値に対して指数関数的に増加する。つまり、極小値が深ければ深いほど、このアルゴリズムはそこに長く留まってしまい、そこから抜け出すのは（深さに対して指数関数的に）難しくなる。

これがマルチカノニカル法を導入する動機である。サンプル時にはコスト関数の極小値の存在を「見えなく」してやることによって、最小値近傍に囚われてしまうのを防ぐというアイデアである。

マルチカノニアルアンサブル

マルチカノニカルアンサンブルとは、次のような重み関数を用いて重みづけサンプリングして得られるような統計集団である。

P({\boldsymbol {r}})={\frac {1}{\rho (F({\boldsymbol {r}}))}}

ここで、

\rho (f)={\frac {1}{V}}\int _{\Omega }\delta (F({\boldsymbol {r}})-f)\,d{\boldsymbol {r}}

はコスト関数に対する状態密度である（ $V$ は位相空間体積）。このような重み関数を選ぶことで、サンプリングされた状態に対応するコスト関数の値が $f$ になるような確率密度関数は以下のようになる。

P(f)={\frac {1}{f_{\max }-f_{\min }}}\int _{\Omega }\delta (f-F({\boldsymbol {r}}))P({\boldsymbol {r}})\,\mathrm {d} {\boldsymbol {r}}={\frac {1}{f_{\max }-f_{\min }}}{\frac {1}{V}}\rho (f)\int _{\Omega }\delta (f-F({\boldsymbol {r}}))\,\mathrm {d} {\boldsymbol {r}}={\frac {1}{f_{\max }-f_{\min }}}={\text{constant}}

このため、「フラットヒストグラム」法と呼ばれることもある。つまり、コスト関数のあらゆる値が均一にサンプリングされ、障壁は存在しなくなる。熱浴に接続された系の場合、次のような重要な利点も存在する。すなわち、サンプリングが温度に依存しないため、一回のシミュレーションにより全ての温度について研究することができる。このゆえに、マルチカノニカル（複数の温度）の名がある。

トンネリング時間と致命的スローダウン

どんなモンテカルロ法もそうであるように、 $P({\boldsymbol {r}})$ に従ってサンプリングされるサンプルには相関がある。相関の尺度として「トンネル時間」がよく使われる。トンネリング時間は、シミュレーションが $F$ のスペクトルの極小値と極大値の間を往復するまでにかかるマルコフ連鎖のステップ数と定義される。トンネリング時間を用いる動機の一つは、スペクトルを横断するならば状態密度の極大領域を通過するため、プロセスが脱相関されることである。また、往復させることでシミュレーションがスペクトルの全域を通過することを保証している。

$F$ に対するヒストグラムは平坦であるため、マルチカノニカル法は $F$ の値という一次元の線上における拡散過程（つまりランダムウォーク）とみることができる。このようなプロセスにおいては、詳細釣り合い条件によりドリフト（英語版）がないことが保証される^[7]。このことから、局所的なダイナミクスにおいてはトンネリング時間は拡散プロセスと同じようにスケールすることになり、したがってトンネリング時間はスペクトルのサイズ $N$ の二乗に比例してスケールする。

\tau _{tt}\propto N^{2}

しかし、イジングモデルを初めとするいくつかの系では、スケーリングが致命的にスローダウンし、 $N^{2+z}$ （ここで $z>0$ は系に依存する数）に比例してスケールするようになる^[4]。

スケーリングを改善し、致命的スローダウンを克服するため、非局所的ダイナミクスが開発されている^[8]（ウォルフのアルゴリズムを参照）。しかし、イジングモデルのようなスピン系において致命的スローダウンに陥いらないような局所的ダイナミクスが存在するのかどうかという疑問はいまだに答えの出ない問題である。

出典

^ ^a ^b ^c ^d ^e ^f Berg, B.; Neuhaus, T. (1992). “Multicanonical ensemble: A new approach to simulate first-order phase transitions”. Physical Review Letters 68 (1): 9–12. arXiv:hep-lat/9202004. Bibcode: 1992PhRvL..68....9B. doi:10.1103/PhysRevLett.68.9. PMID 10045099.
^ ^a ^b ^c Wang, F.; Landau, D. (2001). “Efficient, Multiple-Range Random Walk Algorithm to Calculate the Density of States”. Physical Review Letters 86 (10): 2050–2053. arXiv:cond-mat/0011174. Bibcode: 2001PhRvL..86.2050W. doi:10.1103/PhysRevLett.86.2050. PMID 11289852.
^ ^a ^b ^c ^d ^e Newmann, M E J; Barkema, G T (2002). Monte Carlo Methods in Statistical Physics. USA: Oxford University Press. ISBN 0198517971
^ ^a ^b Dayal, P.; Trebst, S.; Wessel, S.; Würtz, D.; Troyer, M.; Sabhapandit, S.; Coppersmith, S. (2004). “Performance Limitations of Flat-Histogram Methods”. Physical Review Letters 92 (9). arXiv:cond-mat/0306108. Bibcode: 2004PhRvL..92i7201D. doi:10.1103/PhysRevLett.92.097201.
^ Lee, J.; Choi, M. (1994). “Optimization by multicanonical annealing and the traveling salesman problem”. Physical Review E 50 (2): R651. Bibcode: 1994PhRvE..50..651L. doi:10.1103/PhysRevE.50.R651.
^ Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. (1953). “Equation of State Calculations by Fast Computing Machines”. The Journal of Chemical Physics 21 (6): 1087. Bibcode: 1953JChPh..21.1087M. doi:10.1063/1.1699114.
^ Robert, Christian; Casella, George (2004). Monte Carlo statistical methods. Springer. ISBN 978-0-387-21239-5
^ Wolff, U. (1989). “Collective Monte Carlo Updating for Spin Systems”. Physical Review Letters 62 (4): 361–364. Bibcode: 1989PhRvL..62..361W. doi:10.1103/PhysRevLett.62.361. PMID 10040213.

[Berg-1] ^ ^a ^b ^c ^d ^e ^f Berg, B.; Neuhaus, T. (1992). “Multicanonical ensemble: A new approach to simulate first-order phase transitions”. Physical Review Letters 68 (1): 9–12. arXiv:hep-lat/9202004. Bibcode: 1992PhRvL..68....9B. doi:10.1103/PhysRevLett.68.9. PMID 10045099.

[Landau-2] Wang, F.; Landau, D. (2001). “Efficient, Multiple-Range Random Walk Algorithm to Calculate the Density of States”. Physical Review Letters 86 (10): 2050–2053. arXiv:cond-mat/0011174. Bibcode: 2001PhRvL..86.2050W. doi:10.1103/PhysRevLett.86.2050. PMID 11289852.

[Newmann-3] Newmann, M E J; Barkema, G T (2002). Monte Carlo Methods in Statistical Physics. USA: Oxford University Press. ISBN 0198517971

[Dayal-4] Dayal, P.; Trebst, S.; Wessel, S.; Würtz, D.; Troyer, M.; Sabhapandit, S.; Coppersmith, S. (2004). “Performance Limitations of Flat-Histogram Methods”. Physical Review Letters 92 (9). arXiv:cond-mat/0306108. Bibcode: 2004PhRvL..92i7201D. doi:10.1103/PhysRevLett.92.097201.

[Lee-5] Lee, J.; Choi, M. (1994). “Optimization by multicanonical annealing and the traveling salesman problem”. Physical Review E 50 (2): R651. Bibcode: 1994PhRvE..50..651L. doi:10.1103/PhysRevE.50.R651.

[Metropolis-6] Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. (1953). “Equation of State Calculations by Fast Computing Machines”. The Journal of Chemical Physics 21 (6): 1087. Bibcode: 1953JChPh..21.1087M. doi:10.1063/1.1699114.

[Robert-7] Robert, Christian; Casella, George (2004). Monte Carlo statistical methods. Springer. ISBN 978-0-387-21239-5

[Wolff-8] Wolff, U. (1989). “Collective Monte Carlo Updating for Spin Systems”. Physical Review Letters 62 (4): 361–364. Bibcode: 1989PhRvL..62..361W. doi:10.1103/PhysRevLett.62.361. PMID 10040213.

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]