#前言
继续之前的Theano学习，本次主要学习混合蒙特卡洛(Hybrid Monte-Carlo Sampling)采样算法。
国际惯例，参考网址
Hybrid Monte-Carlo Sampling
Hybrid Monte Carlo
#理论
能量模型所使用的极大似然学习需要鲁棒的算法进行反向阶段的采样。比如玻尔兹曼机中负对数似然函数的梯度为：
−∂logp(x)∂θ=∂F(x)∂θ−∑x^p(x^)∂F(x^)∂θ-\frac{\partial \log p(x)}{\partial \theta}=\frac{\partial F(x)}{\partial \theta}-\sum_{\hat{x}}p(\hat{x})\frac{\partial F(\hat{x})}{\partial \theta} −∂θ∂logp(x)=∂θ∂F(x)−x^∑p(x^)∂θ∂F(x^)
上式包含两项，分别为正相和负相，不是说它们各项的符号，而是说它们对模型概率密度的作用，第一项通过降低对应的自由能来增大训练数据的概率密度，而第二项则是降低模型所生成的样本的概率密度，其实整体来说是让右边大点。
当我们训练RBM的时候，常采用了对比散度(CD)算法和持续性对比散度(PCD)算法，其实就是块吉布斯采样，条件分布p(h∣v)p(h|v)p(h∣v)和p(v∣h)p(v|h)p(v∣h)被用于执行马尔科夫链的转移(transition)操作。
在某些场景中，条件概率分布可能比较难采样(比如均值-协方差RBM需要很昂贵的矩阵逆操作)。而且，即使吉布斯采样能够被有效执行，它也不过是随机游走，这可能无法有效地满足某些分布。在此背景下，当从连续变量中采样时，混合蒙特卡洛采样(Hybrid Monte Carlo,HMC) 就是一种强有力的算法，主要通过哈密顿动力学(Hamiltonian dynamics)分治物理仿真系统避免了随机游走，这个过程具有预防比较奇怪的条件分布的潜力。
在HMC中，模型样本通过仿真物理学系统得到，粒子绕着高维场景运动，映射到隐藏动能中。粒子代表一个位置向量或者状态s∈RDs\in R^Ds∈RD，粒子的状态合成起来就是χ=(s,ϕ)\chi=(s,\phi)χ=(s,ϕ)，哈密顿就是势能(谢谢评论区@weixin_41669936指正)E(s)E(s)E(s)(这个和与能量模型的能量函数相同)与动能K(ϕ)K(\phi)K(ϕ)的和：
H(s,ϕ)=E(s)+K(ϕ)=E(s)+12∑iϕi2H(s,\phi)=E(s)+K(\phi)=E(s)+\frac{1}{2}\sum_i \phi_i^2 H(s,ϕ)=E(s)+K(ϕ)=E(s)+21i∑ϕi2
与直接从p(s)p(s)p(s)中采样不同的是，HMC从正则分布(canonical distribution)中采样，表达式为
p(s,ϕ)=1Zexp(−H(s,ϕ))=p(s)p(ϕ)p(s,\phi)=\frac{1}{Z}\exp(-H(s,\phi))=p(s)p(\phi) p(s,ϕ)=Z1exp(−H(s,ϕ))=p(s)p(ϕ)
由于两个变量相互独立，边缘化ϕ\phiϕ是平凡(可行)的，并且能够恢复原始的兴趣分布。
##哈密顿动力学
状态sss和速度ϕ\phiϕ 是改进过的，因此H(s,ϕ)H(s,\phi)H(s,ϕ)在整个仿真过程中保持为常量：
dsidt=∂H∂ϕi=ϕidϕidt=−∂H∂si=−∂E∂si\frac{d s_i}{d t}=\frac{\partial H}{\partial \phi_i}=\phi_i\\ \frac{d\phi_i}{d t}=-\frac{\partial H}{\partial s_i}=-\frac{\partial E}{\partial s_i} dtdsi=∂ϕi∂H=ϕidtdϕi=−∂si∂H=−∂si∂E
上述变换能够保持物理仿真系统的体积并且可逆，证明戳[这篇](Probabilistic Inference Using Markov Chain Monte Carlo Methods)文章。上述的动态可以被用于马尔科夫链的转移操作，并且保持p(s,ϕ)p(s,\phi)p(s,ϕ)不变。链本身不是各态遍历的，是因为动态仿真(或者动力学仿真)保持着一个固定的哈密顿函数H(s,ϕ)H(s,\phi)H(s,ϕ)，HMC交替执行哈密顿动力学步骤，使用吉布斯采样速度。由于p(s)p(s)p(s)和p(ϕ)p(\phi)p(ϕ)是独立的，因此ϕnew∼p(ϕ)\phi_{new}\sim p(\phi)ϕnew∼p(ϕ)是平凡的，因为p(ϕ∣s)=p(ϕ)p(\phi| s)=p(\phi)p(ϕ∣s)=p(ϕ)，其中p(ϕ)p(\phi)p(ϕ)经常采用单变量高斯分布。
##Leap-Frog算法
在实际中，由于时间离散(time discretization)问题，我们并无法精确模拟哈密顿动力学，有很多方法可以做到这一点。为了保证马尔科夫链的不变性，必须保持体积守恒和时间可逆性，而Leap-Frog算法就通过如下三步做到了这一点：
ϕi(t+ϵ2)=ϕi(t)−ϵ2∂E(s(t))sisi(t+ϵ)=si(t)+ϵϕi(t+ϵ2)ϕi(t+ϵ)=ϕi(t+ϵ2)−ϵ2∂E(s(t+ϵ))∂si\phi_i\left(t+\frac{\epsilon}{2}\right)=\phi_i(t)-\frac{\epsilon}{2}\frac{\partial E(s(t))}{s_i}\\ s_i(t+\epsilon)=s_i(t)+\epsilon\phi_i(t+\frac{\epsilon}{2})\\ \phi_i(t+\epsilon)=\phi_i\left(t+\frac{\epsilon }{2}\right)-\frac{\epsilon}{2}\frac{\partial E(s(t+\epsilon))}{\partial s_i} ϕi(t+2ϵ)=ϕi(t)−2ϵsi∂E(s(t))si(t+ϵ)=si(t)+ϵϕi(t+2ϵ)ϕi(t+ϵ)=ϕi(t+2ϵ)−2ϵ∂si∂E(s(t+ϵ))
因此可以执行半步速度更新，时刻是t+ϵ2t+\frac{\epsilon}{2}t+2ϵ，然后被用于计算s(t+ϵ)s(t+\epsilon)s(t+ϵ)和ϕ(t+ϵ)\phi(t+\epsilon)ϕ(t+ϵ)
##接受或者拒绝
在实际中，使用有限步长ϵ\epsilonϵ不会精确保存H(s,ϕ)H(s,\phi)H(s,ϕ)，在仿真中会引入偏置。而且由于浮点数的使用所导致的舍入误差表明变换将不会被完美保存下来。
HMC精确消除了这些影响，主要通过在经过n次leapfrog步骤后，增加梅特罗波利斯的接收/拒绝状态，新的状态χ′=(s′,ϕ′)\chi '=(s',\phi')χ′=(s′,ϕ′)被接受的概率是pacc(χ,χ′)p_{acc}(\chi,\chi')pacc(χ,χ′):
pacc(χ,χ′)=min(1,exp(−H(s′,ϕ′exp(−H(s,ϕ)))p_{acc}(\chi,\chi')=min\left(1,\frac{\exp(-H(s',\phi'}{\exp(-H(s,\phi))}\right) pacc(χ,χ′)=min(1,exp(−H(s,ϕ))exp(−H(s′,ϕ′)
HMC算法
三个步骤：
- 从单变量高斯分布中采样一个新的速度
- 执行n次
leapfrog，获取新的状态χ′\chi'χ′ - 判断新状态的引动是拒绝还是接受
动力学仿真
为了执行n次leapfrog，首先使用Scan定义可以迭代的一个函数，与直接实现Leap-Frog算法公式不同的是，我们需要注意的是通过执行ϕ\phiϕ的初始半步更新获得s(t+nϵ)s(t+n\epsilon)s(t+nϵ)与ϕ(t+nϵ)\phi(t+n\epsilon)ϕ(t+nϵ)，然后再对s,ϕs,\phis,ϕ进行n次完整的更新，以及对于ϕ\phiϕ再进行完整的和半步的更新，循环的形式如下：
- 先执行初始的半步更新
ϕi(t+ϵ2)=ϕi(t)−ϵ2∂E(s(t))∂sisi(t+ϵ)=si(t)+ϵϕi(t+ϵ2)\phi_i\left(t+\frac{\epsilon}{2}\right)=\phi_i(t)-\frac{\epsilon}{2}\frac{\partial E(s(t))}{\partial s_i}\\ s_i(t+\epsilon)=s_i(t)+\epsilon \phi_i\left(t+\frac{\epsilon}{2}\right) ϕi(t+2ϵ)=ϕi(t)−2ϵ∂si∂E(s(t))si(t+ϵ)=si(t)+ϵϕi(t+2ϵ)
-
再执行n次完整的更新
ϕi(t+(m−12)ϵ)=ϕi(t+(m−32ϵ))−ϵ∂E(s(t+(m−1)ϵ))∂sisi(t+mϵ)=si(t)+ϵϕi(t+(m−12)ϵ)\phi_i\left(t+\left(m-\frac{1}{2}\right)\epsilon\right)=\phi_i \left(t+\left(m-\frac{3}{2}\epsilon\right)\right)-\epsilon\frac{\partial E(s(t+(m-1)\epsilon))}{\partial s_i}\\ s_i(t+m\epsilon)=s_i(t)+\epsilon\phi_i\left(t+\left(m-\frac{1}{2}\right)\epsilon\right) ϕi(t+(m−21)ϵ)=ϕi(t+(m−23ϵ))−ϵ∂si∂E(s(t+(m−1)ϵ))si(t+mϵ)=si(t)+ϵϕi(t+(m−21)ϵ) -
最后一次更新ϕi\phi_iϕi
ϕi(t+nϵ)=ϕi(t+(n−12)ϵ)−ϵ2∂E(s(t+nϵ))∂si\phi_i(t+n\epsilon)=\phi_i\left(t+\left(n-\frac{1}{2}\right)\epsilon\right)-\frac{\epsilon}{2}\frac{\partial E(s(t+n\epsilon))}{\partial s_i} ϕi(t+nϵ)=ϕi(t+(n−21)ϵ)−2ϵ∂si∂E(s(t+nϵ))
#在Theano中实现HMC
在Theano选中，更新字典和共享变量为实现采样算法提供了较好的实现。采样器的当前状态可以使用Theano的共享变量表示，HMC更新可以使用Theano函数中的更新列表实现。
##搭建HMC
首先将HMC划分为以下子成分：
- 动力学仿真：符号python函数给出了初始位置和速度，可以执行n次
leapfrog更新，并可以为提议状态χ′\chi'χ′返回符号变量 - HMC移动：给定了初始位置的符号python函数，通过随机采样一个速度向量生成χ\chiχ，然后被称为动力学仿真，并且决定了χ→χ′\chi\to \chi'χ→χ′的转移是否被接受
- HMC更新：给定了HMC移动输出的python函数，为HMC的一次迭代生成更新列表
- HMC采样器：一个帮助类，将所有的东西集中起来
接下来我们就逐步开始了：
-
哈密顿动力学方程
H(s,ϕ)=E(s)+K(ϕ)=E(s)+12∑iϕi2H(s,\phi)=E(s)+K(\phi)=E(s)+\frac{1}{2}\sum_i \phi_i^2 H(s,ϕ)=E(s)+K(ϕ)=E(s)+21i∑ϕi2
注意，我们使用pospospos代表sss，使用velvelvel代表ϕ\phiϕ实现:def kinetic_energy(vel):return 0.5*(vel**2).sum(axis=1) def hamiltonian(pos,vel,energy_fn):return energy_fn(pos)+kinetic_energy(vel) -
动力学仿真
仿真阶段包含leapfrog阶段，用于循环采样，得到最终的pospospos和velvelveldef simulate_dynamics(initial_pos,initial_vel,stepsize,n_steps,energy_fn):def leapfrog(pos,vel,step):dE_dpos=TT.grad(energy_fn(pos).sum(),pos)new_vel=vel-step*dE_dpos;new_pos=pos+step*new_vel;return [new_pos,new_vel],{}initial_energy=energy_fn(initial_pos)dE_dpos=TT.grad(initial_energy.sum(),initial_pos)vel_half_step=initial_vel-0.5*stepsize*dE_dpos#半步pos_full_step=initial_pos+stepsize*vel_half_step#一步(all_pos,all_vel),scan_updates=theano.scan(leapfrog,outputs_info=[dict(initial=pos_full_step),dict(initial=vel_half_step),],non_sequences=[stepsize],n_steps=n_steps-1)final_pos=all_pos[-1]final_vel=all_vel[-1]assert not scan_updatesenergy=energy_fn(final_pos)final_vel=final_vel-0.5*stepsize*TT.grad(energy.sum(),final_pos)return final_pos,final_vel可以发现在动力学仿真系统
simulate_dynamics中，首先定义了leapfrog函数，执行次数为step，其实就是动力学仿真小结中的第一步，只不过把步数ϵ\epsilonϵ换成了step，然后分别写出经过step次迭代以后新的pospospos和velvelvel,然后进行第二步，也就是执行nnn次Leap-Frog算法，得到新的数据，这里需要注意的是，需要预先计算一下初始的sss和ϕ\phiϕ，而且ϕ\phiϕ是半步初始化，sss是一步初始化。 -
HMC移动
主要是计算当前状态到下一个状态的转移，包含metropolis_hastings的转移提议概率。给定了初始状态s∈RN×Ds\in R^{N\times D}s∈RN×D(位置)，能量函数E(s)E(s)E(s)(energy_fn)，HMC定义了一个符号图去计算n次HMC移动，步长是给定的stepsize
先写出metropolis_hastings转移接受概率def metropolis_hasting_accept(energy_prev,energy_next,s_rng):ediff=energy_prev-energy_next;return (TT.exp(ediff)-s_rng.uniform(size=energy_prev.shape))>=0然后定义HMC的移动
def hmc_move(s_rng,positions,energy_fn,stepsize,n_steps):initial_vel=s_rng.normal(size=positions.shape)final_pos,final_vel=simulate_dynamics(initial_pos=positions,initial_vel=initial_vel,stepsize=stepsize,n_steps=n_steps,energy_fn=energy_fn)accept=metropolis_hasting_accept(energy_prev=hamiltonian(positions,initial_vel,energy_fn),energy_next=hamiltonian(final_pos,final_vel,energy_fn),s_rng=s_rng)return accept,final_pos首先是从单变量高斯分布中采样一个随机速度，然后在动力学仿真系统中利用
metropolis_hastings转移概率执行n_steps次Leap-Frog算法，最终返回的是新状态final_pos被接受的符号变量accept，也就是说是否被接受。 -
HMC更新
无论HMC采样何时被调用，整个模型都必须有参数更新，包含position,stepsize,avg_acceptance_rate，这些参数被用于计算下一次的新状态。
先看看整个更新过程需要输入的参数：def hmc_updates(positions,stepsize,avg_acceptance_rate,final_pos,accept,target_acceptance_rate,stepsize_inc,stepsize_dec,stepsize_min,stepsize_max,avg_acceptance_slowness):首先依据接受概率计算新的位置
new_positionsaccept_matrix=accept.dimshuffle(0,*(('x',)*(final_pos.ndim-1))) new_positions=TT.switch(accept_matrix,final_pos,positions)在书写HMC参数更新的时候，首先计算HMC转移提议的平均接受率，使用具有时间长亮的指数移动平均：1−avgacceptanceslowness1-avg_acceptance_slowness1−avgacceptanceslowness，实现如下：
mean_dtype=theano.scalar.upcast(accept.dtype,avg_acceptance_rate.dtype) new_acceptance_rate=TT.add(avg_acceptance_slowness*avg_acceptance_rate,(1.0-avg_acceptance_slowness)*accept.mean(dtype=mean_dtype) )如果平均接受概率大于
target_accpetance_rate，那么就增加stepsize，增加量是stepsize_inc，反正按照stepsize_dec降低步长_new_stepsize=TT.switch(avg_acceptance_rate>target_acceptance_rate,stepsize*stepsize_inc,stepsize*stepsize_dec) new_stepsize=TT.clip(_new_stepsize,stepsize_min,stepsize_max)最终就可以返回更新参数了：
return [(positions,new_positions),(stepsize,new_stepsize),(avg_acceptance_rate,new_acceptance_rate)] -
定义采样器
HMC_sampler
主要成分是：new_from_shared_positions：构造函数，用于将hmc_move和hmc_updates所需要的不同的共享变量和字符放一起做分配，同时建立simulate函数，唯一的目标是执行hmc_updates所产生的更新。draw：调用simulate比较简单的方法，返回的是共享变量positions内容的副本。
class HMC_sampler(object):def __init__(self,**kwargs):self.__dict__.update(kwargs)@classmethoddef new_from_shared_positions(cls,shared_positions,energy_fn,initial_stepsize=0.01,target_acceptance_rate=.9,n_steps=20,stepsize_dec=0.98,stepsize_min=0.001,stepsize_max=0.25,stepsize_inc=1.02,avg_acceptance_slowness=0.9,seed=12345):stepsize=sharedX(initial_stepsize,'hmc_stepsize')avg_acceptance_rate=sharedX(target_acceptance_rate,'acg_acceptance_rate')s_rng=theano.sandbox.rng_mrg.MRG_RandomStreams(seed)accept,final_pos=hmc_move(s_rng,shared_positions,energy_fn,stepsize,n_steps)simulate_updates=hmc_updates(shared_positions,stepsize,avg_acceptance_rate,final_pos=final_pos,accept=accept,stepsize_min=stepsize_min,stepsize_max=stepsize_max,stepsize_inc=stepsize_inc,stepsize_dec=stepsize_dec,target_acceptance_rate=target_acceptance_rate,avg_acceptance_slowness=avg_acceptance_slowness)simulate=function([],[],updates=simulate_updates)return cls(positions=shared_positions,stepsize=stepsize,stepsize_min=stepsize_min,stepsize_max=stepsize_max,avg_acceptance_rate=avg_acceptance_rate,target_acceptance_rate=target_acceptance_rate,s_rng=s_rng,_updates=simulate_updates,simulate=simulate)def draw(self,**kwargs):self.simulate()return self.positions.get_value(borrow=False)
测试HMC
从一个多元高斯分布中做采样，先生成一个随机均值mu和协方差矩阵cov，这样可以用于定义对应的高斯分布的能量函数:gaussian_energy，然后通过分配一个共享变量position来初始化采样器的状态，最后将它和目标能量函数一起传递到HMC构造函数HMC_sampler
在生成大量的样本后，将实验均值和误差与真实值作比较：
def sampler_on_nd_gaussian(sampler_cls,burnin,n_samples,dim=10):batchsize=3rng=numpy.random.RandomState(123)mu=numpy.array(rng.rand(dim)*10,dtype=theano.config.floatX)cov=numpy.array(rng.rand(dim,dim),dtype=theano.config.floatX)cov=(cov+cov.T)/2.cov[numpy.arange(dim),numpy.arange(dim)]=1.0cov_inv=numpy.linalg.inv(cov)def gaussian_energy(x):return 0.5*(theano.tensor.dot((x-mu),cov_inv)*(x-mu)).sum(axis=1)position=rng.randn(batchsize,dim).astype(theano.config.floatX)position=theano.shared(position)sampler=sampler_cls(position,gaussian_energy,initial_stepsize=1e-3,stepsize_max=0.5)garbage=[sampler.draw() for r in range(burnin)]_samples=numpy.asarray([sampler.draw() for r in range(n_samples)])samples=_samples.T.reshape(dim,-1).Tprint('******TARGET VALUES*********')print('taget mean:',mu)print('taget conv:\n',cov)print('******EMPIRICAL MEAN/COV USING HMC********')print('empirical mean:',samples.mean(axis=0))print('empirical_cov:\n',numpy.cov(samples.T))print('******HMC INTERNALS*********')print('final stepsize',sampler.stepsize.get_value())print('final acceptance_rate',sampler.avg_acceptance_rate.get_value())return sampler
测试代码：
def test_hmc():sampler=sampler_on_nd_gaussian(HMC_sampler.new_from_shared_positions,burnin=1000,n_samples=1000,dim=5)assert abs(sampler.avg_acceptance_rate.get_value()-sampler.target_acceptance_rate)<.1assert sampler.stepsize.get_value()>=sampler.stepsize_minassert sampler.stepsize.get_value()<=sampler.stepsize_max
本人电脑的运行结果：
******TARGET VALUES*********
taget mean: [6.9646916 2.8613935 2.2685146 5.513148 7.1946898]
taget conv:[[1. 0.6619711 0.71141255 0.5576664 0.35753822][0.6619711 1. 0.310532 0.45455486 0.37991646][0.71141255 0.310532 1. 0.62800336 0.3800454 ][0.5576664 0.45455486 0.62800336 1. 0.5080787 ][0.35753822 0.37991646 0.3800454 0.5080787 1. ]]
******EMPIRICAL MEAN/COV USING HMC********
empirical mean: [6.947451 2.826787 2.2515843 5.518336 7.20614 ]
empirical_cov:[[1.00443031 0.67283693 0.73674632 0.59897923 0.35239996][0.67283693 0.97687277 0.31994575 0.45047961 0.31960069][0.73674632 0.31994575 1.03849032 0.71499541 0.43953283][0.59897923 0.45047961 0.71499541 1.05124182 0.55184436][0.35239996 0.31960069 0.43953283 0.55184436 0.99196057]]
******HMC INTERNALS*********
final stepsize 0.43232968
final acceptance_rate 0.91514903
结语
这个理论貌似有点小复杂啊，其实主要涉及到另一种循环采样方法，先构建了一个动力学系统，然后利用梅特罗波利斯-哈斯廷斯循环采样，整个过程还是蛮严谨的，参数太多了，光看理论，我自己估计都不一定能实现出来，虽然博客出来了，里面的细节还是得琢磨，另外贴一下个人的实验代码：链接：https://pan.baidu.com/s/1VUG41qQZHR3QWFPiKE_q0w 密码：saej
)
