【Python】モンテカルロ法の応用｜二次元イジング模型

本記事では、マルコフ連鎖モンテカルロ法の応用の一つであるイジング模型のシミュレーションを行っていきます。

イジング模型は統計力学という学問でよく出てくるモデルですが、近年量子計算や最適化計算でも注目されてきています。

本記事では統計力学をよく知らない方でも理解できるようにイジング模型のシミュレーションを紹介していきます。

＊未完成のため随時更新します（間違いがあったらコメントよろしくお願いします）

イジング模型の基礎知識
- エネルギーベースモデル
- イジング模型の確率分布
二次元イジングモデルのシミュレーション
二次元イジングモデルの実装（Python）
参考資料
- 参考文献
- 参考URL
まとめ

イジング模型の基礎知識

まずは、イジング模型（強磁性イジング模型）の基礎知識を簡単に説明していきます。

エネルギーベースモデル

機械学習の文脈でもよく出てくるエネルギーベースモデルは、エネルギー関数$E(\mathbf{x} ; \theta)$を用いて以下のように確率分布を定義するモデルです。

$$P(\mathbf{X} = \mathbf{x} \mid \theta) = \frac{e^{- \beta E(\mathbf{x} \mid \theta)}}{Z(\theta, \beta)}$$

ここで、規格化定数$Z(\theta)$は分配関数呼ばれ以下のように定義されます。

$$Z(\theta, \beta) = \sum_{\mathbf{x}} e^{- \beta E(\mathbf{x} \mid \theta)}$$

$\sum_{\mathbf{x}} \cdots $は$\mathbf{x}$の可能な状態全てに関する和を表しています（確率変数が連続値を取る場合は積分に変換してください）

$\beta$は一般に正の値を取り、逆温度パラメータと呼ばれたりします。

統計物理学を勉強したことがある方は、教科書通り$\beta \equiv 1/k_{B}T$で定義される量だと思ってください。

イジング模型の確率分布

（強磁性）イジングモデルは、確率変数が離散値$x_{i} \in \{-1, 1\}$の値を取り格子状にエネルギー関数が以下のように定義されるエネルギーベースモデルの一つです。

$$E(\mathbf{x} ; \mathbf{J}) = – J \sum_{\langle i, j \rangle} x_{i} x_{j} – h \sum_{i} x_{i}$$

各$i$は格子点を整数でラベルしたものを表し、各格子点上に確率変数$x_{i}$を定義します。

二次元イジング模型では、確率変数のことを物理の応用を視野に入れてスピンと呼ぶことも多いので、本記事では以降スピンと呼ぶことにします。

また、$\langle i, j \rangle$は隣接$i, j$の組みを表し、総和は隣接する変数同士の和を取ることを意味します。

格子の形状によって様々な振る舞いが見られます。

二次元イジングモデルのシミュレーション

今回具体的にシミュレーションするモデルは有名な正方格子上のイジングモデルです。

使用する手法は、詳細つりあい条件を満たすモンテカルロ法である『メトロポリス法』と『ギブスサンプリング』を説明します。

二次元正方格子の2次元イジングモデル

二次元格子点上に確率変数を定義するため格子点を$(i, j)$と表し、格子点上に定義される確率変数を$x_{ij}$とあらわあすことにします。

$n \times n$の格子点なら$i, j$はそれぞれ$i=j=\{1, \ldots, n\}$の整数値となります。

二次元正方格子イジングモデルは、以下のエネルギー関数で定義されるエネルギーベースモデルです。

\begin{align}E(\mathbf{x};J) &= -J \sum_{\langle i,j \rangle} x_{i, j} x_{i,j} \\ &= -J \sum_{ij} x_{i,j} (x_{i-1,j} + x_{i+1,j} + x_{i,j-1} + x_{i,j+1}) \end{align}

具体的な確率分布はエネルギー関数を代入して以下のように表せます。

$$P(\mathbf{x} \mid \beta, J) = \frac{1}{Z(\beta, \theta)} \exp \left(\beta J \sum_{\langle i,j \rangle} x_{i, j} x_{i,j} \right)$$

分配関数を具体的に表すと以下のようになります。

$$Z(\beta, J) = \sum_{\mathbf{x}} \exp \left(\beta J \sum_{\langle i,j \rangle} x_{i, j} x_{i,j} \right)$$

ここで、$\sum_{\mathbf{x}}$は$\mathbf{x}$のとりうる全ての状態に関する和なので、一般に計算するは困難です。

しかし、後述のマルコフ連鎖モンテカルロ法を使用すれば、この問題を回避することができます。

二次元イジングモデルの目的

イジングモデルは統計物理学という分野で磁性体の性質を明らかにするために導入されたモデルです。

二次元イジングモデルの場合、厳密に統計量を評価することができ、その統計量に特異性が現れる相転移という現象を見ることができます。

より詳しく知りたい方は下記で紹介している統計力学の参考書で勉強してみてください。

メトロポリス法

メトロポリス法の受理確率は、エネルギーベースモデルの場合以下のように計算できます。

\begin{align}A(\mathbf{x}, \mathbf{x}^{\prime}) &= \min\left(1, \frac{P(\mathbf{x}^{\prime})}{P(\mathbf{x})}\right) \\ &= \min \left(1, \exp \left(\beta (E(\mathbf{x};J) ~- E(\mathbf{x}^{\prime};J)) \right) \right) \end{align}

この受理確率からもわかるようにマルコフ連鎖を構成する際に分配関数$Z(\beta, J)$の計算を行う必要がなくなります！

また、ランダムに格子点$(i,j)$を選んで$x_{i,j}^{\prime} \to ~- x_{ij}$にフリップするような提案確率を選ぶと受理確率は以下のように表せます（実際は、ランダムではなくあらかじめ決まった順番を取るように提案しても良い）

$$A(\mathbf{x}, \mathbf{x}^{\prime}) = \min \left(1, \exp\left(-2J\beta x_{i,j}(x_{i-1,j}+x_{i+1,j}+x_{i, j-1}+x_{i,j+1}) \right) \right)$$

遷移確率を求めることができたのであとはマルコフ連鎖を実際に構成するだけです。

具体的には、以下のようにアルゴリズムを構成します。

メトロポリス法

Input : 任意の初期状態$\mathbf{x}^{0}$

ランダムに$(i, j)$を選び、$x_{ij}^{\prime} \to ~-x_{ij}$のようにフリップする
受理確率$A(\mathbf{x}, \mathbf{x}^{\prime})$を計算する
$0 \le r <1$の一様乱数$r$を発生し、$r < A(\mathbf{x}, \mathbf{x}^{\prime})$なら$x_{i}^{\prime}$を受理し、そうでなければ棄却する
①〜③を繰り返す

ギブスサンプリング

正方格子二次元イジングの場合は、各確率変数$x_{i}$の条件付き確率が以下のように簡単に求めることができるのでギブスサンプリングを利用することができます。

$$P(x_{ij} \mid \mathbf{x} \backslash{x_{ij}}) = \frac{\exp\left(J\beta x_{i,j}(x_{i-1,j}+x_{i+1,j}+x_{i, j-1}+x_{i,j+1})\right)}{2 \cosh \left(J\beta x_{i,j}(x_{i-1,j}+x_{i+1,j}+x_{i, j-1}+x_{i,j+1}) \right)}$$

あとは、この条件付き確率を利用してマルコフ連鎖を構築するだけです。

アルゴリズム的には以下のように行います。

ギブスサンプリング

ランダムに$i$を選ぶ
$P(x_{ij}=1 \mid \mathbf{x} \backslash{x_{ij}})$を計算する
$0 \le r <1$の一様乱数$r$を発生し、$r < P(x_{ij}=1 \mid \mathbf{x} \backslash{x_{ij}})$ならば$x_{ij} \leftarrow 1$とし、そうでなければ$x_{ij} \leftarrow -1$とする
①〜③を繰り返す

二次元イジングモデルで評価する量

二次元イジングモデルでは、以下で定義される磁化という量を評価します。

$$m = \frac{1}{n^{2}} \sum_{i, j} x_{ij}$$

この量がある逆温度の値で0から有限値になる性質があります。

また、磁化が０から有限値になる際に以下で定義される比熱が$n \to \infty$の極限で発散することが知られています。

$$C = \beta^{2}(\langle E^{2} \rangle – \langle E \rangle^{2})$$

比熱の発散は$n \to \infty$でのみ生じますが、有限の$n$でもその片鱗を見ることができます。

このような量が発散することを相転移といいます。

二次元イジングモデルの実装（Python）

ここからは、二次元イジングモデルの実装方法を紹介します。

具体的には、以下の手順で実装を行います。

最近接格子のスピンの和を計算する関数を定義
エネルギーを計算する関数を定義
メトロポリス法・ギブスサンプリングを定義
マルコフ連鎖モンテカルロ法を実装
各逆温度ごとにマルコフ連鎖モンテカルロ法を実装

今回は$x$方向の格子の数を$N_{x}$、$y$方向の格子の数を$N_{y}$としました。

最近接格子のスピンの和を計算する関数を定義

二次元イジングモデルにマルコフ連鎖モンテカルロ法を適用する際、最近接格子のスピンの和を計算することがよくあります。

そのため、最近接格子の和を計算する関数を定義しておくと便利です。

具体的には、以下のコードになります。

def neighor_spin_sum(s, x, y):
    """近接スピンの和を計算

    Parameters
    --------------------
    s : スピン配位
    x : 格子点のx座標
    y : 格子点のy座標

    Returns
    --------------------
    neighor_spin_sum : 近接スピンの和
    """
    x_right=x+1
    x_left=x-1
    y_up=y+1
    y_down=y-1

    # 周期境界条件
    if x_right>=Nx:
        x_right-=Nx
    if x_left<0:
        x_left+=Nx
    if y_up>=Ny:
        y_up-=Ny
    if y_down<0:
        y_down+=Ny

    neighor_spin_sum=s[x_right][y]+s[x_left][y]+s[x][y_up]+s[x][y_down]
    return neighor_spin_sum

格子の形が異なる場合は、この部分を変更すればOKです。

また、境界部分は周期境界条件を課しました。

エネルギーを計算する関数を定義

エネルギー期待値を評価するためにエネルギーを計算するためのコートを定義しておきます。

def calc_energy(s, h=0.01):
    """エネルギーを計算

    Parameters
    --------------------
    s : スピン配位
    h : 外部磁場(バイアス)

    Returns
    --------------------
    energy : エネルギー
    """
    energy = 0
    for x in range(Nx):
        for y in range(Ny):
            # dobule count
            energy += - neighor_spin_sum(s, x, y)/2
    energy += h*np.sum(s)
    return energy

先ほど定義した『neighor_spin_sum()』を使用すれば簡単に計算することができます。

しかし、単純に和を取るだけだと二重に項を計算してしまうので２で割るこを忘れないようにしましょう。

メトロポリス法・ギブスサンプリングを定義

一回のメトロポリス法・ギブスサンプリングによるアップデートを定義します。

メトロポリス法のアップデートは以下のように定義することができます。

def metropolis(s, beta=1, h=0.001):
    """メトロポリス法

    Parameters
    --------------------
    s : スピン配位
    beta : 逆温度パラメータ
    h : 外部磁場(バイアス)

    Returns
    --------------------
    s : 更新後のスピン配位
    """
    xs=list(range(Nx))
    random.shuffle(xs)
    ys=list(range(Ny))
    random.shuffle(ys)
    for x in xs:
        for y in ys:
            k=neighor_spin_sum(s, x, y)+h
            trans_prob = np.exp(-2*beta*s[x][y]*k)
            if np.random.random()<=trans_prob:
                s[x][y] = -s[x][y]
    return s

ギブスサンプリングの場合は、以下のように定義できます。

def gibbs_sampling(s, beta=1.0, h=0.001):
    """ギブスサンプリング

    Parameters
    --------------------
    s : スピン配位
    beta : 逆温度パラメータ
    h : 外部磁場(バイアス)

    Returns
    --------------------
    s : 更新後のスピン配位
    """
    xs=list(range(Nx))
    random.shuffle(xs)
    ys=list(range(Ny))
    random.shuffle(ys)
    for x in xs:
        for y in ys:
            k=neighor_spin_sum(s, x, y)-h
            trans_prob = np.exp(beta*k) / (np.exp(beta*k)+np.exp(-beta*k))
            if np.random.random()<=trans_prob:
                s[x][y]=1
            else:
                s[x][y]=-1
    return s

メトロポリス法もギブスサンプリングも今回はランダムにフリップするスピン配位を選び採択するか棄却するかを決定しています。

マルコフ連鎖モンテカルロ法を実装

ここまで定義してきた関数を使用してマルコフ連鎖を実装します。

具体的には以下のように実装できます。

def mcmc(s, steps=1000, beta=1.0, h=0.001, method='metropolis', interval=10, burn_in=10**2, animation=True, ani_interval=10):
    """マルコフ連鎖モンテカルロ法を実行
    
    Parameters
    -------------------
    s : 初期スピン配位
    steps : モンテカルロステップ
    beta : 逆温度パラメータ
    h : 外部磁場(バイアス)
    method : モンテカルロ法の手法, メトロポリス法 : metropolis, ギブスサンプリング : gibbs
    interval : サンプルを取得するインターバル
    burn_in : バーンインタイム

    Returns
    -------------------
    ms : 各ステップの磁化のリスト
    energies : 各ステップのエネルギーのリスト
    square_energies : 各ステップのエネルギーの二乗のリスト
    """
    ms=[]
    energies=[]
    square_energies=[]
    # Metropolis 
    if method=='metropolis':
        for step in range(steps):
            s = metropolis(s, beta=beta, h=h)
            # 自己相関が切れる部分かつburn inを捨てる
            if (step%interval==0)&(step>=burn_in):
                # 磁化を計算
                m = np.sum(s) / (Nx*Ny)
                # エネルギーを計算
                energy = calc_energy(s, h=h)
                # エネルギーの二乗を計算
                square_energy = energy**2 
                ms.append(m)
                energies.append(energy)
                square_energies.append(square_energy)

    # Gibbs
    elif method=='gibbs':
        for step in range(steps):
            s = gibbs_sampling(s, beta=beta, h=h)
            if (step%interval==0)&(step>=burn_in):
                # 磁化を計算
                m = np.sum(s) / (Nx*Ny)
                # エネルギーを計算
                energy = calc_energy(s, h=h) 
                # エネルギーの二乗を計算
                square_energy = energy**2
                ms.append(m)
                energies.append(energy)
                square_energies.append(square_energy)
    return ms, energies, square_energies

『method』の引数を変えることでメトロポリス法, ギブスサンプリングどちらを利用するかを選択できるようにしました。

マルコフ連鎖モンテカルロ法を実行

以下のコードを実行することでマルコフ連鎖モンテカルロ法を実行することができます。

# 格子の大きさ
Nx = 32
Ny = 32
steps=10**4
# 初期配位
s = np.ones((Nx, Ny)).tolist()
ms, energies, square_energies=mcmc(s, steps=steps, beta=0.1, h=0.0)

実際に得られた磁化の結果を表示してみましょう。

# 磁化
time = np.arange(1, len(ms)+1)
fig, ax=plt.subplots()
ax.plot(time, ms)
ax.set_ylim([-1.1, 1.1])
ax.set_xlabel('MCS', fontsize=20)
ax.set_ylabel('Magnetization', fontsize=20)
plt.tick_params(axis='y', width=1, length=4, pad=5, color='k', direction='in', labelsize=15, labelcolor='k')
plt.tick_params(axis='x', width=1, length=4, pad=5, color='k', direction='in', labelsize=15, labelcolor='k')
plt.show()

下記のコードでアニメーションをみることができるので逆温度の値を変えて遊んでみてください。

Nx = 256
Ny = 256
steps=100
# 初期配位
s = np.random.randint(0, 2, (Nx, Ny)).tolist()
fig, ax = plt.subplots()
ims = []
im = ax.imshow(s, animated=True)
ims.append([im])
for step in range(steps):
    s = gibbs_sampling(s, beta=0.44, h=0)
    im = ax.imshow(s, animated=True)
    ims.append([im])

# ArtistAnimationにfigオブジェクトとimsを代入してアニメーションを作成
anim = animation.ArtistAnimation(fig, ims)
# Google Colaboratoryの場合必要
rc('animation', html='jshtml')
plt.close()
anim

各逆温度ごとにマルコフ連鎖モンテカルロ法を実行

最後に各逆温度ごとにマルコフ連鎖モンテカルロ法を実行するコードを紹介します。

def mcmc_beta(s, betas, total_steps=10**5, burn_in=10**2, interval=20, h=0.01):
    """各温度毎にMCMCを実行

    Returns
    --------------------
    s : 初期配位
    betas : MCMCを走らせる温度のリスト
    total_steps : モンテカルロステップ
    burn_in : バーンインタイム
    interval : サンプルを取得するインターバル
    h : 外部磁場

    Returns
    -------------------
    mags : 各温度の磁化
    mags_err : 各温度の各磁化の標準誤差
    energies : 各温度のエネルギー
    energies_err : 各温度の各エネルギーの標準誤差
    capacities : 各温度の比熱
    """
    b_mags=[]
    b_mags_err=[]
    b_energies=[]
    b_energies_err=[]
    b_capacities=[]
    for beta in betas:
        ms, energies, square_energies=mcmc(s, steps=total_steps, beta=beta, h=h, interval=interval, burn_in=burn_in)
        # 磁化
        b_mag=np.mean(ms)
        b_mag_err=np.std(ms)/np.sqrt(len(ms)-1)
        # エネルギー
        b_energy = np.mean(energies)
        b_energy_err = np.std(energies)/(np.sqrt(len(ms)-1))
        # 比熱
        b_square_energy = np.mean(square_energies)        
        b_capacity = (beta**2)*(b_square_energy - b_energy**2)

        b_mags.append(b_mag)
        b_mags_err.append(b_mag_err)
        b_energies.append(b_energy)
        b_energies_err.append(b_energy_err)
        b_capacities.append(b_capacity)
    return b_mags, b_mags_err, b_energies, b_energies_err, b_capacities

各温度ごとに磁化・エネルギー・比熱を計算しています。

具体例として$64 \times 64$の正方格子の場合を紹介します（多少時間がかかる場合は格子のサイズをもう少し小さくして実行してみてください）

Nx = 64
Ny = 64
steps=10**4
# 初期状態
s = np.ones((Nx, Ny)).tolist()
interval=20
# サンプルする逆温度の値 
betas=np.linspace(0.05, 1.0, 20).tolist()
# mcmcを実行
b_mags, b_mags_err, b_energies, b_energies_err, b_capacities = mcmc_beta(s, betas, total_steps=10**3)

得られた磁化は以下のようになります。

fig, ax=plt.subplots()
betas = np.array(betas)
b_mags = np.array(b_mags)
b_mags_err = np.array(b_mags_err)
ax.errorbar(betas, b_mags, yerr=b_mags_err, capsize=4)
ax.set_xlabel(r'$\beta$', fontsize=20)
ax.set_ylabel('Magnetization', fontsize=20)
plt.tick_params(axis='y', width=1, length=4, pad=5, color='k', direction='in', labelsize=15, labelcolor='k')
plt.tick_params(axis='x', width=1, length=4, pad=5, color='k', direction='in', labelsize=15, labelcolor='k')
plt.show()

磁化が0から有限の値に変化する様子が見られます。

次に比熱の値を図示します。

fig, ax = plt.subplots()
betas = np.array(betas)
b_capacities = np.array(b_capacities)
ax.plot(betas, b_capacities)
ax.set_xlabel(r'$\beta$', fontsize=20)
ax.set_ylabel('Capacity', fontsize=20)
plt.tick_params(axis='y', width=1, length=4, pad=5, color='k', direction='in', labelsize=15, labelcolor='k')
plt.tick_params(axis='x', width=1, length=4, pad=5, color='k', direction='in', labelsize=15, labelcolor='k')
plt.show()