ねぇ。Binder Parameterって知ってる？

はじめに

　突然ですが、イジングモデルのシミュレーションのデバッグってしたことありますか？このデバッグ方法として一般にBinder Parameterが用いられます（？）「臨界温度を見ておけばいいんだろ？」と思われたかもしれません。たしかに、無限系（無限個の粒子を想定）のイジングモデルでは、臨界温度パラメータ（特徴的な温度・沸点とか融点みたいなやつ）が解析解として一般に知られています。一方、シミュレーションは有限系であり、臨界温度は知られていません。十分スケールを大きくすればいいという話もありますが、それって正しいですか？かっこ悪くないですか？

対象読者

　解析解が得られていない問題に対して、シミュレーションを用いてそれを求めるにはどうしたら良いでしょうか？十分スケールを大きくしますか？それとも知恵を振り絞りますか？この記事は、後者の人に向けた内容です。チェッカーボードアルゴリズムとか効率的なシミュレーションの話ではないです。

モチベーション

　もしも系のスケールに依存しないパラメータが分かればと思いませんでしたか？無限系の解析解と有限系のシミュレーションでスケール非依存のパラメータが一致するならば、正しいシミュレーションができたと安心できます。そんなパラメータが（臨界温度における）Binder Parameterです。私が探した範囲では、これのスケール依存性の計算が間違っているようにみえる証明・論文が多く、なかなか正しいものに巡り会えなかったです。（引用先が無関係だったり、引用たらい回しとか、、、正しそうに見えても、機械系出身には正しいか見分けがつかない鬼ハード証明だったり、、、）唯一納得できたのが、ドイツ語の資料でした。後続の人が困らないように思い出を書き残しておきます。

ドイツ語の資料

"Statistische Mechanik ungeordneter Systeme, SoSe 2005" でググると、いい感じのものがいっぱい出てくるんで参考にしてみてください。本記事は、以下を参考にしています。

考え方

　無限系における臨界温度をシミュレーション結果から間接的に求める。「スケールサイズに依らず」「臨界温度直上」で一定値を取る関数 B(L, T) を用いる。複数のスケールでのシミュレーション結果 B(L, T) に対して最小二乗法を用いて交点である臨界温度を求める。（本件では最小二乗法については触れない）解析解と最小二乗法によって得られた臨界温度を比較する。比較結果が妥当であれば、デバッグよし！！とする。
この2つ大事ね。「スケールサイズに依らず」「臨界温度直上」。わたしは頭が悪いのでここでハマった。

計算

Binder Parameterと呼ばれる関数値は、(1)の式で表される。この値が臨界温度直上で一定値を取ることを証明する。

\begin{eqnarray}
B(L,T) = 1- \frac{\langle M^4 \rangle}{3 \langle M^2 \rangle ^2} 
\end{eqnarray}
\tag{1}

ここで、$M$は磁化の和であり、全スピンの集合$D$を用いて表される変数

\begin{eqnarray}
M = \sum_{i,j\in{D}}^{N}{\phi _{i,j}}
\end{eqnarray}
\tag{2}

でありブラケット$\langle \hspace{10pt} \rangle$は、シミュレーション時間$\tau_\infty$を用いて

\begin{eqnarray}
\langle A^\omega(T) \rangle = \frac{1}{\tau_\infty} \sum_{\tau=0}^{\tau_\infty} A_\tau^\omega(T)
\end{eqnarray}
\tag{3}

と表される計算である。次に、このBinder Parameterが臨界温度$T_c$をでシミュレーション系のスケールサイズに依存しないことを示すため、いくつかの変数を定義する。$T_c$を用いて温度に関する無次元量$t$を以下で定義する。

\begin{eqnarray}
t = \left| \frac{T-T_c}{T_c} \right|
\end{eqnarray}
\tag{4}

また、臨界温度臨界温度$T=T_c$付近で以下の関係式

\begin{eqnarray}
\xi_{L,T} \sim t^{-\nu} \hspace{15pt} ( 0 < \nu)
\end{eqnarray}
\tag{5}

\begin{eqnarray}
M \sim t^\beta \hspace{15pt} (0 < \beta)
\end{eqnarray}
\tag{6}

が成り立つとき、$\kappa=\frac{\beta}{\nu}$を用いて

\begin{eqnarray}
M \sim \xi_{L,T} ^{-\kappa}
\end{eqnarray}
\tag{7}

となる。ここで、$\xi_{L,T}$は、シミュレーション系の代表長さ$L$と温度$T$を引数としてスピン変数の相関長をである。この式(7)に磁化の特徴的な値$m$を用いて

\begin{equation}
M = m \xi_{L,T} ^{-\kappa}
\end{equation}
\tag{8}

とする。また、シミュレーション系において、$t \sim 0$で$\xi_{L,T} \rightarrow L$が成り立つとみなすと

\begin{equation}
M = m L ^{-\kappa}
\end{equation}
\tag{9}

となる。これを用いるとシミュレーション系における$M$の実現確率$P_{L,T}(M)$は

\begin{equation}
P_{L,T}(M) = C_L P_{L,T}(mL^{-\kappa})
\end{equation}
\tag{10}

となる。ここで$C_L$は規格化定数であり

\begin{eqnarray}
\int_{- \infty}^{\infty} P_{L,T}(m)dm = 1
\end{eqnarray}
\tag{11}

から、$C_L=L^{-\kappa}$となる。この実現確率$P_{L,T}$の添字は、$P_{L,T}$の分布が$L,T$に依存することを表している。この$T$を無限系における相関長$\xi_{\infty,T}$で置き換えることを考える。無限系の相関長がシミュレーション系の式(5)と同様に$T=T_c$付近で

\begin{eqnarray}
\xi_{\infty,T} \sim t^{-\nu}
\end{eqnarray}
\tag{12}

が成り立つとすれば、$T$を$\xi_{\infty,T}$の関数とみなせるので

\begin{equation}
P_{L,T}(mL^{-\kappa}) = P_{L,\xi_{\infty, T}}(mL^{-\kappa}) 
\end{equation}
\tag{13}

となり、式(10)を用いると

\begin{eqnarray}
P_{L,T}(m) = C_L P_{L,\xi_{\infty, T}}(mL^{-\kappa}) 
\end{eqnarray}
\tag{14}

となる。この実現確率を用いて$M^\omega$の期待値を求める。まず、期待値$M^\omega$を変数$m$に対して連続化を行い、実現確率$P_{L,T}$を用いた積分で表すと

\begin{eqnarray}
\langle M^\omega \rangle =\int_{- \infty}^{\infty}m^\omega P_{L,T}(m)dm
\end{eqnarray}
\tag{15}

これに実現確率を式(15)に用いると

\begin{eqnarray}
\langle M^\omega \rangle = C_L \int_{- \infty}^{\infty}m^\omega P_{L,\xi_{\infty, T}}(mL^{-\kappa})dm
\end{eqnarray}
\tag{16}

となる。ここで、$\frac{L}{\xi_{\infty,T}}$のみによって$P_{\xi_{\infty,T}}(mL^{-\kappa})$の分布形状が決定するという仮定を用いれば

\begin{equation}
P_{L,\xi_{\infty, T}}(mL^{-\kappa}) = P \left( mL^{-\kappa}, \frac{L}{\xi_{\infty, T}} \right)
\end{equation}
\tag{17}

となるから、式(16)は

\begin{eqnarray}
\langle M^\omega \rangle = C_L \int_{- \infty}^{\infty}m^\omega P \left( mL^{-\kappa}, \frac{L}{\xi_{\infty, T}} \right)　dm
\end{eqnarray}
\tag{18}

となる。さらに、臨界温度直上で、$\xi_{\infty, T} \rightarrow \infty$であることを利用すれば

\begin{eqnarray}
\langle M^\omega \rangle = C_L \int_{- \infty}^{\infty}m^\omega P(mL^{-\kappa}, 0)　dm
\end{eqnarray}
\tag{19}

となり、分布の形状がスケールサイズに依存しない形式になった。$x=mL^{-\kappa}$とおいて置換積分をすると

\begin{eqnarray}
\langle M^\omega \rangle = C_L \int_{- \infty}^{\infty}m^\omega P(mL^{-\kappa}, 0)dm \nonumber \\
= C_L \int_{- \infty}^{\infty} (L^{\kappa} x)^\omega P(x, 0) L^\kappa dx \nonumber \\
=\int_{- \infty}^{\infty} (L^{\kappa} x)^\omega P(x, 0) dx \nonumber \hspace{27pt} \\
= L^{\kappa\omega} \int_{- \infty}^{\infty} x^\omega P(x, 0) dx \nonumber \hspace{27pt} \\
= L^{\kappa\omega} F(\omega) \hspace{80pt}
\end{eqnarray}
\tag{20}

となる。この式(20)を用いて、式(1)の臨界温度直上の$T=T_c$におけるBinder Parameterを求めると

\begin{eqnarray}
B(L,T_c) = \left . 1- \frac{\langle M^4 \rangle}{3 \langle M^2 \rangle ^2} \right |_{T=T_C}  \nonumber \\
= 1 - \frac{L^{4\kappa} F(4)}{\left( L^{2\kappa} F(2) \right)^2} \nonumber \hspace{14pt} \\
= 1-\frac{F(4)}{F^2(2)} \hspace{37pt}
\end{eqnarray}
\tag{21}

となる。以上より、Binder Parameterは$T=T_c$において、スケールサイズに依存しない定数値を取ることがわかる。

まとめ

個人的に分かりやすく納得のできる証明だったので、記事にしてみました。もっとシンプルな方法とかあるんでしょうか。。。
数値計算は、チェッカーボードアルゴリズムでやりました。使い慣れてるpython+tensorflowからGPUを活用しました。

おまけ

私のクソ論文からの引用画像（シミュレーション結果）

先行研究[21]は以下
Kun Yang, Yi-Fan Chen, Georgios Roumpos, Chris Colby, John Anderson.High Performance Monte Carlo Simulation of Ising Model on TPU Clusters. In Proc. 2019 ACM/IEEE Int. Conf. for High Performance Computing, Networking, Storage, and Analysis (SC’19). IEEE, 2019.
ほんとにそうか？と思ったんで詳しそうな人にtwitterで聞いてみた。

この時は外場がわずかにある条件で計算していたのがよくなかったみたいで、外場0でやったらある程度システムサイズ大きいとこでは割と綺麗に交わった記憶があります

— まそらったー (@masorata) November 30, 2021

境界条件が結構効くかもｗｗ 確かに、固定端にするのは熱浴に晒しているのと同じようなもんなんだしなぁ。 固定端の効果がどれくらいか気になってきたけど、筆を置きます。

雰囲気ソースコード

import numpy as np
import tensorflow as tf


class boundary_fix:
    def __init__(self, width, height):
        self.width = width
        self.height = height
        ones_map = np.ones((width, height))
        self.boundary_map_0 = self.fixed_boundrizer(ones_map)
        self.boundary_map_1 = np.ones((1, width, height, 1)) - self.boundary_map_0
        # self.boundary_map_1 = np.reshape(ones_map, (1, width, height, 1)) - self.boundary_map_0

    def boundry_tensor_maker(self, free_tensor, boundary_value=0.0):
        free_tensor = free_tensor * self.boundary_map_0
        free_tensor = free_tensor + boundary_value * self.boundary_map_1
        return free_tensor

    def fixed_boundrizer(self, free_map):
        free_map[0] = 0
        free_map[self.width - 1] = 0
        free_map[..., 0] = 0
        free_map[..., self.height - 1] = 0
        # free_map = np.reshape(free_map, (self.width, self.height))
        free_tensor = np.reshape(free_map, (1, self.width, self.height, 1))
        return free_tensor


def spins_initilizer(bound_fixer, map_shape=(1, 10, 10, 1)):
    using_initilizer = tf.initializers.random_uniform(minval=0.9, maxval=0.99)
    # using_initilizer = tf.initializers.random_uniform(minval=0.2, maxval=0.4)
    initiale_spins = tf.Variable(using_initilizer(map_shape), trainable=False, name="initiale_spiN")
    initiale_spins = bound_fixer.boundry_tensor_maker(initiale_spins, 1.0)
    # initiale_spins = tf.reshape(initiale_spins, (1, map_shape[0], map_shape[1], 1))
    initiale_spins = value01_2_sign(initiale_spins)
    return initiale_spins


def value01_2_sign(value01):
    round01 = tf.math.round(value01 + 1e-5)
    sign_plus_minus = 2 * round01 - 1
    return sign_plus_minus


def checker_flip_builder(flipper, checker_filter):
    tmp_flipper = tf.multiply(flipper, checker_filter)
    checker_one = tf.ones(shape=checker_filter.shape, dtype=tf.float64) - checker_filter
    checker_flipper = tf.math.add(tmp_flipper, checker_one)
    return checker_flipper


class spin_updater:
    def __init__(self, spins, energy_filter, kbt):
        self.spins = spins
        self.energy_filter = energy_filter
        self.kbt = kbt

    def get_new_spins(self, checker_filter):
        self.checker_filter = checker_filter
        energy_map = self.calc_energy()
        new_spins = self.spins_flipper(energy_map)
        return new_spins
        # return tf.nn.conv2d(self.spins, self.energy_filter, strides=(1, 1, 1, 1), padding="SAME")

    def calc_energy(self):
        interaction_map = tf.nn.conv2d(self.spins, self.energy_filter, strides=(1, 1, 1, 1), padding="SAME")
        energy_map = -tf.multiply(self.spins, interaction_map)
        return energy_map

    def spins_flipper(self, energy_map):
        energy_diff = -2.0 * energy_map
        accept_rate = tf.math.divide(tf.math.exp(-energy_diff / self.kbt),
                                     tf.random.uniform(shape=tf.shape(energy_map), minval=1e-8, maxval=1.0, dtype=tf.float64))
        tmp_flipper = tf.sigmoid(accept_rate - 1.0)
        flipper = -value01_2_sign(tmp_flipper)
        checker_flipper = checker_flip_builder(flipper, self.checker_filter)
        new_spins = tf.multiply(self.spins, checker_flipper)
        new_spins = tf.math.round(new_spins)
        return new_spins