Single Flip Metropolis算法的失效

20世纪50年代，随着计算机的发明和发展，Monte Carlo算法推动课了整个计算物理学的开创和发展。当时的统计物理学热点问题，除了非常火爆的量子统计场论，还有复杂性系统和非平衡问题。

Ising模型于1920年提出，1943年Lars Onsager给出了二维零外场模型的解析解，轰动一时。在后来的二十余年中，无数研究者试图解析解三维Ising模型，但是都以失败告终。与此同时，Ising模型的数值模拟快速发展。

Single Flip Metropolis算法是最简单直接的一种算法。如果你熟悉这个算法，并且熟悉Frenkel的记号体系[1]，可以跳过下面的说明。SF算法的大致思想是：

每一个循环中，随机选择一个格点；
随机选择一个移动方向（对于Ising模型来说，只有一种选择就是flip，所以就用不着随机了）；按照这个选择得到的构型成为提议构型；
计算提议构型和原构型的能量差$\Delta E=E_1-E_0$；
Metropolis算法：产生一个$[0,1)$随机数，假如该随机数小于$p=\exp(-\beta\Delta E)$，则提议构型被接受，成为新构型；否则提议构型被拒绝，原构型成为新构型。这一步通常称为“按照概率$p$接受”。

数学原理上我们要证明这个算法能够在统计的样本池中产生符合Boltzmann分布的采样。这里我们假设算法满足细致平衡，来验证概率$p$确实如此表达。首先，依据细致平衡，

$\begin{equation} N(o)\pi(o\to n)=N(n)\pi(n\to o) \end{equation}$

这里$N(o)$和$N(n)$表示样本池中就构型和新构型的数量比例，通常就等价于对应构型的Boltzmann权重。$\pi(o\to n)$和$\pi(n\to o)$是转移矩阵的矩阵元，表示一个条件概率：假如当前处于某一态a，下一步处于另一态b的概率记为$\pi(a\to b)$。

现在Boltzmann权重知道了，求转移矩阵的矩阵元。根据上述算法，处于构型o时，产生试探构型n的概率记为$\alpha(o\to n)$，试探构型被接收的概率记为是$acc(o\to n)$。显然我们待求的$p$就是$acc(o\to n)$。

$\begin{equation} \text{由构型o转移到构型n的概率}=\text{由构型o产生试探构型n的概率}\times \text{接受构型n作为o的下一步的概率}\nonumber \end{equation}$

即

$\begin{equation} \pi(o\to n)=\alpha(o\to n) acc(o\to n) \end{equation}$

那么代入细致平衡的公式，得到

$\begin{equation} \frac{acc(o\to n)}{acc(n\to o)}=\frac{N(n)\alpha(n\to o)}{N(o)\alpha(o\to n)} \end{equation}$

然而任何一步随机选择待翻转的格点的过程都是与构型无关的，所以必然在任何时候都相等，即$\alpha(n\to o)=\alpha(o\to n)$。所以

$\begin{equation} \frac{acc(o\to n)}{acc(n\to o)}=\frac{N(n)}{N(o)}=\exp[-\beta(E(n)-E(o)] \end{equation}$

得到这一步之后，通常按照以下方式给$acc(o\to n)$赋值

$\begin{equation} p=acc(o\to n) = \min\{1, \exp[-\beta(E(n)-E(o))]\} \end{equation}$

仔细分析一步就可以看出这个$p$的取值就是上面给出的接受概率。

这个算法在高温（高于相变温度）时可以得到很好的结果，效率也不错。但是在接近相变点以及相变点以下，就会变得非常缓慢，原因一是因为“临界降速”（critical slowing down），临界点的Single Flip难以达到平衡；二是因为高自由能垒导致的对称性自发破缺，低温下相界面等的自由能垒难以跨越，表现为$acc(o\to n)$极低。为了解决这些问题，人们提出了很多算法，其中最有效的是cluster算法。

Cluster算法核心思想

1987年，Robert Swendsen和Jian-Sheng Wang提出第一个cluster算法。这个算法给出了cluster一类算法的核心思想：

每一步创建一个尝试构型，这个构型被创建的概率（频率）为对应的Boltzmann权重；
该尝试构型被接受的概率为1。

我们可以将这个思想用公式表达，从而探讨创建构型的方法。从细致平衡条件开始

$\begin{equation}\label{eqn:detailed_balance} N(o)P_\text{Gen}^o(\{\text{cluster}\})P^{\{\text{cl}\}}(o\to n) acc(o\to n) =N(n)P_\text{Gen}^n(\{\text{cluster}\})P^{\{\text{cl}\}}(n\to o) acc(n\to o) \end{equation}$

$N(o)$是原构型的Boltzmann权重，$P_{\text{Gen}}^o(\\{\text{cluster}\\})$是在构型o下产生特定的cluster的概率，$P^{\\{\text{cl}\\}}(o\to n)$是将整个cluster移动到指定位置得到提议构型的概率（对于Ising这样的spin-1/2模型别无选择，就是flip），最后$acc(o\to n)$是这个提议构型被接受的概率。

通常可取$P^{\\{\text{cl}\\}}(o\to n)=P^{\\{\text{cl}\\}}(n\to o)$，即只要选定了相同的cluster，来回移动/翻转的概率是一样的。此外，算法要求$acc(o\to n)=1$。这样我们得到了

$\begin{equation} \frac{P_\text{Gen}^o(\{\text{cluster}\})}{P_\text{Gen}^n(\{\text{cluster}\})}=\frac{N(n)}{N(o)}=\exp(-\beta\Delta E)\label{eqn:14.3.3} \end{equation}$

Swendsen-Wang算法

下面的问题是怎么设计算法、求与算法对应的各个概率。
在以上核心思想指导下，Swendsen和Wang针对Ising模型提出了如下的构型创建方法：
由于Ising模型仅在相邻两格点之间有相互作用，那么我们称相邻两格点之间的空间为键位。每个键位有两种状态，为成键和不成键。

如果相邻两格点自旋反平行，那么它们不成键；这两个格点之间的键位成为不可成键键位；
如果相邻两格点自旋平行，那么它们成键的概率是$p$，不成键的概率是$(1-p)$，这里$p$待求；这两个格点之间的键位成为可成键键位；

（显然，系统能量$E$ = (反平行的相邻格点对数目 - 平行的相邻格点对)数目$\times J$ = (不可成键键位数 - 可成键键位数)$\times J$）

对系统所有键位考察确定每个键位成键或不成键之后，得到该构型下的一个划分，这样的一个划分生成的概率为

$\begin{equation}\label{eqn:cluster_o} P_{\text{Gen}}^o(\{\text{cluster}\})=p^{n_c}(1-p)^{n_b} \end{equation}$

其中$n_c$为所有可成键键位中成键的数目，$n_b$为所有可成键键位中不成键的数目。

这样一个划分中，我们称成键的所有格点为一个团簇。系统被划分称许多个团簇，只包含一个格点也可以称为一个团簇。我们就从所有团簇中选一部分移动（翻转）。由于要求相同划分下团簇翻转来回的可逆$P^{\\{\text{cl}\\}}(o\to n)=P^{\\{\text{cl}\\}}(n\to o)$，那么可以设每个团簇被翻转的概率都是0.5。原则上这个概率可以任取，但是太小对构型的改变不够显著，太大则是整体翻转，所以最合适的值是0.5[2]。这样得到提议构型。提议构型被接受的概率为1。

下面，我们求各个反向概率，即从提议构型返回到原构型的各个概率。

首先假设上面的过程之后有$\Delta$个反平行的不可成键键位变成了平行的可成键键位，那么系统的能量就变为

$\begin{equation} E(n)=E(o)-2J\Delta \end{equation}$

这样Boltzmann权重之比就确定了。而反向的团簇翻转概率和正向的相同，只要生成同一个划分。所以就剩下求反向得到与正向相同的划分的概率$P_{\text{Gen}}^n(\\{\text{cluster}\\})$。

为了产生相同的划分，可成键键位中成键的数目依然为$n_c$，这些键都在团簇内部。但是可成键键位总数变化，增加了$\Delta$个，于是可成键键位中不成键的数目变为$(n_b+\Delta)$。这样

$\begin{equation}\label{eqn:cluster_n} P_{\text{Gen}}^n(\{\text{cluster}\})=p^{n_c}(1-p)^{n_b+\Delta} \end{equation}$

由式$\eqref{eqn:14.3.3}$得到

$\begin{equation} \frac{p^{n_c}(1-p)^{n_b}}{p^{n_c}(1-p)^{n_b+\Delta}}=\exp(2\beta J\Delta) \end{equation}$

$\Delta$随构型变化而变化，对于任意$\Delta$，上式都应该成立，于是解得$p$

$\begin{equation} p=1-\exp(-2\beta J) \end{equation}$

综上，Swendsen-Wang算法表述为[2]：

每步确定所有的键位是成键还是不成键，其中可成键键位的成键的概率为$p=1-\exp(-2\beta J)$，从而得到一个划分；
在划分中的每一个团簇有0.5的概率翻转，得到提议构型；
提议构型按照概率1翻转；
所有键位的状态清空。

Wolff cluster算法

从以上的Swendsen-Wang算法，Wolff改进得到以下等价的算法[2]

每步随机选择一个格点，加入到种子栈和团簇中
从种子栈栈出栈一个格点，检查其所有相邻格点，若其自旋与当前出栈的格点平行，且未在团簇中，则按照概率$p=1-\exp(-2\beta J)$加入种子栈和团簇中
若种子栈非空，则重复2
将团簇中的所有格点都翻转。

我自己的改进是将双栈改成单队列：

每步随机选择一个格点加入到团簇队列中作为起始元素，队列head指向起始元素的下一个，tail指向起始元素（整个队列指的是从起始元素到末尾元素head前一个的所有元素，而不是从tail到head的一段）
检查队列tail指向的元素格点的最近邻格点，如果最近邻格点与tail元素格点平行，且未在团簇整个队列中，则按概率$p=1-\exp(-2\beta J)$加入到head（之后head指向加入元素的下一个）
tail指向下一个（刚才指向的元素不删除）
若tail和head指向同一个元素，则生成结束，否则goto 2
将整个队列的元素格点翻转。

对点阵模型通用的Wolff cluster算法

参看文献[3]。

广义cluster算法

见[1]第14章第3节General Cluster Moves。

通常，对于一般的多粒子模型，不会总能找到方法构造接受概率100%的构型。但是cluster算法通常可以用来提高提议构型的被接受概率。

减小cluster的尺寸

有的时候过高的连通性会导致过大的cluster。这时候我们希望牺牲一点接受概率来控制cluster的尺寸。由式$\eqref{eqn:detailed_balance}$开始，继续要求$P^{\\{\text{cl}\\}}(o\to n)=P^{\\{\text{cl}\\}}(n\to o)$，但是不要求$acc=1$，那么得到

$\begin{equation} \frac{P_\text{Gen}^o(\{\text{cluster}\})}{P_\text{Gen}^n(\{\text{cluster}\})}=\frac{N(n)acc(n\to o)}{N(o)acc(o\to n)} \end{equation}$

式$\eqref{eqn:cluster_o}$、式$\eqref{eqn:cluster_n}$继续成立，那么

$\begin{equation} \frac{acc(o\to n)}{acc(n\to o)}=\exp(2\beta J\Delta)(1-p)^\Delta \end{equation}$

其中$\Delta$是这一步cluster翻转中可成键键位增加的数目。通过这个式子，我们就可以调节加入团簇的概率$p$和提议构型的接受概率$acc(o\to n)$寻找平衡。例如，若设定$p=1-\exp(-\beta J)$，相当于在高温下构造团簇；这样可以求出来对应的提议构型接受概率为

$\begin{equation} acc(o\to n)=\min[1,\exp(\beta J\Delta)] \end{equation}$

显然当$\Delta>0$时，提议构型的能量更低，一定被接受；当$\Delta<0$时，提议构型的能量更高。我们可以做一个估计：要求

$\begin{equation} \exp(\beta J\Delta)>\exp(-1)\approx0.368 \end{equation}$

则$\beta J\Delta>-1$，设$J=k_B=1$，那么$-\Delta<T$。由于$T$基本上是个位数，那么$-\Delta$就不能太大。

当然，这只是一种设计，遇到特殊的情况还是需要特殊的设计来处理。

参考文献

Daan Frenkel, and Berend Smit. Understanding molecular simulation: from algorithms to applications. Elsevier, 2001.
E. Luijten. Interoduction to Cluster Monte Carlo Algorithms. http://csml.northwestern.edu/resources/Reprints/lnp_color.pdf
Kent-Dobias, Jaron, and James P. Sethna. “Cluster representations and the Wolff algorithm in arbitrary external fields.” Physical Review E 98.6 (2018): 063306.