HMC

课上没听懂，看了一些博客再结合老师的课件之后对一些符号更加迷糊了，于是本强迫症决定重新整理以下整个理解的思路。肯定有错误或者不足之处。欢迎指出，后续改正。

Hamiltonian Monte Carlo (HMC) 算法，这是一种基于哈密顿力学（Hamiltonian mechanics）的抽样算法，它的效率比一般的MH算法要更高。通俗说就是更快达到逼近真实分布的采样。 看图说话：(盗图) 黑点是从真实分布获得的采样值；红圈是模拟采样，折线应该连接的上一步与下一步模拟的采样值。 HMC：几步后就基本接近真实的分布啦~ MH：要走十多步才接近真实分布的采样值，并且后续一直不太稳定（采样值老是跑到边边角角等出现概率较小的地方）。

参考知乎 哈密顿动力系统是一种由哈密顿方程驱动的动力系统，由哈密顿于1833年建立。 哈密尔顿动力学利用给定物体在时刻 t t t的位置 x x x和速度 v v v来描述其运动过程。位置和速度分别对应着势能和动能，故势能和动能分别为关于位置和速度的函数，记势能为 U ( x ) U(x) U(x)，动能为 K ( v ) K(v) K(v).

注1：物理学中，动能公式为: 1 2 m v 2 \frac{1}{2}mv^2 21​mv2（m代表质量）.速度与位置的关系: v = d x d t v=\frac{dx}{dt} v=dtdx​

系统的总能量为势能与动能之和,为了后续求出的结果更漂亮，把动能由 K ( v ) K(v) K(v)改写为 K ( p ) K(p) K(p)， p = m v p=mv p=mv代表动量。系统总能量记为: H ( x , p ) = U ( x ) + K ( p ) (1) H(x,p)=U(x)+K(p) \tag{1} H(x,p)=U(x)+K(p)(1) 根据能力守恒定理，没有外力干扰下，总能量（关于时间 t t t）不变即恒为常量。进一步求偏导可得到哈密尔顿动力学的偏微分方程。 0 = ∂ H ( x , p ) ∂ t = ∂ H ∂ x d x d t + ∂ H ∂ p d p d t (2) 0=\frac{\partial H(x,p)}{\partial t}=\frac{\partial H}{\partial x}\frac{d x}{dt}+\frac{\partial H}{\partial p}\frac{d p}{d t} \tag{2} 0=∂t∂H(x,p)​=∂x∂H​dtdx​+∂p∂H​dtdp​(2) ∂ H ( x , p ) ∂ p = ∂ K ( p ) ∂ p = d x d t (3) \frac{\partial H(x,p)}{\partial p}=\frac{\partial K(p)}{\partial p}=\frac{dx}{dt} \tag{3} ∂p∂H(x,p)​=∂p∂K(p)​=dtdx​(3) 由(2),(3)可得哈密顿方程： d p d t = − ∂ H ∂ x = − ∂ U ( x ) ∂ x , d x d t = ∂ K ( p ) ∂ p (4) \frac{d p}{d t}=-\frac{\partial H}{\partial x}=-\frac{\partial U(x)}{\partial x}, \quad \frac{d x}{d t}=\frac{\partial K(p)}{\partial p} \tag{4} dtdp​=−∂x∂H​=−∂x∂U(x)​,dtdx​=∂p∂K(p)​(4)

注2：由注1可得 K ( p ) = 1 2 m v 2 = 1 2 m p 2 , K ′ ( p ) = 1 m p = v K(p)=\frac{1}{2}mv^2=\frac{1}{2m}p^2,K'(p)=\frac{1}{m}p=v K(p)=21​mv2=2m1​p2,K′(p)=m1​p=v.

哈密顿方程描述了动能和势能随着时间相互转换的过程。

首先将后验密度写为如下形式： π ( θ ∣ X ) ∝ exp ⁡ ( − U ( θ ) ) \pi (\theta|X) \propto \exp(-U(\theta)) π(θ∣X)∝exp(−U(θ)). U \quad U U即势能函数， θ \theta θ即为位置。 接着将动能定义为： K ( p ) = 1 2 p T M − 1 p K(p)=\frac{1}{2}p^TM^{-1}p K(p)=21​pTM−1p，则总能量 H ( θ , p ) = U ( θ ) + K ( p ) = − l n ( π ( θ ∣ X ) ) + 1 2 p T M − 1 p H(\theta,p)=U(\theta)+K(p)=-ln(\pi(\theta|X))+\frac{1}{2}p^TM^{-1}p H(θ,p)=U(θ)+K(p)=−ln(π(θ∣X))+21​pTM−1p

注3：这样定义实则是将势能函数定义为后验密度的负对数似然函数，而动能函数则是正态分布密度函数的负对数形式（即把动量 p p p看作服从正态分布的随机变量，至于为什么是正态，应该只是一种简单且符合系统物理意义的选择），因此总能量即为后验密度与动量 p p p服从的正态密度的乘积再取负对数。希望读者能从这些解释中更好理解上面的公式。 故结合上述解释以及公式，可得 θ , p \theta, p θ,p的联合密度为： π ( θ , p ∣ X ) = π ( θ ∣ X ) ϕ ( p ; M ) ∝ exp ⁡ ( − H ( θ , p ) ) \pi(\theta,p|X)=\pi(\theta|X)\phi(p;M) \propto \exp(-H(\theta,p)) π(θ,p∣X)=π(θ∣X)ϕ(p;M)∝exp(−H(θ,p)). 由此也可见 θ , p \theta, p θ,p的独立性。独立性的用处后文会说明。

根据(4)式，可得： d θ = ∂ K ( p ) ∂ p d t = M − 1 p d t , d p = − ∂ U ( θ ) ∂ θ d t (5) d\theta = \frac{\partial K(p)}{\partial p}dt=M^{-1}pdt, \quad dp=-\frac{\partial U(\theta)}{\partial \theta}dt \tag{5} dθ=∂p∂K(p)​dt=M−1pdt,dp=−∂θ∂U(θ)​dt(5) 此即联系上了哈密顿动力系统中的动能、势能转化过程。

类比构造马尔科夫链使得平稳分布为目标分布从而使用转移过程中的状态作为采样模拟值的想法，HMC是构造哈密顿动力系统使得势能为目标分布（的负对数形式），总能量为联合密度（的负对数形式），再根据该系统能量转化过程（哈密顿方程）来转移目标状态（即参数更新）。 更近一步，HMC的更新过程中还可以使用类似MH采样的接受拒绝采样方法（以联合密度作为马尔科夫链的平稳分布）。接受率公式如下： A = m i n { 1 , π ( θ ∗ , p ∗ ∣ X ) π ( θ , p ∣ X ) } (6) A = min\{1, \frac{\pi(\theta^*,p^*|X)}{\pi(\theta,p|X)} \} \tag{6} A=min{1,π(θ,p∣X)π(θ∗,p∗∣X)​}(6) 其中 θ ∗ , p ∗ \theta^*,p^* θ∗,p∗代表继 θ , p \theta,p θ,p之后下一步的参数更新值。

注4：HMC实施的是对联合分布的采样，其参数更新通常指的是目标参数 θ \theta θ与动量参数 p同时更新（采用(5)式更新），具体更新方法就在后文。 然而我们并不关心“虚构 ”出来的动量参数 p p p，由注3中讲到的独立性，我们进行联合采样后可以直接扔掉 p p p，保留我们想要的参数 θ \theta θ即可。

在(5)式中，如果把 d t dt dt看作步长，那么 d θ , d p d\theta, dp dθ,dp即为对应参数的更新步长。那么一个常见的想法是： θ ( t + ϵ ) = θ ( t ) + ϵ d θ = θ ( t ) + M − 1 p ( t ) (7) \theta(t+\epsilon)=\theta(t)+\epsilon d\theta=\theta(t)+M^{-1}p(t) \tag{7} θ(t+ϵ)=θ(t)+ϵdθ=θ(t)+M−1p(t)(7) p ( t + ϵ ) = p ( t ) + ϵ d p = p ( t ) − ∇ U ( θ ( t ) ) (8) \quad p(t+\epsilon)=p(t)+\epsilon dp=p(t)-\nabla U(\theta(t)) \tag{8} p(t+ϵ)=p(t)+ϵdp=p(t)−∇U(θ(t))(8) 此即利用欧拉折线法将连续时间的哈密顿系统离散化。 但这没有充分利用信息，也可以将其改进为： p ( t + ϵ ) = p ( t ) − ϵ ∇ U ( θ ( t ) ) (9) \quad p(t+\epsilon)=p(t)-\epsilon \nabla U(\theta(t)) \tag{9} p(t+ϵ)=p(t)−ϵ∇U(θ(t))(9) θ ( t + ϵ ) = θ ( t ) + ϵ M − 1 p ( t + ϵ ) (10) \theta(t+\epsilon)=\theta(t)+\epsilon M^{-1}p(t+\epsilon) \tag{10} θ(t+ϵ)=θ(t)+ϵM−1p(t+ϵ)(10)但是这还不够，有名的Leapfrog方法是这样的： p ( t + ϵ / 2 ) = p ( t ) − ( ϵ / 2 ) ∇ U ( θ ( t ) ) (11) \quad p(t+\epsilon/2)=p(t)-(\epsilon/2)\nabla U(\theta(t)) \tag{11} p(t+ϵ/2)=p(t)−(ϵ/2)∇U(θ(t))(11) θ ( t + ϵ ) = θ ( t ) + M − 1 p ( t + ϵ / 2 ) (12) \theta(t+\epsilon)=\theta(t)+M^{-1}p(t+\epsilon/2) \tag{12} θ(t+ϵ)=θ(t)+M−1p(t+ϵ/2)(12) p ( t + ϵ ) = p ( t + ϵ / 2 ) − ( ϵ / 2 ) ∇ U ( θ ( t + ϵ ) ) (13) p(t+\epsilon)=p(t+\epsilon/2)-(\epsilon/2)\nabla U(\theta(t+\epsilon)) \tag{13} p(t+ϵ)=p(t+ϵ/2)−(ϵ/2)∇U(θ(t+ϵ))(13) Leapfrog方法的大概含义是先走半步更新 p p p，再走一步更新 θ \theta θ，最后走半步更新 p p p. 每一步都利用到了最新更新的信息。

详细的采样步骤参看此处.

优缺点：HMC遍历参数空间的速度一般比MH要快，但也存在一些缺点，例如有更多的参数要调、单步计算复杂度比MH高等。

全称Stochastic Gradient Hamiltonian Monte Carlo。 想法是引入二次抽样的思想来提高算法效率。我的理解是后验密度中的X本身就是抽样值，如果要利用全部样本X计算更新参数需要的梯度 ∇ U \nabla U ∇U会很耗时耗力（特别是样本量较大的时候）。这时候就需要从已得到的样本X中进行二次抽样计算梯度。计算公式如下：（即按照子样本的容量比例放大梯度） ∇ U ~ ( θ ) = − ∣ D ∣ ∣ D ~ ∣ ∑ X ∈ D ~ ∇ ln ⁡ π ( θ ∣ X ) , D ~ ⊂ D \nabla \tilde{U}(\theta)= -\frac{|D|}{|\tilde{D}|}\sum_{X\in \tilde{D}}\nabla \ln\pi(\theta|X), \quad\tilde{D}\subset D ∇U~(θ)=−∣D~∣∣D∣​X∈D~∑​∇lnπ(θ∣X),D~⊂D π ( θ ∣ X ) ∝ π ( X ∣ θ ) π ( θ ) \pi(\theta|X) \propto \pi(X|\theta)\pi(\theta) π(θ∣X)∝π(X∣θ)π(θ), 故 ∇ ln ⁡ π ( θ ∣ X ) = ∇ ln ⁡ π ( X ∣ θ ) + ∇ ln ⁡ π ( θ ) \nabla \ln\pi(\theta|X)=\nabla \ln\pi(X|\theta)+\nabla \ln\pi(\theta) ∇lnπ(θ∣X)=∇lnπ(X∣θ)+∇lnπ(θ).

此时，由于二次抽样的随机性，在引入随机梯度之后，整个系统发生了变化，导致参数更新有了新的随机影响因素。相当于在光滑地面上运动的小球受到随机方向的风的影响，导致了我们的目标 π ( θ , p ∣ X ) \pi(\theta,p|X) π(θ,p∣X)不再是该系统的平稳分布。（这里基本是老师上课讲的内容，我只能大概理解） 系统变为： d θ = M − 1 p d t , d p = − ∇ U ( θ ( t ) ) d t + N ( 0 , 2 B d t ) d\theta =M^{-1}pdt, \quad dp=-\nabla U(\theta(t))dt+N(0,2Bdt) dθ=M−1pdt,dp=−∇U(θ(t))dt+N(0,2Bdt)

具体证明可以参考Chen, Tianqi, Emily Fox, and Carlos Guestrin. “Stochastic gradient hamiltonian monte carlo.” International conference on machine learning. PMLR, 2014.

因此SGHMC利用“摩擦力”来面对疾风。引入摩擦力后的系统变为： d θ = M − 1 p d t , d p = − ∇ U ( θ ( t ) ) d t + N ( 0 , 2 B d t ) − B M − 1 p d t d\theta =M^{-1}pdt, \quad dp=-\nabla U(\theta(t))dt+N(0,2Bdt)-BM^{-1}pdt dθ=M−1pdt,dp=−∇U(θ(t))dt+N(0,2Bdt)−BM−1pdt 摩擦力降低了系统的能量，从而使“风”的影响相对变小。

SGHMC有许多参数可调，且在实际应用中有很多近似计算方法。具体案例可以参见上述参考文献。

在标准HMC和SGHMC中，虽然 π ( θ , p ∣ X ) \pi(\theta,p|X) π(θ,p∣X)是平稳分布，但是由于连续过程离散化，导致模型在实际应用中存在误差，除非随着迭代次数增加 ϵ → 0 \epsilon \rightarrow 0 ϵ→0，否则理论上都应该在算法的最后加上MH步骤来进行修正。 但是在实际中，递减的 ϵ \epsilon ϵ 会影响马尔可夫链的移动速度，因此可以以一定的误差为代价，使用常数 ϵ \epsilon ϵ 且不加MH修正。在标准HMC中，即使不加MH修正，对计算效率的提升也比较有限。

注：这里的意思应该是对于标准HMC（使用固定的 ϵ \epsilon ϵ），即使不加 MH步骤，因其使用所有样本计算梯度故其每次迭代都很慢。

上图！各种方法的比较如下 ： 此处盗图自老师课件 说明：“Standard”的意思应该是使用固定的 ϵ \epsilon ϵ ，更新参数的方法使用Leapfrog；“with MH”,"no MH"分别意味着加与不加MH步骤；"Naive"意味着SGHMC中不加入“摩擦力”。

在图中对应的案例中，不加“摩擦力”的SGHMC只要加入MH步骤即可有良好的效果。只有不加“摩擦力”以及MH步骤的方法效果很差，其余方法效果相近。

尝试利用HMC和SGHMC对正态分布参数 ( μ , σ 2 ) (\mu,\sigma^2) (μ,σ2)进行抽样。

有空更吧~