现代Monte Carlo算法的起源:Metropolis算法
Monte Carlo算法的思想起源很早,一个很著名的例子是18世纪的布丰投针实验。人们很早就发现Monte Carlo算法能够用于求数值积分。
统计物理中的Monte Carlo算法始于20世纪40年代,由波兰数学家斯塔尼斯拉夫·乌拉姆(Stanisław Marcin Ulam)提出,用来数值求解曼哈顿计划中中子散射问题的微分方程。1945年10月,第一台通用电子计算机ENIAC诞生,乌拉姆将这一思想告知冯·诺依曼(John von Neumann),随后冯·诺依曼将Monte Carlo算法编写成程序,利用计算机求解。而Monte Carlo这个名字,则由尼古拉斯·梅特罗波利斯(Nicholas Metropolis)提出,来源与乌拉姆的叔叔经常去赌博的摩纳哥蒙特卡罗赌场。
1953年,梅特罗波利斯与罗森布卢特(Rosenbluth)夫妇、泰勒(Teller)夫妇共同发表论文Equation of State Calculations by Fast Computing Machines[1],提出了一种利用计算机求解强关联多体系统的统计性质的随机算法。在这篇文章中,他们计算了二维硬球(rigid sphere)模型。文章核心思路是将正则系综的相空间看做概率空间。如果给定N个粒子的坐标,那么可以求这一构型的能量(势能)。因为势能只与位置有关,与速度无关,于是相空间积分当中的对动量自由度积分便可约去,所以系统的统计性质只与势能部分有关
上式中d表示维度,N表示粒子数量,q表示(广义)坐标,A为待求物理量,$\la A \ra$表示对对物理量A的统计平均。
一个直观的思路是随机给定粒子坐标,求这一构型的物理量A及能量E,然后给定其统计权重为$e^{-E/K_BT}$,最后求物理量A的加权平均。这一方法对稀疏系统可行,但是对稠密系统而言,由于绝大多数构型的权重极小(由于重叠而使得能量极大),而权重大的构型难以被采样,算法效率会极其低下。作者提出了修正方法,即举世闻名的Metropolis算法:与其随机均匀地抽取构型再赋以权重$e^{-E/k_BT}$,不如抽取概率为$e^{-E/k_BT}$的构型再按照平均。具体地,
- 随机地初始化N个粒子的坐标,求其能量。
- 随机选择一个粒子,将其移动$q_i\to q_i+\alpha\xi_i$,这里$\alpha$是步长,可以人为设定,$xi_i$是一个$[-1,1]范围内均匀分布的随机数。为了减小界面效应,作者采用了周期性边界条件。
- 求新构型与原构型的能量差($\Delta E$)。如果新构型能量更低,则新构型被“接受”(Accepted)。如果新构型能量更高,那么取$[0,1]$均匀分布的随机数$\varepsilon$,若其小于$e^{-\Delta E/k_BT}$,则新构型被“接受”;否则被“拒绝”(rejected),即恢复到本次移动之前的构型。
- 【重要!】无论新构型是否被接受,经过2、3两步之后的构型被视为一个抽样构型。
- 重复步骤2-4多次,直到抽样构型足够多以控制统计误差到接受范围内。样本容量记为$M$。
- 对待求物理量做无权重的样本平均,平均值就是一次“试验”的结果。有时会需要求物理量的高阶矩,同样这个高阶矩是一次“试验”的结果。
以上就是Metropolis算法,其中第4步是其核心。显然,这是一个Markov过程,每一次转移(步骤2-4)得到一个状态,从第$n$个状态得到第$n+1$个状态的转移概率$P(X_{n+1}=x|X_{n}=x_n)$与$n$无关。
我们可以进行多次“试验”,得到多次“试验”的结果,再取平均得到物理量的计算结果;对多次“试验”结果求样本标准差得到物理量的计算误差。这是一个完整的计算过程。
对Metropolis算法采得的构型权重都正比于$e^{-E/k_BT}$,即
作者证明如下:
首先,一个粒子允许在整个空间中移动,只要步数足够多,那么这个粒子一定能够到达空间的任意位置。这就是说,算法满足遍历性(ergocity)。当然,这只是申明算法能够充分采样,并不需要在实际计算过程中遍历整个空间。
其次,需要证明如果样本容量足够大,样本池中状态为r的样本数量$\nu_r$应当正比于$e^{-E_r/k_BT}$。关于这点,首先明确这个过程中涉及到
- 先验(a priori)概率$P_{rs}$:当构型为s时,选择下一个构型为r的(条件)概率。显然$P_{rs}=P_{sr}$,因为粒子从位置1能够移动到位置2,那么单纯的看空间位置,当然也能从位置2移动到位置1,而且在步骤2中,选择移动位置是均匀随机的。
- 接受(acceptance)概率$\alpha_{rs}$:即从状态r到状态s能够被“接受”的概率。若$E_r>E_s$,那么$\alpha_{rs}=1$,$\alpha_{sr}=e^{-(E_r-E_s)/k_BT}$。
这样,Markov链的转移矩阵为
由随机过程理论可知,当Markov过程达到平稳时,即采样充分时,满足细致平衡条件(detailed balance)
即
证毕。
这里特别说明:之所以步骤4非常重要,即必须在移动被拒绝的时候将原构型加入样本,可以通过计算$p_{r\to r}$看出:由于转移矩阵的性质$\sum_s p_{r\to s}=1$,
即:状态r回到自身,不仅可以通过直接在第2步选择自身并直接移回($e^{-\Delta E/k_BT}=1$),还能够通过第2步选择其他状态但是被拒绝之后移回。
Metropolis算法的局限性
Metropolis算法能够保证在足够多的Monte Carlo采样下得到足够精确的物理量,但是也有重要的不足:
- 一般地,这个算法不能控制达到足够精度的步数,或者说不能告诉我们样本容量需要达到多大才能满足精度要求。达到足够的精度又称为“收敛”(converge)。Markov过程需要达到平稳才能保证抽样的正确分布,也就是说通常需要弃掉一开始的多个构型使得系统达到“平衡”,再进行采样。但是系统需要多久才能达到“平衡”并非已知。有的时候(例如在临界相变点附近),会发生临界降速(critical slowing-down),达到平衡的时间变为无穷长。
- 对于连续空间的Metropolis算法,步长$\alpha$取值尚未确定。显然,$\alpha$与收敛速度有关,如何调节$\alpha$值以使Monte Carlo算法最有效率,是一个很重要的问题。
经过维尔弗雷德·哈斯廷(Wilfred Keith Hastings)的总结推广[2],Metropolis算法中接受概率可以不由波尔兹曼分布得出,而可以是任意要求的概率。这样推广的Metropolis算法成为一种通用的计算数学方法,称为Metropolis-Hastings算法。哈斯廷在总结中还讨论了Metropolis算法的误差估计等。
Metropolis算法是Markov Chain Monte Carlo(MCMC)算法的一种,也是重要性采样(Importance Sampling)方法的一例。除了Metropolis-Hastings算法外,由Metropolis算法还可以导出处理高维采样、计算联合分布的Gibbs抽样法。这些方法在计算统计学、机器学习等领域中有很重要的应用。关于这些算法,可以参考:
随机采样和随机模拟:吉布斯采样Gibbs Sampling
Gibbs采样
参考资料
- Equation of State Calculations by Fast Computing Machines, Metropolis N., Rosenbluth A. W., Rosenbluth M. N., Teller A. H. and Teller E., Journal of Chemical Physics, 21, 1087-1092 (1953) PDF
- Hastings, W. Keith. “Monte Carlo sampling methods using Markov chains and their applications.” (1970): 97-109.
【持续更新】