初投稿になります！本記事では，著者自身が最近勉強している「イジングモデル」と「古典アニーリング」に関連する概念をまとめています．
概要
アニーリングとは何たるかをつかむため，数式を使って理論的背景を追っていきます．
1. スピンの導入
古典力学の世界では，ポテンシャルエネルギー$\,U\,$を受ける質量$\,m\,$の質点は運動方程式$\,m\ddot{\boldsymbol{r}}=-\nabla U(\boldsymbol{r})\,$に従って運動します．これは「エネルギーの低いほうに，物体は力を受けて動く」ことを意味していました．統計力学の文脈でも同様に，
系はエネルギー$\,\mathcal{H}\,$を最小化するように時間発展する
というのは，一つの大事なポイントだと思います．この概念を念頭に置いて，次にスピン変数を導入します．スピン変数は
$$
\sigma \in \{+1,-1\}\equiv\{\uparrow,\downarrow\}
$$の2値をとる変数として定義されます，$\,\sigma=+1\,(\uparrow)\,$ を上向きスピンと呼び，それと符号が逆の$\,\sigma=-1\,(\downarrow)\,$を下向きスピンと呼びます．このスピン変数が$\pm1$どちらの値をとるか，言い換えるとスピンがどちらの方向を向くかは，系のエネルギーがどのように設定されるかに依ります．
1.1. 具体的なエネルギー設定の例その1
例えば，$\,h\,$を定数(磁場)としてエネルギー $\ \mathcal{H}=-h\cdot \sigma\,$ の場合を考えます．このとき，スピンがとりうる状態とエネルギーの値の組み合わせは以下の2通りになります．
| $\,\sigma\,$ | $\mathcal{H}$ |
|---|---|
| $+1$$\,(\,\uparrow \,)$ | $\ -h\,$ |
| $-1$$\,(\,\downarrow \,)$ | $\ +h\,$ |
この系において，スピンが$\,\uparrow,\downarrow\,$のどちらの状態の方がエネルギーが最小になるかは磁場$\,h\,$の符号によって決まります．符号について場合分けを行うと
磁場 $ h>0 $ のとき$\,\sigma\,$は$\,\uparrow\,$を向きたがる．[最小エネルギー値：$-h\,$]
磁場 $h <0 $ のとき$\,\sigma\,$は$\,\downarrow\,$を向きたがる．[最小エネルギー値：$+h\,$]
の2パターンに絞られます．上の項目で「向きたがる」という言葉は，スピンがその方向を向くと，エネルギーが最低値をとることを意味しています．ちなみに，$\mathcal{H}$の値が系がとりうるエネルギー$\mathcal{H}$の最小値と一致するとき，この状態を基底状態と呼びます．
上記の設定で，磁場$h$は人が入力できるパラメータです．$h$の正の向きを座標軸上向きにとると，「この$\,\mathcal{H}\,$のもとでは磁場$\,h\,$とスピン$\,\sigma\,$が同じ方向を向く(と基底状態になる)」ことが主張されます．以下画像のようなイメージです．
1.2. 具体的なエネルギー設定の例その2
次にスピン変数が2つになった$\,\mathcal{H}=-J\sigma_1\cdot \sigma_2\,$の場合を考えていきます．$J$は2スピン間の相互作用(結合)を表すパラメータです．この系では，2体それぞれが$\uparrow,\downarrow$の２値をとるため，合計$\,2\times 2=4\,$パターンの状態が考えられます．
| $\,\sigma_1\,$ | $\,\sigma_2\,$ | $\mathcal{H}$ |
|---|---|---|
| $+1$$\,(\,\uparrow \,)$ | $+1$$\,(\,\uparrow \,)$ | $-J$ |
| $+1$$\,(\,\uparrow \,)$ | $-1$$\,(\,\downarrow \,)$ | $+J$ |
| $-1$$\,(\,\downarrow \,)$ | $+1$$\,(\,\uparrow \,)$ | $+J$ |
| $-1$$\,(\,\downarrow \,)$ | $-1$$\,(\,\downarrow \,)$ | $-J$ |
同様に場合分けを行うと,
結合 $ J>0 $ のとき$\,\sigma_1,\sigma_2\,$は同じ方向を向きたがる．[最小エネルギー値：$-J\,$]
結合 $ J<0 $ のとき$\,\sigma_1,\sigma_2\,$は逆の方向を向きたがる．[最小エネルギー値：$+J\,$]
スピン同士が同じ方向を向く$J>0\,$の状態を強磁性結合と呼び，反対に$J<0\,$の互いに逆向きを向く状態を反強磁性結合と呼びます．「この$\,\mathcal{H}\,$のもとでは相互作用$\,J\,$の値によってスピン同士の向きを同じ向き・逆向きにできる」ことになります．
以上の例はスピン変数が１つ，２つでしたが，もう少し一般的な複数のスピンを考えると$\,\mathcal{H}\,$は次式の形になります．この一般形のエネルギーの式はイジングハミルトニアンと呼ばれます．
\mathcal{H}=-\sum_{i\in V}h_i\sigma_i-\sum_{(i,j)\in E}J_{ij}\sigma_i\sigma_j
$i,j$はスピン変数を識別するインデックスです．また，式中の$V$は磁場がかかっているスピンの集合，$E$は相互作用の設定されたスピンの集合を意味します．例えば，以下図のようにスピンが配置されているとしましょう．
この場合，磁場は1番のスピンにのみ印加されているので$V=\{1\}$となり，相互作用はそれぞれ1-2，2-3の間に設定されているので，$E=\{(1,2),\ (2,3)\}$と表されるわけです．
ここで，具体的に$h_1=+2,\, J_{12}=+1,\, J_{23}=-3$と与えた状況を考えてみましょう．ハミルトニアンは以下のようになります．
\begin{align*}
\mathcal{H}&=-h_1\sigma_1-J_{12}\sigma_1\sigma_2-J_{23}\sigma_2\sigma_3\\[7pt]
&=-2\,\sigma_1-\sigma_1\sigma_2+3\,\sigma_2\sigma_3
\end{align*}
まず，1番のスピンは磁場と向きが揃うほうがエネルギー値が低くなるので上を向きます．次に，スピン(1,2) 間は$J_{12}>0$であり強磁性的なので1番のスピンと同じ方向($\,\uparrow\,$)になります．最後に，スピン(2,3) 間は反強磁性相互作用が設定されているので2番のスピンと逆向き($\,\downarrow\,$)を向くと基底状態になります．参考に，全状態の組み合わせパターンを表として示しておきます．
| $\,\sigma_1\,$ | $\,\sigma_2\,$ | $\,\sigma_3\,$ | $\mathcal{H}$ |
|---|---|---|---|
| $+1$$\,(\,\uparrow \,)$ | $+1$$\,(\,\uparrow \,)$ | $+1$$\,(\,\uparrow \,)$ | $0$ |
| $+1$$\,(\,\uparrow \,)$ | $+1$$\,(\,\uparrow \,)$ | $-1$$\,(\,\downarrow \,)$ | $-6$ |
| $+1$$\,(\,\uparrow \,)$ | $-1$$\,(\,\downarrow \,)$ | $+1$$\,(\,\uparrow \,)$ | $-4$ |
| $+1$$\,(\,\uparrow \,)$ | $-1$$\,(\,\downarrow \,)$ | $-1$$\,(\,\downarrow \,)$ | $2$ |
| $-1$$\,(\,\downarrow \,)$ | $+1$$\,(\,\uparrow \,)$ | $+1$$\,(\,\uparrow \,)$ | $6$ |
| $-1$$\,(\,\downarrow \,)$ | $+1$$\,(\,\uparrow \,)$ | $-1$$\,(\,\downarrow \,)$ | $0$ |
| $-1$$\,(\,\downarrow \,)$ | $-1$$\,(\,\downarrow \,)$ | $+1$$\,(\,\uparrow \,)$ | $-2$ |
| $-1$$\,(\,\downarrow \,)$ | $-1$$\,(\,\downarrow \,)$ | $-1$$\,(\,\downarrow \,)$ | $4$ |
基底状態は，最もエネルギーが低くなる，表の上から2番目の状態に対応します．
ここら辺になってくると表を書いて基底状態を求めるのはしんどくなってきます...
上の例ではスピンが３つなので全状態のパターンは８通りですが，一般に$\,N\,$個のスピンがとりうる状態の数は$\,2^N\,$通りなので指数増加的し，莫大な組み合わせ数になります．アニーリングは，この多数の組み合わせから
ハミルトニアン$\,\mathcal{H}\,$を最小化するようなスピンの向きの組を近似的に求める
手段です．
アニーリングを必要とする場面として，例えば物理の世界では，2次元格子点にスピンが配置された次の図のような模型がよく登場します．
このような模型をイジング模型と呼びます．スピン同士の結合がすべて強磁性結合であればスピンの向きはすべて同じ方向を向き，すべて反強磁性結合ならば隣り合うスピン同士が逆を向く配置を取ります．では，次の図の場合はどうでしょうか？
赤線は強磁性結合で，青線は反強磁性結合を表しています．このように相互作用の符号がランダムに配置されると，途端に基底状態を求めるのが困難になります．上図のような強磁性と反強磁性がバラバラに入り混じった構造をスピングラスと呼びます．アニーリングはこのスピングラス構造において，基底状態を(近似的に)求めることを得意とする手法です．
2. 統計力学における確率の記述
2.1. 定常状態での確率
統計力学で記述される系において，各状態に通し番号$\,n\,$をふることを考えます．$n$番目の状態におけるスピンの組を$\,z_n\,$で表します．例えば，上述した$\mathcal{H}=-J\sigma_1\cdot\sigma_2$でハミルトニアンが設定される2体スピン系では，以下のように順に番号を与えればokです．
| $z_n$ | $\,\sigma_1\,$ | $\,\sigma_2\,$ | $\mathcal{H}$ |
|---|---|---|---|
| $z_1$ | $+1$$\,(\,\uparrow \,)$ | $+1$$\,(\,\uparrow \,)$ | $-J$ |
| $z_2$ | $+1$$\,(\,\uparrow \,)$ | $-1$$\,(\,\downarrow \,)$ | $+J$ |
| $z_3$ | $-1$$\,(\,\downarrow \,)$ | $+1$$\,(\,\uparrow \,)$ | $+J$ |
| $z_4$ | $-1$$\,(\,\downarrow \,)$ | $-1$$\,(\,\downarrow \,)$ | $-J$ |
統計力学では，系の温度$\,T\,$で状態$\,z_n\,$をとる確率は以下のように記述されます．
p(z_n) = \frac{1}{Z}e^{-\beta \mathcal{H}(z_n)}\qquad Z=\sum_n e^{-\beta \mathcal{H}(z_n)}\quad [\,\beta=1/T\,]
ここで，$\mathcal{H}(z_n)$は状態$\,z_n\,$におけるエネルギー値であり，表の一番右の列に対応します．状態$z_n$をとる確率が$\,e^{-\beta \mathcal{H}(z_n)}\,$の因子で与えられる確率分布をボルツマン分布と呼びます．また，$Z$は確率の総和が１になるための規格化因子です
2.1.1. 具体的な確率分布描画
単一スピンで定常磁場がかかっているときの確率分布を考えてみましょう．ハミルトニアンは$\ \mathcal{H}=-h\cdot \sigma\,$で与えられ，全状態は以下の表のようになるのでした．
| $\ z_n $ | $\,\sigma\,$ | $\mathcal{H}$ |
|---|---|---|
| $\ z_1$ | $+1$$\,(\,\uparrow \,)$ | $\ -h\,$ |
| $\ z_2$ | $-1$$\,(\,\downarrow \,)$ | $\ +h\,$ |
規格化因子を計算すると
\begin{align}
&Z=\sum_ne^{-\beta \mathcal{H}(z_n)}=e^{\beta h}+e^{-\beta h}=2\cosh(\beta h)
\end{align}
となるので，$\ \sigma=\uparrow,\downarrow$ が得られる確率はそれぞれ
\begin{align}
&p(z_1)=\frac{e^{\beta h}}{2\cosh(\beta h)},\quad p(z_2)=\frac{e^{-\beta h}}{2\cosh(\beta h)}
\end{align}
と与えられます．磁場$\ h=1$としてそれぞれの状態をとりうる確率をプロットしてみます．
上向きの磁場$\,h$とスピン$\sigma$がそろう方がエネルギーが低くなり，取りうる確率が高くなるので，$\ \beta$が上がるにつれて$\ z_1$状態に遷移していく様子が確認できます．
2.2. 確率の時間発展
時刻$\ t\ $で状態$\,z_n\,$をとる確率を$\,p(z_n,t)\,$と書くとき，時間発展する確率微分方程式は
\begin{align}
p(z_n,t+\Delta t)=-\sum_{m(\neq n)}p(z_n,t)\,\phi_{n\to m}\,\Delta t+\sum_{m(\neq n)}p(z_m,t)\,\phi_{m\to n}\,\Delta t
\end{align}
と定義されます．$\ \phi_{n\to m}\ $は状態$\,z_n\to z_m\,$への単位時間あたりの遷移確率です．第1項は状態$\,z_n\,$からほかの状態に遷移する確率で，第2項はほかの状態から$\,z_n\,$へ遷移する確率を表しています．ここで，記号の簡略化のため
\begin{align}
\mathcal{L}_{nm} \equiv \phi_{m \to n},\quad \mathcal{L}_{nn}=-\sum_{m(\neq n)}\mathcal{L}_{mn}
\end{align}
を定義すると，元の確率微分方程式は
\begin{align}
\frac{p(z_n,t+\Delta t)}{\Delta t}=\mathcal{L}_{nn}\ p(z_n,t)+\sum_{m(\neq n)}\mathcal{L}_{nm}\, p(z_m,t)
\end{align}
と書き直すことができます．$\Delta t\to 0$の極限をとると，以下の「マスター方程式」が導出されます．
\begin{equation}
\frac{dp(z_n,t)}{dt}=\sum_{m} \mathcal{L}_{nm}\ p(z_m,t)
\end{equation}
時間に依存しない定常状態(つまり左辺=0となる状態)では，系はボルツマン分布
\begin{equation}
p_{\mathrm{eq}}(z_n)=\frac{1}{Z} e^{-\beta \mathcal{H}(z_n)}
\end{equation}
に従うことが知られているので，
\begin{equation}
0=-\sum_{m(\neq n)}\mathcal{L}_{mn}\ p_{\mathrm{eq}}(z_n)+\sum_{m(\neq n)}\mathcal{L}_{nm}\ p_{\mathrm{eq}}(z_m)
\end{equation}
が成り立ちます．上式はつりあい条件と呼ばれます．移項して和の記号の中身を比較すると
\begin{align}
\frac{\mathcal{L}_{mn}}{\mathcal{L}_{nm}}=\frac{\phi_{n\to m}}{\phi_{m\to n}}=\exp\{-\beta [\mathcal{H}(z_m)-\mathcal{H}(z_n)]\}
\end{align}
が成立します．この式は「詳細つりあい条件」と呼ばれます．
3. 古典アニーリングのアルゴリズム
以下に具体的なシミュレーティッド・アニーリングのアルゴリズムを示します．
3.1. MetroPolis法
- スピン状態を初期化．$\ \beta=\beta_{\text{min}}\,$と設定．
- 以下の受容率$\ A(\sigma_k\ \to-\sigma_k)\,$を計算．
\begin{align}
\color{red}{
A(\sigma_k\ \to -\sigma_k)=\mathrm{min}\left(1, e^{-\beta\,\left[\,\mathcal{H}(-\sigma_k)\,- \,\mathcal{H}(\sigma_k)\,\right]}\right)
}
\end{align}
3. $\,[0,1]$間の一様乱数$\,r\,$を発生させ，$\,r\leq A(\sigma_k\ \to -\sigma_k)\,$なら，スピン反転操作をコミット
4. 逆温度$\beta$を上げてステップ２. へ戻る．$\,\beta=\beta_{\text{max}}\,$まで到達した場合は終了．
上にアルゴリズムをまとめました．$-\sigma_k\,$は「ランダムに選択された$\,k\,$番目のスピンを反転させる」 操作になります．受容率の計算部分で，詳細つりあい条件が出現していることが分かります．この式中の$\ \mathrm{min}(1,x)$は１と$\,x$で小さいほうを取る演算子です．これは，乱択されたスピンを反転した後の状態が，前の状態に比べてエネルギーが小さくなる場合，後の状態に確実に(確率１で)遷移させるため，このような形になっています．
3.2. 熱浴法
熱浴法のアルゴリズムは，Metropolis法のアルゴリズムのうち，受容率の計算部分が異なります．具体的には，受容率の式を
\begin{equation*}
A(\sigma_k\to-\sigma_k)=\frac{e^{-\beta\mathcal{H}(-\sigma_k)}}{e^{-\beta\mathcal{H}(-\sigma_k)}+e^{-\beta\mathcal{H}(\sigma_k)}}
\end{equation*}
と差し替えれば，熱浴法のアルゴリズムになります．上式も，もちろん詳細つりあい条件を満たしています．
まとめ
本記事では，次の項目についてまとめ，シミュレーティッド・アニーリング(SA)のアルゴリズムを踏襲しました．
- スピン変数&ハミルトニアンの導入
- 統計力学的における確率の記述
- MetroPolis法&熱浴法
次回は，本記事で紹介したSAのアルゴリズムをもとに古典アニーリングマシンの実装を行います．