Contents

概要

最終的には、古典ハイゼンベルグスピン模型を解くマルコフ連鎖モンテカルロプログラムを書いて、 それを最適化することです。 これを通して、疑似乱数、メモリアクセスの順序、メモリアロケーションの性能への影響を学びます。

マルコフ連鎖モンテカルロ法を知らない人は、川島先生の講義ノート物性物理におけるモンテカルロ法の1.1〜1.8章を参考にしましょう。

2次元イジングスピン模型（正方格子・最近接相互作用のみ)に関しては、黒木玄さんによるコードが公開されています。

以下、モテンカルロ法を実装する際に必要な基礎知識をまとめておきます。

マルコフ連鎖モンテカルロ

一般論

マルコフ連鎖モンテカルロを使うことで、 狙った分布関数に従うような状態の時系列を構成することが出来ます。 詳しいことは、川島先生の講義ノートを参考にしてもらうこととして、コードを読み解くのに必要な知識のみをまとめておきます。

例として、ある分布に従う古典的状態を生成したいとします。 ここで、時刻\(t\) (整数)における状態を\(\Sigma_t\)と書きます。 狙っている分布を\(W(\Sigma)\)とすると、 ボルツマン分布の場合には、

次に、天下り的ですが、以下の条件を満たす、\(\Sigma\)から\(\Sigma'\)への遷移確率を\(T(\Sigma|\Sigma')\)を考えます。

  • \(\sum_{\Sigma'} T(\Sigma|\Sigma') W(\sigma') = \sum_{\Sigma'} T(\Sigma'|\Sigma) W(\sigma)\)（釣り合い条件)

  • 任意の2つの状態間\(\Sigma\), \(\Sigma'\)の間を有限回数の状態更新で結ぶことができる（エルゴード性）

釣り合い条件は、ある時刻\(t\)における分布関数\(P_t(\Sigma)\)が狙っている\(W(\Sigma)\)と一致しているとき、 状態\(\Sigma\)に流れ込む確率の流れ (左辺）と、流れ出す確率の流れ (右辺)が釣り合うべしという意味です。 一方、エルゴード性が満たされていない場合、適当な初期状態から出発したとき、どんなに更新を重ねても到達できない状態が存在するので、 期待される分布が得られない可能性があります。 適当な初期分布関数\(P_0(\Sigma)\)に従って初期状態\(\Sigma_0\)を生成し、 遷移確率に従って更新を重ねていると、\(\lim_{t\rightarrow +\infty} P_t(\Sigma) = W(\Sigma)\)となることが示されます。

実際の計算では、釣り合い条件よりも強い詳細釣り合い条件\(T(\Sigma|\Sigma') W(\sigma') = T(\Sigma'|\Sigma) W(\sigma)~(\forall \Sigma, \Sigma')\)を課すことが多いです。 詳細釣り合い条件は、任意の2つの状態間で確率の流れが釣り合うという条件を表しています。

イジングスピン模型の場合

欲しい分布、つまり、ボルツマン分布は、

\[ W(\Sigma) = Z^{-1} \exp(-\beta E(\Sigma)) \]

と書けます。ここで、\(\beta\)は逆温度、\(E(\Sigma)\)は状態\(\Sigma\)のエネルギー、\(Z\)は配位関数です。 この場合、以下で示す「シングルスピンアップデート」が釣り合い条件と、エルゴード性を満たすことを示します。 シングルスピンフリップでは、\(N_\mathrm{s}\)個のスピンから毎回ランダムに一個スピンを選び、 そのスピンに対して以下の更新を行います。 今、\(i\)番目のスピンが選ばれていて、現在のスピン状態を\(\Sigma_t = (S_1^t, \cdots, S_{N_\mathrm{s}}^t)\)としましょう。 遷移先の候補として、\(i\)番目のスピン以外を固定したすべての状態を考えます。 イジングスピンは離散自由度しか持たないので、 \(\Sigma_\uparrow = (S_1^t, \cdots, \uparrow, \cdots, S_{N_\mathrm{s}}^t)\)、 \(\Sigma_\downarrow = (S_1^t, \cdots, \downarrow, \cdots, S_{N_\mathrm{s}}^t)\)の2通りしかありません（片方は現在の状態と一致する)。 ここで、次の状態\(\Sigma^{t+1}\)を、確率\(P_\uparrow = 1/(1+\exp(\beta \Delta E))\)で\(\Sigma_\uparrow\), 確率\(P_\downarrow = 1-P_\uparrow\)で\(\Sigma_\downarrow\)を選びます。

シングルスピンアップデート法は、以下のように解釈することも出来ます。 現在のスピン状態\(\Sigma_t\)から、注目しているスピン以外を固定した「全てのスピン状態」を遷移先の候補として考えます。 その状態の集合（＝状態空間の部分集合）の中で、ボルツマン分布に従って、次のスピン状態を選びます。

演習問題

\(P_\uparrow = E(\Sigma_\uparrow)/(E(\Sigma_\uparrow) + E(\Sigma_\downarrow))\)を使って、上の状態更新が詳細釣り合い条件を満たしていることを証明しましょう。なお、\(i\)番目のサイトを選ぶ確率も考慮に入れましょう。

実際の計算では、毎回スピンをランダムに選ばず、\(N_\mathrm{s}\)個のスピンを順番に探査することが一般的です。 この場合、厳密な証明を書きませんが、詳細釣り合い条件は破れますが、釣り合い条件は満たされています (証明は、今見ているスピン番号を状態に含めた拡張された状態空間を考えます）。

疑似乱数

疑似乱数とは、決定的論的に生成される数列ではあるが、実質「乱数」と見なせるようなものです。 長年の研究から様々な手法が提案されていますが、Juliaにはメルセンヌツイスタと呼ばれるものが標準で実装されています。乱数の質（乱雑さ）、計算速度のバランスに優れており、この方法を使っていれば問題ないと思います。 より最新のアルゴリズム (xorshiftなど)が提案されていますが、物性物理の用途において有意な違いはでないと個人的には考えています。

下では、先ず1行目で疑似乱数生成器を、乱数の種 (seed=4649)を使って生成しています。 疑似乱数では同じseedを使うことで、デバッグ時の再現性を担保することが出来ます。 異なるseedをつかって作った疑似乱数列は、互いに無相関だと見なしてかまいません。

2行目では、1行目で作った乱数生成器rngを使って、{-1, 1}から同じ確率で選んだ10個の数列を配列として返します。

using Random
rng = MersenneTwister(4649)
rand(rng, Int8[-1, 1], 10)

10-element Array{Int8,1}:
  1
  1
 -1
 -1
 -1
  1
  1
  1
 -1
  1
Preparation 1D Ising model