投稿: 哈密顿蒙特卡洛算法为什么会发散或停止 徐晓鹏

content/post/2020-05-26-hmc.md

@@ -0,0 +1,320 @@
+---
+title: 哈密顿蒙特卡洛算法为什么会发散或停止
+author: 徐晓鹏
+date: '2020-05-26'
+slug: hmc
+categories:
+  - 统计计算
+  - 机器学习
+---
+
+本文讨论哈密顿蒙特卡洛算法的稳定性问题及其解决方法。
+
+这里的主要思想是：一个哈密顿仿真的总能量守恒，则仿真是稳定的。
+
+比如如下这个仿真，红色的是位置-势能轨迹，黑色是位置轨迹，而蓝色是势能（对应概率）的三倍标准差界限。
+如果总能量守恒，则最大势能是一个常数，则仿真可以持续进行下去。
+
+![总能量守恒的哈密顿仿真](https://statisticalcomputing.github.io/ltximg/hmc5.gif)
+
+反之,如果总能量不守恒:
+
+1.  总能量持续变大，则最大势能也会持续变大，导致仿真发散;
+
+![总能量增大的哈密顿仿真](https://statisticalcomputing.github.io/ltximg/hmc5_1.gif)   
+
+1.  总能量持续变小，则最大势能也会持续变小，导致仿真停止。
+
+![总能量减小的哈密顿仿真](https://statisticalcomputing.github.io/ltximg/hmc5_2.gif)
+
+因此保持能量守恒对于哈密顿蒙特卡洛方法的正确和稳定地运行非常重要。
+
+如果针对一个系统（比如一群气体粒子），将总能量守恒原则进行展开，则意味着：
+
+1. 如果系统的势能之和为常数，则系统的动能之和为常数。
+
+1. 如果系统的势能之和与位置有关，则系统的动能之和也与位置有关。
+
+在统计计算中，我们一般研究第二种情况。这意味着不同位置对应的动能之和不同，对应的动量分布也不同。因此在仿真算法中，如果在不同位置，以相同分布（比如标准正态）生成动量初始值，那么就破坏了能量守恒原则，也违反了统计（力学）的基本假设。
+
+受到krauth的statis  tical mechanics algorithm and computation的启发，我们采用扔皮球的方式进行讨论。
+
+# Metropolis
+
+## Metropolis算法
+
+此时从位置a扔一个皮球，假设皮球只能落在位置a,b,c，而且皮球将等可能地落在这三个位置上。
+
+系统处于a,b等位置的静态概率
+
+$$\pi(a),\pi(b),\dots$$
+
+假设从位置a移动到位置a,b和c以后，根据一个概率决定是否接受这个移动。
+
+![Metropolis的假设](https://statisticalcomputing.github.io/ltximg/mh.png)
+
+从位置a到位置b的接受概率，等等
+$$p(a\to b),p(a\to c),\dots$$
+
+假如从位置a跳转到位置a,b,c，那么有：
+$$p(a \to a)+p(a \to b) + p(a \to c) =1$$
+
+考虑到位置a只能从位置a,b,或者c跳转而来:
+$$\pi(a)=\pi(b)p(b \to a) + \pi(c)p(c \to a) + \pi(a)p(a \to a)$$
+
+得到：
+$$\pi(a)(1-p(a \to a) = \pi(b)p(b \to a) + \pi(c)p(c \to a)$$
+
+由于
+$$1-p(a\to a)=p(a \to b)+p(a \to c)$$
+
+代入上式，得到
+$$\pi(a)p(a\to b)+\pi(a)p(a\to c)=\pi(c)p(c\to a)+\pi(b)p(b\to a)$$
+
+从而得到详细平衡方程
+$$\pi(a)p(a \to b)=\pi(b)p(b \to a)$$
+
+$$\pi(a)p(a \to c)=\pi(c)p(c \to a)$$
+
+根据详细平衡方程，可以得到Metropolis算法，即跳转的接受概率与跳转前后位置的静态概率有关，具体讨论请见krauth的书。
+$$p(a \to b)=\min[1,\frac{\pi(b)}{\pi(a)}]$$
+
+
+## Metropolis-Hastings算法
+
+如果位置之间的跳转概率不等，那么总的转移概率$\mathcal P(a \to b)$分为两部分，即从皮球从位置a出发，落在b的概率，以及接受这个跳转的概率。
+$$\mathcal P(a\to b)=\mathcal A(a \to b)\cdot p(a \to b)$$
+
+
+![Metropolis-Hastings的假设](https://statisticalcomputing.github.io/ltximg/mh-with-prior.png)
+
+这时候的细致平衡方程为
+$$\pi(a)\mathcal P(a\to b)=\pi(b)\mathcal P(b \to a)$$
+
+于是得到
+$$\frac{p(a\to b)}{p(b\to a)} = \frac{\pi(b)}{\mathcal A(a\to b)} \frac{\mathcal A(b \to a)}{\pi(a)}$$
+
+于是得到广义的Metropolis算法:
+$$p(a\to b)=\min[1,\frac{\pi(b)}{\mathcal A(a\to b)}\frac{\mathcal A(b \to a)}{\pi(a)}]$$
+
+其中$\mathcal A(a\to b)$是从位置a移动到位置b的概率,而$p(a\to b)$是接受这个移动的接受概率。
+
+## 两种Metropolis算法之间的关系
+
+显然如果位置a和b之间的相互跳转概率相同（满足如下公式）：
+$$\mathcal A(a \to b)=\mathcal A(b \to a)$$
+
+，则广义Metropolis算法会简化为Metropolis算法。
+特别需要指出的是，在哈密顿动力学中，如果两个状态能够相互转换，并且总能量相同，则二者的相互跳转概率亦相同。因此能量守恒的哈密顿系统可以应用Metropolis算法来计算状态的概率。
+
+
+# 哈密顿蒙特卡洛算法
+
+## 哈密顿动力学
+
+哈密顿系统中，总能量等于动能和势能之和：
+$$H=K+U$$
+
+哈密顿公式指出了变量对时间的导数和能量对变量的偏导数之间具有如下关系
+$$\dot x=\frac{dK}{dp}$$
+
+$$\dot p=-\frac{dU}{dx}$$
+
+
+
+## 存在的问题
+
+假设未知概率 $\pi(x)$ 参数的维度很高，则普通的蒙特卡洛方法无法进行抽样。
+这时可以其用势能来表示此概率：
+$$U(x)=-\log \pi(x)$$
+
+动能一般取二次函数:
+$$K(p)=-\frac{p^T p}{2}$$
+
+通过哈密顿仿真得到轨迹，在轨迹上使用蒙特卡洛进行抽样，可以获得  $\pi(x)$ 的参数。
+以下是哈密顿蒙特卡洛算法[information theory, inference and learning,p388]。
+该算法以下有两方面的问题:
+
+1.  接受概率公式
+2.  能量守恒
+
+，我们逐一进行讨论。 
+
+![img](https://statisticalcomputing.github.io/ltximg/hmc.png)
+
+
+### 接受概率
+
+算法第5行到第7行运行哈密顿仿真。算法第3行为每次仿真前的总能量，第8行为仿真后的总能量。
+
+由于哈密顿仿真保持总能量不变，因此（从原则上）仿真前后的总能量相同。
+
+根据以上算法，可知，接受概率的公式为:
+$$p(x \to x_*)=\min(1,\frac{e^{-U(x_*)-K(p_*)}}{e^{-K(p)-U(x)}})=\min(1,1)=1$$
+
+
+由于接受概率总是1，该算法实际上是在哈密顿轨迹上等概率选取样本点，而没有使用Metropolis算法。
+
+如果认为哈密顿与蒙特卡洛属于不同的技术，那么哈密顿蒙特卡洛采样器实际上应当称为哈密顿采样器，因为它只使用了哈密顿技术。
+
+### 能量守恒
+
+哈密顿蒙特卡洛算法有一个特点，即非常容易发散或者停止运行。
+
+发散和停止现象是由于该算法破坏了能量守恒法则。算法的第三行会重新生成一个随机的动量，这时会改变动能，因此总能量也随之改变了。
+
+不断在不同仿真之间改变总能量可能导致总能量发散或者消失现象，这也是该算法不稳定的根本原因。
+
+
+# 更正后的HMC
+
+如果将哈密顿蒙特卡洛接受概率公式与广义Metropolis公式相比，该算法有两个假设:
+
+1.  $\pi(x)=e^{-U(x)}$
+
+2.  $\frac{\mathcal A(x \to x_*)}{\mathcal A(x_* \to x)} = \frac{e^{-K(p)}}{e^{-K(p_*)}}$
+
+第一个假设是目标假设，可以认为是真。问题出在第二个假设。
+
+为了揭示其问题，考虑一个特殊情况：如果在仿真之间保持总能量守恒，则两个状态的相互跳转概率相同：
+
+$$\mathcal A(x\to x_*)=\mathcal A(x_*\to x)$$
+
+则第二个假设变为以下形式：
+$$1 = \frac{e^{-K(p)}}{e^{-K(p_*)}}$$
+
+但由于每次仿真的初始动量p和末尾动量p\*显然不相同，因此上述假设的错误就非常明显了。
+
+
+## 另一个例子
+
+为了说明其错误，我们还可以举一个更具体的例子。
+假设皮球可以在地面上弹起，并且没有能量损失。
+
+可以设置两个状态：
+
+1.  地面
+2.  弹跳的最高点处
+
+显然这两个状态相互转移:
+
+1.  所需时间相同
+2.  两个状态的总能量相同
+3.  两个状态跳转概率相同，因为是依次跳转
+4.  两个状态动量不同，动能也不同
+
+则哈密顿第二个假设的左边仍然为1，右边仍然不等于1，再次显示这个假设有问题。
+
+因此在不同仿真之间保持总能量恒定的前提下，正确的接受概率的公式应当为：
+$$p(x \to x_*)=\min(1,\frac{e^{-U(x_*)}}{e^{-U(x)}})$$
+
+我们也可以换一种方式推论。由于两个状态的总能量相同，按照哈密顿蒙特卡洛原来的接受概率公式，两个状态的相互接受概率都是１，因此两个状态的采样概率相同，都是0.5。
+
+然而，根据势能和概率关系的假设（即假设１），势能小的状态采样概率应当更高。因此哈密顿蒙特卡洛产生了错误的采样结果。
+
+为了保持能量守恒，每次重新初始化动量以后，应当重新设置其向量半径。
+
+根据恒定总能量和当前的势能，可以算出所需的动能：
+$$K(p)=\frac{p^T \cdot p}{2} = \frac{r^2}{2}=H-U(x)$$
+
+
+由此得到动量向量的半径应当为：
+$$r=\sqrt{2(H-U(x))}$$
+
+
+因此动量应当在半径为 $r$ 的超球上均匀分布。
+
+顺便提一下，对于N个粒子总动能恒定的系统，物理学家或者数学家称之为（半径为r的）3N-1球 $\mathbb{S}_R^{3N-1}$ 。
+
+综上所述，得到更正后的算法：
+![遵守能量守恒的哈密顿蒙特卡洛算法](https://statisticalcomputing.github.io/ltximg/chmc.png)
+
+
+# 实现代码
+
+本文介绍如何实现最简单的哈密顿蒙特卡洛。
+
+1.  首先定义势能函数及其导数，这里使用了一个二次函数，对应二维正态分布
+
+        SIGMA = np.array([[1,.95],[.95,1]])
+        U = lambda x: np.dot(x,np.dot(np.linalg.inv(SIGMA),x))/2
+        dU = lambda x: np.dot(np.linalg.inv(np.array([[1,.8],[.8,1]])),x)
+
+2.  然后定义动能函数及其导数
+
+        K = lambda p: np.dot(p,p)/2
+        dK = lambda p: p
+
+3.  定义初始位置和动量
+
+        xStar = np.random.randn(2)
+        pStar = np.random.randn(2)
+
+4.  定义步长和步数超参数
+
+        STEPS = 100
+        delta = 0.1
+
+
+
+5.  保持能量守恒
+    为了保持总能量恒定，仿真的每步都重新设置一下动量幅值大小，可参阅[半径计算公式](hmc1.md)
+    1.  记录总能量
+
+            H = U(xStar) + K(pStar)
+
+            H = U(xStar) + K(pStar)
+
+    2.  重新设置动量半径
+
+            pn = np.sqrt(2*(H - U(xStar)))
+            pStar = pStar/np.linalg.norm(pStar,2) * pn
+
+6.  Metropolis
+    1.  接受概率
+
+        接受概率大当且仅当新位置的能量小
+
+            alpha = min(1,np.exp(U(x0)  - U(xStar)))
+
+    2.  Metropolis 算法
+
+            if alhpa<rand():
+                xStar = x0
+                pStar = p0
+
+7.  开始仿真
+
+        for i in rangee(STEPS):
+            xStar = xStar + delta*dK(pStar)
+            pStar = pStar - delta*dU(xStar)
+8.  整个代码
+
+        import numpy as np
+        from numpy.random import rand
+
+        SIGMA = np.array([[1,.95],[.95,1]])
+        U = lambda x: np.dot(x,np.dot(np.linalg.inv(SIGMA),x))/2
+        dU = lambda x: np.dot(np.linalg.inv(np.array([[1,.8],[.8,1]])),x)
+        K = lambda p: np.dot(p,p)/2
+        dK = lambda p: p
+        xStar = np.random.randn(2)
+        pStar = np.random.randn(2)
+        STEPS = 10
+        delta = 0.1
+        #总能量
+        H = U(xStar) + K(pStar)
+        SIMULATIONS  = 100
+        for i in range(SIMULATIONS):
+            #重新设置动量幅值
+            pn = np.sqrt(2*(H - U(xStar)))
+            pStar = np.random.randn(2)
+            pStar = pStar/np.linalg.norm(pStar,2) * pn
+            x0 = xStar
+            for i in range(STEPS):
+                xStar = xStar + delta*dK(pStar)
+                pStar = pStar - delta*dU(xStar)
+            alpha = min(1,np.exp(U(x0)  - U(xStar)))
+            if alpha < rand():
+                xStar = x0
+