Hamiltonian Monte Carlo抽样算法的初步理解

吐血学习
HMC 算法算是代码实现没几行，背后原理让人吐血的一个算法。

接受拒绝采样算法

我们的目标分布 $\tilde{p}$解析形式已知，例如
$$\tilde{p}(z)=0.3exp(-(z-0.3{)}^2)+0.7exp(-(z-2{)}^2/0.3)$$
我们要面临的问题是如何从上述复杂分布中采样？

拒绝采样算法通过使用一个已知的易于采样的分布 $q(z)$ 去逼近我们要采样的目标分布$\tilde{p}(z)$，

---

MCMC回顾

平稳分布的定义

$$\pi =P\pi$$
称 $\pi$ 为马尔可夫链的平稳分布

直观上，如果马尔科夫链的平稳分布存在，那么以该平稳分布作为初始分布，而未来进行随机状态转移，之后任何一个时刻的状态分布都是该平稳分布。

细致平衡条件：
$$p_{ji}{\pi }_j=p_{ij}{\pi }_i,i,j=1,2,\cdots$$

MCMC中最重要的一个定理
满足细致平衡方程的状态分布$\pi$就是该马尔可夫链的平稳分布。即
$$\pi =P\pi$$

MH算法是怎样设计转移矩阵$P$的？

MH算法中，细致平衡条件有如下表述：

$\pi (x)$ 为要抽样的分布， 转移核为 $p(x^{\prime }|x)$
$$p(x^{\prime }|x)=q(x^{\prime }|x)\alpha (x,x^{\prime })$$

$q(x,x^{\prime })$ ,$\alpha (x,x^{\prime })$ 分别称为建议分布和接受分布

$$\alpha (x,x^{\prime })=\min \left\{1,\frac{\pi (x^{\prime })q(x|x^{\prime })}{\pi (x)q(x^{\prime }|x)}\right\}$$

且满足细致平衡条件：
$$\pi (x)p(x^{\prime }|x)=\pi (x^{\prime })p(x|x^{\prime })$$

---

我们可以自由的设计MH算法的转移核（transition kernel），对其进行的不同变化形成了不同的MCMC算法。从这个角度看，后来发展的HMC，slice sampling，Langevin Monte Carlo等等实际上都可以看作是MH的一种。很自然我们会思考有什么指导思想去自由设计转移核？

在实际应用MCMC过程中，可能会遇到下面的问题，比如：

1. MH的转移核（transition kernel）非常可能会造成随机游走（random walk behavior），使得状态在状态空间中进行小范围的探索而不能走出去，以至于退化成单纯的随机游走建议点，这样会造成有效采样个数（ESS）较小，抽样效率过低，还有就是当目标函数是多峰（multi-modal）的时候不能很好的找到所有的模式。

2. MH算法的拒绝率和转移核（transition kernel）有关，我们希望拒绝率尽量低。

设计转移核（transition kernel）的指导思想一般就是解决上面两个问题，而后提出的HMC算法就由于能很好的缓解上面的两个问题而出名。HMC的出现开辟了一条MCMC算法研究的新路。

---

Hamiltonian dynamics

拉格朗日方程

牛顿力学的运动微分方程:

$$m\frac{d^2\vec{r}}{d{t}^2}=\vec{F}$$

拉格朗日方程的特点是避开矢量力，而利用标量动能和势能来描述运动。

从牛顿方程出发推导拉格朗日方程

1、 单个质点在保守力场中运动: $U(\vec{r})$ —势能函数

直角坐标系中：$$\vec{r}=x\vec{i}+y\vec{j}+z\vec{k}$$
$$\vec{F}=-\nabla U=-\frac{\partial U}{\partial \vec{r}}=-\frac{\partial U}{\partial x}\vec{i}-\frac{\partial U}{\partial y}\vec{j}-\frac{\partial U}{\partial z}\vec{k}$$
由牛顿第二定律，质点的运动方程为: $$m\ddot{\vec{r}}=\vec{F}=-\frac{\partial U}{\partial \vec{r}}$$
分量形式：$$\{\begin{array}{l}m\ddot{x}=F_x\\ m\ddot{y}=F_y\\ m\ddot{z}=F_z\end{array}$$
又记$x,y,z$为$x_1,x_2,x_3$,上式又写为：
$$\{\begin{array}{c}m{\ddot{x}}_1={\vec{F}}_1=-\frac{\partial U}{\partial x_1}\\ m{\ddot{x}}_2={\vec{F}}_2=-\frac{\partial U}{\partial x_2}\\ m{\ddot{x}}_3={\vec{F}}_3=-\frac{\partial U}{\partial x_3}\end{array}$$
上式合写为：
$$m{\ddot{x}}_i=F_i=-\frac{\partial U}{\partial x_i}(i=1,2,3)$$

2、 直角坐标系中质点的动能为
$$T=\frac{1}{2}m\left({\dot{x}}^2+{\dot{y}}^2+{\dot{z}}^2\right)=\frac{1}{2}m\sum _{i=1}^3{\dot{x}}_i^2$$
动能对$\dot{x_i}$求偏导
$$\frac{\partial T}{\partial \dot{x_i}}=\frac{1}{2}m\times 2\dot{x_i}=m\dot{x_i}$$

上式再对时间求微分得:

$$\frac{d}{dt}\left(\frac{\partial T}{\partial \dot{x_i}}\right)=m\ddot{x_i}=F_i$$

由$\frac{d}{dt}\left(\frac{\partial T}{\partial \dot{x_i}}\right)=F_i$和$F_i=-\frac{\partial U}{\partial x_i}$二式相加得：
$$\frac{d}{dt}\left(\frac{\partial T}{\partial {\dot{x}}_i}\right)+\frac{\partial U}{\partial x_i}=0$$
3、引入拉格朗日函数$L$
$$L=T-U$$

动能$T$仅是速度$\dot{x_i}$的函数，势能$U$仅是坐标$x_i$的函数， 因此

$$\begin{array}{l}\frac{d}{dt}\left(\frac{\partial T}{\partial {\dot{x}}_i}\right)=\frac{d}{dt}\left(\frac{\partial (T-U)}{\partial {\dot{x}}_i}\right)=\frac{d}{dt}\left(\frac{\partial L}{\partial {\dot{x}}_i}\right)\end{array}$$
$$\begin{array}{l}\frac{\partial U}{\partial x_i}=\frac{\partial (U-T)}{\partial x_i}=-\frac{\partial L}{\partial x_i}\end{array}$$
以上两式相加得：
$$\frac{d}{dt}\left(\frac{\partial L}{\partial {\dot{x}}_i}\right)-\frac{\partial L}{\partial x_i}=0$$

此式即为用拉格朗日函数表示牛顿运动定律的拉格朗日方程。

可以证明，将$x_1,x_2,x_3$换成广义坐标$q_1,q_2,\dots ,q_s$，即可得到用广义坐标表示的具有$s$个自由度的系统的般形式的拉格朗日方程。
$$\frac{d}{dt}\left(\frac{\partial L}{\partial {\dot{q}}_i}\right)-\frac{\partial L}{\partial q_i}=0(i=1,2,\cdots ,$$