ワン・ランダウ法のソースを表示

'''ワン・ランダウ法'''（ワン・ランダウほう、{{lang-en-short|Wang and Landau algorithm}}）は、Fugao Wangと{{仮リンク|デビッド・P・ランダウ|en|David_P._Landau}}により発案された<ref name="WangLandau"/>、系の[[状態密度]]を計算するために用いられる[[モンテカルロ法]]のひとつである。このアルゴリズムでは状態密度の計算に必要な、系の取り得る全ての状態のエネルギーを迅速に計算するため、非[[マルコフ連鎖]][[ランダムウォーク]]を行う。[[マルチカノニカル法]]の実行に必要となる、状態密度を計算するために重要である。

ワンランダウ法はコスト（エネルギー）関数により特徴づけられるどのような系にも応用可能である。たとえば、[[数値積分]]<ref name="Belardinelli_Integrals" />や[[タンパク質フォールディング]]<ref name="Ojeda1"/><ref name="Ojeda2"/>への応用が知られている。ワン・ランダウサンプリングは[[メタダイナミクス法]]とも関連する<ref>Junghans, Christoph, Danny Perez, and Thomas Vogel. </ref>。

== 概要 ==
ワン・ランダウ法により、コスト関数で特徴づけられる系の状態密度を得ることができる。非マルコフ連鎖[[確率過程]]を用いて、漸近的に[[マルチカノニカル法|マルチカノニカルアンサンブル]]を得る（サンプル分布が状態密度の逆数となるような[[メトロポリス・ヘイスティングス法]]に漸近する）<ref name="WangLandau"/>。このことの帰結として、[[活性化エネルギー|エネルギー障壁]]を無視したシミュレーションを行うことができる。つまり、このアルゴリズムは（コスト関数的に良いものも悪いものも含めた）系の取り得る全ての状態をメトロポリス法よりも非常に速く網羅することができるのである<ref name="Berg"/>。

== アルゴリズム ==
[[位相空間 (物理学)|位相空間]] <math>\Omega</math> 上に定義される系と、そのコスト関数（例えばエネルギー） {{Mvar|E}} を考える。 {{Mvar|E}} のスペクトルは <math>E\in\Gamma = [E_\min,E_\max]</math> の範囲に限定されるものとし、付随する状態密度 <math>\rho(E) \equiv \exp(S(E))</math> を計算する。 ワン・ランダウ法は離散スペクトルを扱うため<ref name="WangLandau"/>、スペクトル <math>\Gamma</math> をそれぞれ <math>\Delta</math> だけ異なる {{Mvar|N}} 個の離散値で分割する。
: <math> N = \frac{E_\max-E_\min}{\Delta}</math>
このような離散スペクトルに対し、まず次の初期化を行う。
* 全ての区間のミクロカノニカルエントロピーをゼロとする。
*:<math>S(E_i) = 0\ \ i=1,2,...,N</math>
* {{Math|1 = ''f'' = 1}} とする。
* 系の初期配置 <math>\boldsymbol{r}\in\Omega</math> をランダムに選ぶ。
次に、[[マルチカノニカル法|マルチカノニカルアンサンブル]]シミュレーションを行う<ref name="WangLandau"/>。つまり、位相空間上でランダムウォークを行い、[[確率分布]]
:<math>P(\boldsymbol{r}) = 1/\rho(E(\boldsymbol{r})) = \exp (-S(E(\boldsymbol{r})))</math>
を再現するような遷移確率 <math>g(\boldsymbol{r} \rightarrow \boldsymbol{r}')</math> で[[メトロポリス・ヘイスティングス法|メトロポリス]]シミューレーションを行う。このとき、訪れた状態のエネルギーをヒストグラム <math>H(E)</math> に記録する。そしてメトロポリス法と同様に新たな配置の生成と採択/棄却プロセスを実行する。
# 新たな配置 <math>\boldsymbol{r}'\in\Omega</math> を任意の遷移確率分布 <math>g(\boldsymbol{r} \rightarrow \boldsymbol{r}')</math> に従って生成する。<br>
# 新たな配置を次の確率で採択/棄却する。
#: <math>A(\boldsymbol{r}\rightarrow \boldsymbol{r}') = \min\left(1,e^{S - S'}\frac{g(\boldsymbol{r}'\rightarrow \boldsymbol{r})}{g(\boldsymbol{r}\rightarrow \boldsymbol{r}')}\right)</math>
#: ここで、 <math>S = S(E(\boldsymbol{r}))</math> および <math>S' = S(E(\boldsymbol{r}'))</math> とする。
このステップを終えたとき、系の遷移先が <math>E_i</math> だったとすると <math>H(E_i)</math> を 1 だけ増やし、エントロピーを以下のように更新する。
: <math> S(E_i) \leftarrow S(E_i) + f</math>
このステップが本アルゴリズムの中心である。また、このステップが存在するため[[確率過程]]が過程の履歴に依存し、ワン・ランダウ法は非マルコフ連鎖となる。したがって、このエネルギー <math>E_i</math> を持つ状態が次に生成された際にはより棄却されやすくなる。このようにして、本アルゴリズムはスペクトル全域を等しく訪れるよう強制する<ref name="WangLandau"/>。この帰結として、ヒストグラム <math>H(E)</math> はシミュレーションを進めるにつれてより平坦となる。しかし、この平坦さは計算されたエントロピーがどれだけ実際のエントロピーを近似できているかに依存する。そして、近似の良さは {{Mvar|f}} に依存する<ref name="Belardinelli_Saturation"/>。実際のエントロピーに近づくため（すなわちヒストグラムを平坦にするため）に、 {{Mvar|f}} を {{Mvar|M}} ステップごとに次のように減らす。
: <math>f \leftarrow f/2</math>
後に、 {{Mvar|f}} を常に半分にし続けると飽和誤差を招くことが明らかとなった<ref name="Belardinelli_Saturation"/>。この問題を避けるため、{{Mvar|t}} をシミュレーションの総ステップ数として、 {{Mvar|f}} を {{Math|1/''t''}} にする変更がアルゴリズムに加えられた<ref name="Belardinelli_Saturation"/>。

== アルゴリズムのテスト ==
今、次のような[[調和振動子]]ポテンシャルの [[状態密度|DOS]] を求めたいとする。
: <math> E(x) =  x^2 </math>
解析的には DOS は次のように与えられる。
: <math> g(E) = \int \delta (E(x)-E) \, dx= \int \delta (x^2-E) \, dx</math>
最後の積分を解くと位相空間の次元によって次のような結果が得られる。
: <math> g(E) \propto 
\begin{cases} 
E^{-1/2}, \text{for } x \in \mathbb{R}^1 \\
\text{const}, \text{for } x \in \mathbb{R}^2 \\
E^{1/2}, \text{for } x \in \mathbb{R}^3 \\
\end{cases}
</math>
一般的に、多次元調和振動子の DOS は {{Mvar|E}} のべき乗で与えられ、その指数は系の次元の関数となる。

このように単純な調和振動子ポテンシャルに対する状態密度は既に分かっているので、このポテンシャルに対してワン・ランダウ法を行い、 得られた状態密度 <math>\rho(E)</math> と <math>g(E)</math> を比べることによりワン・ランダウ法の精度を確かめることができる。

== サンプルコード ==
次に示すのはPythonで実装されたワン・ランダウ法のサンプルコードである。遷移確率密度は対称であることを仮定している。
: <math>g(\boldsymbol{x}'\rightarrow \boldsymbol{x})=g(\boldsymbol{x}\rightarrow \boldsymbol{x}')</math>
このコードにおいて、対象となる系は "system" により表されるものとしている。<syntaxhighlight lang="python">

currentEnergy = system.randomConfiguration() # a random initial configuration

while (f > epsilon):
    system.proposeConfiguration() # a proposed configuration is proposed
    proposedEnergy = system.proposedEnergy() # the energy of the proposed configuration computed

    if (random() < exp(entropy[currentEnergy]-entropy[proposedEnergy])):
        # if accepted, update the energy and the system:
        currentEnergy = proposedEnergy
        system.acceptProposedConfiguration()
    else:
        # if rejected
        system.rejectProposedConfiguration()
    
    H[currentEnergy] += 1
    entropy[currentEnergy] += f
    
    if (isFlat(H)): # isFlat tests whether the histogram is flat (e.g. 95% flatness)
        H[:] = 0
        f *= 0.5 # refine the f parameter
</syntaxhighlight>

== ワン・ランダウ分子動力学法 ==
ワン・ランダウ法はモンテカルロ法のみならず[[分子動力学法]]にも適用することができる。そのためには、系の[[温度]]を次のように制御する。
: <math> T'(E) \rightarrow \left(\frac{\partial S(E)}{\partial E}\right) T(E)</math>
ここで {{Math|''S''(''E'')}} は系のエントロピー、 {{Math|''T''(''E'')}} は[[ミクロカノニカルアンサンブル|ミクロカノニカル]]温度、 {{Math|''T''&prime;(''E'')}} は分子動力学シミュレーションに実際用いられる「スケールされた」温度である。

== 出典 ==
{{Reflist|refs = <ref name=WangLandau>
{{cite journal
| title = Efficient, Multiple-Range Random Walk Algorithm to Calculate the Density of States
| author = Wang, Fugao  and Landau, D. P.
| journal = Phys. Rev. Lett.
|  volume = 86
|  issue = 10
|  pages = 2050–2053
|date=Mar 2001
| pmid = 11289852
|  doi = 10.1103/PhysRevLett.86.2050
|  publisher = American Physical Society
| bibcode=2001PhRvL..86.2050W
|arxiv = cond-mat/0011174 }}
</ref>
<ref name=Belardinelli_Integrals>
{{cite journal
| title = Analysis of the convergence of the 1∕t and Wang–Landau algorithms in the calculation of multidimensional integrals
| author = R. E. Belardinelli and S. Manzi and V. D. Pereyra
| journal = Phys. Rev. E
|  volume = 78
|  page = 067701
|date=Dec 2008
|  doi = 10.1103/PhysRevE.78.067701
|  publisher = American Physical Society
| issue = 6
|arxiv = 0806.0268 |bibcode = 2008PhRvE..78f7701B }}
</ref>
<ref name=Ojeda1>
{{cite journal
| title = Monte Carlo Simulations of Proteins in Cages: Influence of Confinement on the Stability of Intermediate States
| author = P. Ojeda and M. Garcia and A. Londono and N.Y. Chen
| journal = Biophys. Jour.
|  volume = 96
|  issue = 3
|  pages = 1076–1082
|date=Feb 2009
|  doi = 10.1529/biophysj.107.125369
|  publisher = Biophysical Society
| bibcode=2009BpJ....96.1076O
}}
</ref>
<ref name=Ojeda2>
{{cite journal
| title = Electric Field-Driven Disruption of a Native beta-Sheet Protein Conformation and Generation of alpha-Helix-Structure
| author = P. Ojeda and M. Garcia 
| journal = Biophys. Jour.
|  volume = 99
|  issue = 2
|  pages = 595–599
|date=Jul 2010
|  doi = 10.1016/j.bpj.2010.04.040
|  publisher = Biophysical Society
| bibcode=2009BpJ....96.1076O
| pmid=20643079
| pmc=2905109
}}
</ref>
<ref name=Berg>{{Cite journal | last1 = Berg | first1 = B. | last2 = Neuhaus | first2 = T. | doi = 10.1103/PhysRevLett.68.9 | title = Multicanonical ensemble: A new approach to simulate first-order phase transitions | journal = Physical Review Letters | volume = 68 | issue = 1 | pages = 9–12 | year = 1992 | pmid =  10045099| pmc = |arxiv = hep-lat/9202004 |bibcode = 1992PhRvL..68....9B }}</ref>
<ref name=Belardinelli_Saturation>
{{cite journal
| title   = Wang–Landau algorithm: A theoretical analysis of the saturation of the error
| author  = Belardinelli, R. E. and Pereyra, V. D.
| journal = Jour. Chem. Phys.
| volume  = 127
| issue   = 18
| page   = 184105 
| year    = 2007
| doi     = 10.1063/1.2803061 
|arxiv = cond-mat/0702414 |bibcode = 2007JChPh.127r4105B }}
</ref>}}

{{DEFAULTSORT:わんらんたうほう}}
[[Category:計算物理学]]
[[Category:モンテカルロ法]]
[[Category:統計アルゴリズム]]