當(dāng)前位置：首頁(yè) > 运维知识 > windows >内容正文

windows

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

發(fā)布時(shí)間：2023/12/13 windows 27 豆豆

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

#前言

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

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

Hybrid Monte-Carlo Sampling

Hybrid Monte Carlo

#理論

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

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

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

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

##哈密頓動(dòng)力學(xué)

狀態(tài) $s$ 和速度 $?\phi$ 是改進(jìn)過的，因此 $H(s,?)H(s,\phi)$ 在整個(gè)仿真過程中保持為常量：
$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)文章。上述的動(dòng)態(tài)可以被用于馬爾科夫鏈的轉(zhuǎn)移操作，并且保持 $p(s,?)p(s,\phi)$ 不變。鏈本身不是各態(tài)遍歷的，是因?yàn)閯?dòng)態(tài)仿真(或者動(dòng)力學(xué)仿真)保持著一個(gè)固定的哈密頓函數(shù) $H(s,?)H(s,\phi)$ ，HMC交替執(zhí)行哈密頓動(dòng)力學(xué)步驟，使用吉布斯采樣速度。由于 $p (s)$ 和 $p(?)p(\phi)$ 是獨(dú)立的，因此 $?new～p(?)\phi_{new}\sim p(\phi)$ 是平凡的，因?yàn)?span id="ze8trgl8bvbq" class="katex--inline"> $p(?∣s)=p(?)p(\phi| s)=p(\phi)$ ，其中 $p(?)p(\phi)$ 經(jīng)常采用單變量高斯分布。

##Leap-Frog算法

在實(shí)際中，由于時(shí)間離散(time discretization)問題，我們并無法精確模擬哈密頓動(dòng)力學(xué)，有很多方法可以做到這一點(diǎn)。為了保證馬爾科夫鏈的不變性，必須保持體積守恒和時(shí)間可逆性，而Leap-Frog算法就通過如下三步做到了這一點(diǎn)：
$?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í)行半步速度更新，時(shí)刻是 $t+?2t+\frac{\epsilon}{2}$ ，然后被用于計(jì)算 $s(t+?)s(t+\epsilon)$ 和 $?(t+?)\phi(t+\epsilon)$

##接受或者拒絕

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

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算法

三個(gè)步驟：

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

動(dòng)力學(xué)仿真

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

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

##搭建HMC

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

動(dòng)力學(xué)仿真：符號(hào)python函數(shù)給出了初始位置和速度，可以執(zhí)行n次leapfrog更新，并可以為提議狀態(tài) $χ′\chi'$ 返回符號(hào)變量
HMC移動(dòng)：給定了初始位置的符號(hào)python函數(shù)，通過隨機(jī)采樣一個(gè)速度向量生成 $χ\chi$ ，然后被稱為動(dòng)力學(xué)仿真，并且決定了 $χ→χ′\chi\to \chi'$ 的轉(zhuǎn)移是否被接受
HMC更新：給定了HMC移動(dòng)輸出的python函數(shù)，為HMC的一次迭代生成更新列表
HMC采樣器：一個(gè)幫助類，將所有的東西集中起來

接下來我們就逐步開始了：

哈密頓動(dòng)力學(xué)方程
$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$ 實(shí)現(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)
動(dòng)力學(xué)仿真
仿真階段包含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)在動(dòng)力學(xué)仿真系統(tǒng)simulate_dynamics中，首先定義了leapfrog函數(shù)，執(zhí)行次數(shù)為step，其實(shí)就是動(dòng)力學(xué)仿真小結(jié)中的第一步，只不過把步數(shù) $?\epsilon$ 換成了step，然后分別寫出經(jīng)過step次迭代以后新的 $p o s$ 和 $v e l$ ,然后進(jìn)行第二步，也就是執(zhí)行 $n$ 次Leap-Frog算法，得到新的數(shù)據(jù)，這里需要注意的是，需要預(yù)先計(jì)算一下初始的 $s$ 和 $?\phi$ ，而且 $?\phi$ 是半步初始化， $s$ 是一步初始化。
HMC移動(dòng)
主要是計(jì)算當(dāng)前狀態(tài)到下一個(gè)狀態(tài)的轉(zhuǎn)移，包含metropolis_hastings的轉(zhuǎn)移提議概率。給定了初始狀態(tài) $s∈RN×Ds\in R^{N\times D}$ (位置)，能量函數(shù) $E (s)$ (energy_fn)，HMC定義了一個(gè)符號(hào)圖去計(jì)算n次HMC移動(dòng)，步長(zhǎng)是給定的stepsize
先寫出metropolis_hastings轉(zhuǎn)移接受概率
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的移動(dòng)
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
首先是從單變量高斯分布中采樣一個(gè)隨機(jī)速度，然后在動(dòng)力學(xué)仿真系統(tǒng)中利用metropolis_hastings轉(zhuǎn)移概率執(zhí)行n_steps次Leap-Frog算法，最終返回的是新狀態(tài)final_pos被接受的符號(hào)變量accept，也就是說是否被接受。
HMC更新
無論HMC采樣何時(shí)被調(diào)用，整個(gè)模型都必須有參數(shù)更新，包含position,stepsize,avg_acceptance_rate，這些參數(shù)被用于計(jì)算下一次的新狀態(tài)。
先看看整個(gè)更新過程需要輸入的參數(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ù)接受概率計(jì)算新的位置new_positions
accept_matrix=accept.dimshuffle(0,*(('x',)*(final_pos.ndim-1))) new_positions=TT.switch(accept_matrix,final_pos,positions)
在書寫HMC參數(shù)更新的時(shí)候，首先計(jì)算HMC轉(zhuǎn)移提議的平均接受率，使用具有時(shí)間長(zhǎng)亮的指數(shù)移動(dòng)平均： $1-avg_acceptance_slowness$ ，實(shí)現(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降低步長(zhǎng)
_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：構(gòu)造函數(shù)，用于將hmc_move和hmc_updates所需要的不同的共享變量和字符放一起做分配，同時(shí)建立simulate函數(shù)，唯一的目標(biāo)是執(zhí)行hmc_updates所產(chǎn)生的更新。
- draw：調(diào)用simulate比較簡(jiǎn)單的方法，返回的是共享變量positions內(nèi)容的副本。
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)

測(cè)試HMC

從一個(gè)多元高斯分布中做采樣，先生成一個(gè)隨機(jī)均值mu和協(xié)方差矩陣cov，這樣可以用于定義對(duì)應(yīng)的高斯分布的能量函數(shù):gaussian_energy，然后通過分配一個(gè)共享變量position來初始化采樣器的狀態(tài)，最后將它和目標(biāo)能量函數(shù)一起傳遞到HMC構(gòu)造函數(shù)HMC_sampler

在生成大量的樣本后，將實(shí)驗(yàn)均值和誤差與真實(shí)值作比較：

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

測(cè)試代碼：

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

本人電腦的運(yùn)行結(jié)果：

******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

結(jié)語

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

總結(jié)

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

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

上一篇： echart 自定义 formatter
下一篇： [数据集]新浪微博数据集Microblo