兔子走遍世界——物理学视角下的采样方法

曾经有过一个「有志气的兔子」的故事，相信很多朋友都听过：为了找出地球上最高的山，一群有志气的兔子们开始想办法：

兔子朝着比现在高的地方跳去。他们找到了不远处的最高山峰。但是这座山不一定是珠穆朗玛峰。这就是局部搜索，它不能保证局部最优值就是全局最优值。

兔子喝醉了。他随机地跳了很长时间。这期间，它可能走向高处，也可能踏入平地。但是，他渐渐清醒了并朝最高方向跳去。这就是模拟退火。

兔子们吃了失忆药片，并被发射到太空，然后随机落到了地球上的某些地方。他们不知道自己的使命是什么。但是，如果你过几年就杀死一部分海拔低的兔子，多产的兔子们自己就会找到珠穆朗玛峰。这就是遗传算法。

兔子们知道一个兔的力量是渺小的。他们互相转告着，哪里的山已经找过，并且找过的每一座山他们都留下一只兔子做记号。他们制定了下一步去哪里寻找的策略。这就是禁忌搜索。

这几个比喻非常经典，经典到我已经不太能确定这个比喻最早的出处。近些年来，人工智能、机器学习、大数据成为了大家街头巷尾热议的主题，上面提到的这些算法得到了非常广泛的注意。在许多实际问题中，我们经常会需要找到「最少的投入」的方案，「能量最低」或者「熵最大」的状态，「效率最高」的某种设计，在这些问题中，上面所介绍的一些算法，再配合上一些改进，可以帮助我们以较小的代价，较快的速度，找到「地球上最高的山」，即找到体系在满足一些约束条件的情况下，所能取得的最小值或者最大值。

不过仅仅只知道最值也还可能并不够。考虑能量的极小化，我们可以用上面介绍的这些方法求得一个最低的能量，虽然我们知道「能量最低」的状态通常是最稳定的状态，但是因为各种随机因素的影响，系统很可能会偏离这样的能量最低状态。假如体系还有若干个跟最小值很接近的次小值（这在玻璃体系里是常见的），那么稍稍扰动一下（提高一点温度），就会有一些粒子被激发到其它次低的能量状态上去，在测量「平均值」的时候就会必须考虑到那些次小值的影响，随着外界扰动大小（温度）的增加，系统会越来越偏离能量最低的状态。因此，当我们希望能预测出在各种不同温度（或者任何其它参数）条件下系统的宏观状态，单纯地去寻求「能量最低」就很可能是不够的了，如果不能对系统的所有可能取值的情况有所了解，就无法对各种不同环境下系统的一些性质进行估计。因此，有时兔子不只要「为了找出地球上最高的山」，还会需要知道整个地球的地形。

第一类方案是最传统的 Monte Carlo，我们可以直接拿来一个摄像机，对着系统拍摄，如果系统的状态发生了改变，我们就记录下系统的状态以及对应的能量，等过了足够长的时间，系统应该就遍历了所有的状态，这时我们就可以进行统计了。但很明显，「足够长的时间」实在太可怕，只要系统的温度不太高，那么系统演化的速度就可能很慢，在很长的一段时间里，系统可能都陷在某个能量极小的状态里。我们继续用之前的兔子的比喻，如果说求解态密度是需要兔子对全地球的地形情况有所了解的话，那么这种方案就是我们跟着一只兔子或者有限的几只兔子（并行）的脚步去了解整个地形情况，一旦兔子陷入到某些坑里，我们就需要很长的时间来等待兔子们从坑里爬出来，这一方案显然会是低效的。

第二类方案是副本交换的一些方法，我们考虑到之所以需要等待「足够长的时间」，主要是因为系统陷在某个能量极小态里，于是我们考虑用「高温」的模拟，来「加速」系统的演化。例如对于同一个体系，我们可以在 200K, 250K, 300K, ..., 550K, 600K 等不同的温度下进行模拟或者采样，在实际操作中，这一计算（ Parallel tempering ）是用并行的方法来进行的，低温体系中的构象可以与高温体系中的构象按照一定的概率 p 发生交换：

这一形式可以看成是传统的 Metropolis 的「升级版」。如果我们用之前兔子的比喻，对应于高温体系的兔子是磕了药且开了挂的兔子，它们可以在各类不同地形的区域不停狂奔，而我们隔一段时间就按照一定的概率把这种兔子跟（在低温情况下）可能陷入到坑里的兔子交换一下状态，即给掉坑里的兔子灌药开挂，告诉原本开了挂的兔子「差不多行了，正常一点」，让它们回复正常状态……利用这种方式，我们就可以相对更高效地得到系统在不同温度情况下的各统计量。值得一提的是，这里的「温度」还可以替换成其它描述环境的变量，例如 pH 等。如果不希望进行并行计算，那么可以考虑让体系在「温度」这一坐标上进行随机行走，关于这一方法的介绍可以参考论文： Enhanced sampling in generalized ensemble with large gap of sampling parameter: Case study in temperature space random walk 。

第三类方法并不限定体系处在某几个特定的温度上，而是直接关注态密度的计算。在温度为 T 的情况下，系统的平均能量。这个式子里面的i下标代表不同的能级，即为第i能级的能量。因为即使能量相同，系统的状态也可能不同，因此需要引入一个「态密度」，它的大小反映了在同一能量的情况下，系统可能处在多少种不同的状态，之后的 Boltzmann 因子正比于我们通常计算随机变量期望值时候的概率。这里其它的概念大家都比较熟悉，唯一觉得有些难以理解的就是态密度，我们用一个非常简单的实例来说明这一概念：

假如我们研究的体系一共有 N 个空位，有 n 个小球（n。小球的例子显得非常简单，因为我们很容易地就把态密度解了出来，但更多复杂的实际问题里，态密度并不能直接用组合的方法求得。但为了求得系统的平均能量，态密度又是不能不求解的，因而我们需要设计一些方法把态密度求出来。这一类的方法称为多正则系综（ Multicanonical ensemble ），其中最常用的就是Wang-Landau 算法了。

如果上面的例子还是觉得有些难以理解，那就只能用兔子的比喻了。地球表面均匀分布了兔子，因为兔子所在的高度不同，能量也就不同。但因为兔子已经均匀分布在地球表面，我们只需要统计在某个确定的海拔高度（能量）上放了多少个兔子，就可以知道地球表面的地形情况，某个高度上的兔子数就是「态密度」。Wang-Landau 算法考虑：如果某一个高度上已经放了兔子了，那么我们就降低下次在这个高度上放兔子的概率，最终，在所有的高度上放兔子的概率变得差不多的时候，这时每个高度上的兔子数应该就等于对应高度上的「态密度」。有了态密度，其余所有的统计物理学量也就可以求解了。

第四类方法是给粒子加上一些人为的「力」，让粒子可以跨越能垒，最终让粒子可以跑遍整个能量面，而在计算统计量的时候，我们再采用一些分析方法把这些人为的因素扣除掉。

这类方案中被广泛使用、且已经有了一些历史的方法就是伞形抽样（ Umbrella sampling ）的方法。继续用兔子的比喻，为了让兔子走遍世界而不浪费太多的时间在各种深坑里，人为地在每个兔子屁股后面加装了一个可以朝地面喷气的装置，这种装置对跨越能垒能起到很大的帮助，最终，我们再考虑用一些分析方法（例如Weighted Histogram Analysis Method ）把这些喷气装置的影响来消除，于是就可以计算得到各个统计量。

Metadynamics 是新世纪里出现的一种方法。勤劳的兔子们希望了解地球表面的地形，但它们很担心自己会在一些坑里浪费太多的时间，于是兔子们考虑来把地球填平（如题图中的下图所示），它们在坑里填上泥巴，这样它们终于到达了地面，于是可以随意奔跑，直到遇到另一个坑，继续填……最终，地球表面所有的坑都填平了，兔子们可以随意在一个光滑的地球表面奔跑了——这时只要把兔子们往各类坑里填上的泥巴全部翻出来，我们就得到了地球表面的地形。大家可以到这个页面来查看这一算法具体实现的技巧以及演示的动画： http://people.sissa.it/~laio/Research/Res_metadynamics.php。

更多采样的方法可以参考介绍这些方法的论文，这里不再一一列举。这里给出一些中文的材料供参考：

分子模拟中的增强抽样方法

杂谈*能计算，PMF，伞形抽样，WHAM

AlchemistryWiki索引

Simulation Acceleration

值得一提的还有 Metropolis 算法的很多发展，关于现在在采样方面被广泛使用的 Markov Chain Monte Carlo （MCMC）方法可以参考@茉茉老师提供的文章： http://stat.wharton.upenn.edu/~stjensen/stat542/lecture14.mcmchistory.pdf。下次我可以再写一些东西把一些偏数学的描述与一些物理学的形式联系起来讨论这些问题。