當前位置：首頁 > 运维知识 > windows >内容正文

windows

【theano-windows】学习笔记十八——混合蒙特卡洛采样

發(fā)布時間：2023/12/13 windows 32 豆豆

生活随笔收集整理的這篇文章主要介紹了【theano-windows】学习笔记十八——混合蒙特卡洛采样小編覺得挺不錯的,現(xiàn)在分享給大家,幫大家做個參考.

#前言

繼續(xù)之前的Theano學習，本次主要學習混合蒙特卡洛(Hybrid Monte-Carlo Sampling)采樣算法。

國際慣例，參考網(wǎng)址

Hybrid Monte-Carlo Sampling

Hybrid Monte Carlo

#理論

能量模型所使用的極大似然學習需要魯棒的算法進行反向階段的采樣。比如玻爾茲曼機中負對數(shù)似然函數(shù)的梯度為：
$??log?p(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}$
上式包含兩項，分別為正相和負相，不是說它們各項的符號，而是說它們對模型概率密度的作用，第一項通過降低對應的自由能來增大訓練數(shù)據(jù)的概率密度，而第二項則是降低模型所生成的樣本的概率密度，其實整體來說是讓右邊大點。

當我們訓練RBM的時候，常采用了對比散度(CD)算法和持續(xù)性對比散度(PCD)算法，其實就是塊吉布斯采樣，條件分布 $p (h ∣ v)$ 和 $p (v ∣ h)$ 被用于執(zhí)行馬爾科夫鏈的轉移(transition)操作。

在某些場景中，條件概率分布可能比較難采樣(比如均值-協(xié)方差RBM需要很昂貴的矩陣逆操作)。而且，即使吉布斯采樣能夠被有效執(zhí)行，它也不過是隨機游走，這可能無法有效地滿足某些分布。在此背景下，當從連續(xù)變量中采樣時，混合蒙特卡洛采樣(Hybrid Monte Carlo,HMC) 就是一種強有力的算法，主要通過哈密頓動力學(Hamiltonian dynamics)分治物理仿真系統(tǒng)避免了隨機游走，這個過程具有預防比較奇怪的條件分布的潛力。

在HMC中，模型樣本通過仿真物理學系統(tǒng)得到，粒子繞著高維場景運動，映射到隱藏動能中。粒子代表一個位置向量或者狀態(tài) $s∈RDs\in R^D$ ，粒子的狀態(tài)合成起來就是 $χ=(s,?)\chi=(s,\phi)$ ，哈密頓就是勢能(謝謝評論區(qū)@weixin_41669936指正) $E (s)$ (這個和與能量模型的能量函數(shù)相同)與動能 $K(?)K(\phi)$ 的和：
$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$
與直接從 $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)$
由于兩個變量相互獨立，邊緣化 $?\phi$ 是平凡(可行)的，并且能夠恢復原始的興趣分布。

##哈密頓動力學

狀態(tài) $s$ 和速度 $?\phi$ 是改進過的，因此 $H(s,?)H(s,\phi)$ 在整個仿真過程中保持為常量：
$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}$
上述變換能夠保持物理仿真系統(tǒng)的體積并且可逆，證明戳[這篇](Probabilistic Inference Using Markov Chain Monte Carlo Methods)文章。上述的動態(tài)可以被用于馬爾科夫鏈的轉移操作，并且保持 $p(s,?)p(s,\phi)$ 不變。鏈本身不是各態(tài)遍歷的，是因為動態(tài)仿真(或者動力學仿真)保持著一個固定的哈密頓函數(shù) $H(s,?)H(s,\phi)$ ，HMC交替執(zhí)行哈密頓動力學步驟，使用吉布斯采樣速度。由于 $p (s)$ 和 $p(?)p(\phi)$ 是獨立的，因此 $?new～p(?)\phi_{new}\sim p(\phi)$ 是平凡的，因為 $p(?∣s)=p(?)p(\phi| s)=p(\phi)$ ，其中 $p(?)p(\phi)$ 經(jīng)常采用單變量高斯分布。

##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}$
因此可以執(zhí)行半步速度更新，時刻是 $t+?2t+\frac{\epsilon}{2}$ ，然后被用于計算 $s(t+?)s(t+\epsilon)$ 和 $?(t+?)\phi(t+\epsilon)$

##接受或者拒絕

在實際中，使用有限步長 $?\epsilon$ 不會精確保存 $H(s,?)H(s,\phi)$ ，在仿真中會引入偏置。而且由于浮點數(shù)的使用所導致的舍入誤差表明變換將不會被完美保存下來。

HMC精確消除了這些影響，主要通過在經(jīng)過n次leapfrog步驟后，增加梅特羅波利斯的接收/拒絕狀態(tài)，新的狀態(tài) $χ′=(s′,?′)\chi '=(s',\phi')$ 被接受的概率是 $pacc(χ,χ′)p_{acc}(\chi,\chi')$ :
$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)$

HMC算法

三個步驟：

從單變量高斯分布中采樣一個新的速度
執(zhí)行n次leapfrog，獲取新的狀態(tài) $χ′\chi'$
判斷新狀態(tài)的引動是拒絕還是接受

動力學仿真

為了執(zhí)行n次leapfrog，首先使用Scan定義可以迭代的一個函數(shù)，與直接實現(xiàn)Leap-Frog算法公式不同的是，我們需要注意的是通過執(zhí)行 $?\phi$ 的初始半步更新獲得 $s(t+n?)s(t+n\epsilon)$ 與 $?(t+n?)\phi(t+n\epsilon)$ ，然后再對 $s,?s,\phi$ 進行n次完整的更新，以及對于 $?\phi$ 再進行完整的和半步的更新，循環(huán)的形式如下：

先執(zhí)行初始的半步更新

$?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)$

再執(zhí)行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\phi_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}$

#在Theano中實現(xiàn)HMC

在Theano選中，更新字典和共享變量為實現(xiàn)采樣算法提供了較好的實現(xiàn)。采樣器的當前狀態(tài)可以使用Theano的共享變量表示，HMC更新可以使用Theano函數(shù)中的更新列表實現(xiàn)。

##搭建HMC

首先將HMC劃分為以下子成分：

動力學仿真：符號python函數(shù)給出了初始位置和速度，可以執(zhí)行n次leapfrog更新，并可以為提議狀態(tài) $χ′\chi'$ 返回符號變量
HMC移動：給定了初始位置的符號python函數(shù)，通過隨機采樣一個速度向量生成 $χ\chi$ ，然后被稱為動力學仿真，并且決定了 $χ→χ′\chi\to \chi'$ 的轉移是否被接受
HMC更新：給定了HMC移動輸出的python函數(shù)，為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$
注意，我們使用 $p o s$ 代表 $s$ ，使用 $v e l$ 代表 $?\phi$ 實現(xiàn):
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階段，用于循環(huán)采樣，得到最終的 $p o s$ 和 $v e l$
def 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
可以發(fā)現(xiàn)在動力學仿真系統(tǒng)simulate_dynamics中，首先定義了leapfrog函數(shù)，執(zhí)行次數(shù)為step，其實就是動力學仿真小結中的第一步，只不過把步數(shù) $?\epsilon$ 換成了step，然后分別寫出經(jīng)過step次迭代以后新的 $p o s$ 和 $v e l$ ,然后進行第二步，也就是執(zhí)行 $n$ 次Leap-Frog算法，得到新的數(shù)據(jù)，這里需要注意的是，需要預先計算一下初始的 $s$ 和 $?\phi$ ，而且 $?\phi$ 是半步初始化， $s$ 是一步初始化。
HMC移動
主要是計算當前狀態(tài)到下一個狀態(tài)的轉移，包含metropolis_hastings的轉移提議概率。給定了初始狀態(tài) $s∈RN×Ds\in R^{N\times D}$ (位置)，能量函數(shù) $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
首先是從單變量高斯分布中采樣一個隨機速度，然后在動力學仿真系統(tǒng)中利用metropolis_hastings轉移概率執(zhí)行n_steps次Leap-Frog算法，最終返回的是新狀態(tài)final_pos被接受的符號變量accept，也就是說是否被接受。
HMC更新
無論HMC采樣何時被調用，整個模型都必須有參數(shù)更新，包含position,stepsize,avg_acceptance_rate，這些參數(shù)被用于計算下一次的新狀態(tài)。
先看看整個更新過程需要輸入的參數(shù)：
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):
首先依據(jù)接受概率計算新的位置new_positions
accept_matrix=accept.dimshuffle(0,*(('x',)*(final_pos.ndim-1))) new_positions=TT.switch(accept_matrix,final_pos,positions)
在書寫HMC參數(shù)更新的時候，首先計算HMC轉移提議的平均接受率，使用具有時間長亮的指數(shù)移動平均： $1-avg_acceptance_slowness$ ，實現(xiàn)如下：
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)
最終就可以返回更新參數(shù)了：
return [(positions,new_positions),(stepsize,new_stepsize),(avg_acceptance_rate,new_acceptance_rate)]
定義采樣器HMC_sampler
主要成分是：
- new_from_shared_positions：構造函數(shù)，用于將hmc_move和hmc_updates所需要的不同的共享變量和字符放一起做分配，同時建立simulate函數(shù)，唯一的目標是執(zhí)行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和協(xié)方差矩陣cov，這樣可以用于定義對應的高斯分布的能量函數(shù):gaussian_energy，然后通過分配一個共享變量position來初始化采樣器的狀態(tài)，最后將它和目標能量函數(shù)一起傳遞到HMC構造函數(shù)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

結語

這個理論貌似有點小復雜啊，其實主要涉及到另一種循環(huán)采樣方法，先構建了一個動力學系統(tǒng)，然后利用梅特羅波利斯-哈斯廷斯循環(huán)采樣，整個過程還是蠻嚴謹?shù)?#xff0c;參數(shù)太多了，光看理論，我自己估計都不一定能實現(xiàn)出來，雖然博客出來了，里面的細節(jié)還是得琢磨，另外貼一下個人的實驗代碼：鏈接：https://pan.baidu.com/s/1VUG41qQZHR3QWFPiKE_q0w 密碼：saej

總結

以上是生活随笔為你收集整理的【theano-windows】学习笔记十八——混合蒙特卡洛采样的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內容還不錯，歡迎將生活随笔推薦給好友。

上一篇： echart 自定义 formatter
下一篇：【theano-windows】学习笔记