如何理解及运用王和朗道算法（Wang and Landau algorithm）？

7.3 有趣的问题，可以答下

7.31 来自拖延症患者的回答（写的这么辛苦，要不要点个关注･*･:≡(　ε:)

一、为什么需要Wang-Landau算法

在统计物理的正则分布中，大家知道微观态的密度函数是 $P_\beta(\omega) = \frac{1}{Z_\beta}e^{-\beta H(\omega)}$ ，我们回忆一下在Ising模型的模拟中用到过的Metropolis算法，在改变当前系统的位形时，需要考虑接受率 $\min\{e^{-\beta \Delta H}, 1\}$ ，当我们的系统从能量较低的区域想要跳跃到能量较高的区域时将会很容易被拒绝位形变换，而往往系统不同的低能区之间有一些高能的区域隔开，这一区域的接受率比较低，因此，这种稳态的束缚将大大限制Monte-Carlo算法对位形空间的遍历和概率权重估计的效率

Metropolis算法： https://www.zhihu.com/question/24334476/answer/333052138

举个例子，在下面的图中，系统的位形想在两个红色区域之间跳跃需要越过高能区，是比较难发生的，因此，经过了很长时间的数值计算，位形也可能只在某一个红色区域内，导致估计的偏差太大

当我们考虑能量（或者其他物理量）时，就会因为能级简并出现一个权重因子， $P_\beta(\varepsilon) = \frac{1}{Z_\beta}e^{-\beta \varepsilon} \int_{\{\omega \in \Gamma : H(\omega) = \varepsilon\} } d \omega = \frac{1}{Z_\beta} g(\varepsilon)e^{-\beta \varepsilon}$

Note. 这里指出这个是因为后面会用到

因为Wang和Landau的论文里使用不同能级之间的转移来定算法的，所以概率密度用的是上面的形式，也就要出现一个权重因子 $g(\varepsilon) = \int_{\{\omega \in \Gamma : H(\omega) = \varepsilon\} } d \omega = e^{S(\varepsilon)}$ ，实际上就是微观态的数量，也就是熵的函数

按照正常的Metropolis算法，并且不考虑温度（能级差）的影响时，接受率应该是 $\alpha(\varepsilon_1 \rightarrow \varepsilon_2) = \min\{\frac{g(\varepsilon_2)}{g(\varepsilon_1)}, 1\}$ ，同样无法避免从态密度很高的区域到态密度很低的区域将导致的拒绝率过高的问题。这时候，Wang, Landau想到的引导采样系统进入这些低态密度区的办法就是把这个接受率倒过来，变成 $\alpha(\varepsilon_1 \rightarrow \varepsilon_2) = \min\{\frac{g(\varepsilon_1)}{g(\varepsilon_2)}, 1\}$ ，这样，当目前所处的状态是高态密度区时，下一步将会以比较高的接受率保持在高态密度区或者进入低态密度区。通过这样来使得采样系统尽可能遍历所有位形足够多次，提高MC算法的效率和精度

Note. 原论文中的原话也是算法基于把接受率倒过来就能通过平坦直方图实现能量分布的发现

Our algorithm is based on the observation that if we perform a random walk in energy space by flipping spins randomly for a spin system, and the probability to visit a given energy level E is proportional to the reciprocal of the density of states 1/g (E ), then a flat histogram is generated for the energy distribution.

所以，在数学上证明了这个算法的收敛性之前还是一个比较神奇的启发式算法，对参数的选取有一些技巧，最近几年已经有工作证明过收敛性

二、怎么实现Wang-Landau算法

我们假设我们用到的提议分布是对称的 $P(x,y)$ ，这样可以保证接受率是上面的形式

1. 首先对能谱进行剖分，分为N个区域，那么间隔 $\Delta = \frac{E_{max}-E_{min}}{N}$ ，区间 $E_i = [E_{min}+(i-1) \Delta, E_{min}+i \Delta), i=1,2,..., N-1$ , $E_N = [E_{max} - \Delta, E_{max}]$ ，映射 $I$ 把能量映射到对应的区间下标

2. 初始化

（1）由于没有先验信息，将各个区域的初始态密度函数设置成1，为了避免溢出，取对数用熵来计算， $S_i = 0, i=1,2,...,N$

（2）步长f=1

（3）初始化位形 $\omega_0 \in \Gamma$

3. while(f > epsilon && t < max_iteration)进入迭代

（1）提议新的位形 $\omega_{t+1} \sim P(\omega_t, \cdot)$

（2）计算接受率 $\alpha = \min\{\exp(S_{I\circ H(\omega_t)} - S_{I\circ H(\omega_{t+1})}), 1\}$

（3）是否接受 $\mathrm{if(rand < \alpha) ~ accept;}$

（4）修改当前位形和数据， $\mathrm{Hist}_{I\circ H(\omega_{t+1})} ++; ~ S_{I\circ H(\omega_{t+1})}+=f; ~ t++; ~ \omega_t = \omega_{t+1};$

（5）若直方图平坦，则修改步长 $\mathrm{if(Hist.is\_flat()) ~ Hist = 0; f*=0.5;}$

三、Wang-Landau算法只能用在物理上吗

当然不一定，如果你真正理解了Wang-Landua算法的思想，这个算法是一个非常有用的重要性采样方法

举个例子，在高维数值积分中的应用

函数 $y:\mathbb{R}^{n} \rightarrow \mathbb{R}, x \mapsto y(x)$ ，计算积分 $I = \int_a^b y(x) dx$

那么这里的位形空间就是 $[a, b] \subset \mathbb{R}^n$ ，我们还是用前面的Note. 那里提到的，先用一个低维的宏观量来描述系统，很明显， $y$ 就很合适，因此，

$I = \int_{y_{min}}^{y_{max}} y g(y) dy$ ，其中， $g(y) = \int_{\{x \in [a, b] : y(x) = y\}} dx$ 是态密度函数

因此，问题就转化成了 $I = \sum_{i=1}^N g_i \frac{y_{i+1} + y_i}{2}$ ，其中， $g_i = \int_{y_{i-1}}^{y_i} g(y) dy$ ，用Wang-Landau算法把态密度函数估计出来就可以了

实际上，Wang-Landau算法作为自适应MC算法，配合变分贝叶斯，在机器学习里似乎也有应用

四、Wang-Landua算法的原理

我们把微观态的概率测度写成 $P(dx) = \pi(x) \lambda(dx)$ ，这里的 $\lambda$ 是定义在相空间 $\Gamma$ 上的某个测度，比如连续情形的Lebesgue测度， $\pi(x) = e^{-\beta H(x)}$

选定一个低维的宏观量（光滑函数） $\xi : \Gamma \rightarrow [a, b]$ ，利用这个宏观量对相空间做一个剖分， $X_i = \xi^{-1}([a_{i-1}, a_i)), i = 1, 2, ..., N-1$ , $X_N = \xi^{-1}([a_{N-1}, a_N])$

我们要估计定义在剖分上的分布函数 $\theta_*(i) = \int_{X_i} \pi(x) \lambda(dx)$ ，也就是在 $\Theta = \{\theta\in\mathbb{R}^N: \theta(i) \geq 0, \forall i=1,2,...,N ~ \mathrm{and} ~ \sum_{i=1}^N \theta(i) = 1\}$ 中找到一列 $\{\theta_n\} \stackrel{a.s.}{\rightarrow} \theta_*$

我们希望找到一个阶梯形势场 $\{V_1, V_2, ..., V_N\}$ ，使得在这个势场的共同作用下，每一个分块上的平均的概率测度都相同 $\frac{1}{N} = \int_{X_i} \pi(x) e^{- \beta V_i} \lambda(dx) = e^{- \beta V_i} \theta_*(i) \Rightarrow e^{- \beta V_i} = \frac{1}{N \theta_*(i)}$

因此，当概率测度变成 $P_{\theta_*}(dx) = \sum_i \pi(x) e^{-\beta V_i} \lambda(dx) \mathbb{I}(x \in X_i) = \frac{1}{N} \sum_i \frac{\pi(x)}{\theta_*(i)}\lambda(dx) \mathbb{I}(x \in X_i)$

MC求和公式变为重要性采样的形式 $\langle A \rangle = \int_\Gamma A(x) \sum_i e^{\beta V_i}\mathbb{I}(x\in X_i) P_{\theta_*}(dx)$

Note. 到这里，我们说明了一个分布函数相当于给出了一个势场，并且通过得到新的哈密顿量来构造新的正则分布，即通过 $\theta_t$ 所给出来的势场改变的正则分布为

$\pi_{\theta_t}(x) = \frac{1}{N} \sum_i \frac{\pi(x)}{\theta_t(i)}\mathbb{I}(x \in X_i) = e^{-\beta(H(x) + \sum_i V_i \mathbb{I}(x \in X_i))}$

到现在为止都是重要性采样的套路，那么算法的困难之处在于 $\theta_*$ 是未知的，因此，Wang-Landua算法解决的问题就是如何在这个未知的 $\theta_*$ 中采样，即产生采样序列 $\{(\theta_t, x_t)\} \subset \Theta \times \Gamma ~ s.t. ~ \lim_{t \rightarrow + \infty} \theta_t = \theta_*$ 或者说找到势场 $\{V_i\}$ 使得各个分块上的块哈密顿量相同

用Metropolis算法更新位形 $x_t \sim P_{\theta_t}(x_{t-1}, \cdot) s.t. \pi_{\theta_t}P_{\theta_t} = \pi_{\theta_t}$

用某一函数更新分布函数 $\theta_{t+1} = \theta_t + \gamma_{t+1}F(x_{t+1}, \theta_t)$ ，使得当 $x_{t+1} \in X_i$ 时， $\Delta \theta(i) >0,~\Delta \theta(j)<0,\forall j\neq i$

来源：知乎 www.zhihu.com

作者：董玉龙

【知乎日报】千万用户的选择，做朋友圈里的新鲜事分享大牛。
点击下载

此问题还有 6 个回答，查看全部。
延伸阅读：
如何用简单易懂的例子解释隐马尔可夫模型？
哪些控制类的算法惊艳了你？