當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

MCMC算法大统一: Involutive MCMC

發布時間：2023/12/9 编程问答 32 豆豆

生活随笔收集整理的這篇文章主要介紹了 MCMC算法大统一: Involutive MCMC 小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

標準采樣方法

先從最基本的方法說起，可能很多人都知道只要可以對分布函數 $F(x)\displaystyle F( x)$ 求逆，并從均勻分布中采樣u,并將u代進逆函數中就能得到x的樣本，即 $x=F?1(u),u～U(0,1)\displaystyle x=F^{-1}( u) ,u\sim U( 0,1)$ 。他的原理是什么？其實他的出發點是找到一個從均勻分布到目標分布的可逆變換 $g\displaystyle g$ ：

$u)\\ p( x) =p( u)\left| \frac{du}{dx}\right| =p_{u}\left( g^{-1}( x)\right)\left| \frac{dg^{-1}( x)}{dx}\right|$

因為我們要保證x是我們的目標分布，那么他一定要滿足 $∫?∞xp(x)dx=F(x)\displaystyle \int ^{x}_{-\infty } p( x) dx=F( x)$ ，因為 $pu(f?1(x))\displaystyle p_{u}\left( f^{-1}( x)\right)$ 是0到1的均勻分布所以他等于1，于是必須有 $∫?∞x∣dg?1(x)dx∣dx=F(x)?g?1(x)=F(x)?g(x)=F?1(x)\displaystyle \int ^{x}_{-\infty }\left| \frac{dg^{-1}( x)}{dx}\right| dx=F( x) \Longrightarrow g^{-1}( x) =F( x) \Longrightarrow g( x) =F^{-1}( x)$ 。
然而該方法的問題在于，并不是每個概率分布的都可以寫出來并且可以求逆的。對于那些特別復雜的分布這樣方法就無能為力了。

拒絕采樣

現在我們的目標是要采樣 $p(x)\displaystyle p( x)$ 的分布。為了采樣這個分布，實際上，我們只需要知道 $p~(x)\displaystyle \tilde{p}( x)$ 就足夠了， $p~(x)\displaystyle \tilde{p}( x)$ 是沒有標準化的 $p(x)=p~(x)/Z\displaystyle p( x) =\tilde{p}( x) /Z$ , 其中 $Z=∫p~(x)dx\displaystyle Z=\int \tilde{p}( x) dx$ 是不需要知道的。那么拒絕采樣的流程是這樣的，首先隨便找一個proposal distribution $q(x)\displaystyle q( x)$ 使得 $Mq(x)?p~(x)\displaystyle Mq( x) \geqslant \tilde{p}( x)$ ，然后從q中隨機一個樣本 $x～q(x)\displaystyle x\sim q( x)$ ，再從均勻分布中隨機一個樣本 $u～U(0,1)\displaystyle u\sim U( 0,1)$ ，如果 $u?p~(x)Mq(x)\displaystyle u\geqslant \frac{\tilde{p}( x)}{Mq( x)}$ 則拒絕，否則接受。為什么這樣就能從p中采樣？

訂正：圖中的 $p\displaystyle p$ 應該是 $p~\displaystyle \tilde{p}$

首先，因為u是0,1區間，所以一定有 $0?uMq(x)?Mq(x)\displaystyle 0\leqslant uMq( x) \leqslant Mq( x)$ ，那么u的作用其實就是選擇0到 $Mq(x)\displaystyle Mq( x)$ 之間的值，如果 $p~(x)?uMq(x)\displaystyle \tilde{p}( x) \leqslant uMq( x)$ ，那么就拒絕掉樣本，也就是上圖白色部分，否則 $p~(x)?uMq(x)\displaystyle \tilde{p}( x) \geqslant uMq( x)$ 就是落在黑色區域，則接受。顯然，對于x軸上的任意一個點 $x′\displaystyle x'$ ，其對應接受的概率就是 $p~(x′)Mq(x′)\displaystyle \frac{\tilde{p}( x')}{Mq( x')}$ 。不過每個點，因為q，抽到的概率都不一樣，所以我們要把q抽中的概率也算上去。

那么用這個方法我們可以得到一個序列的樣本 $(x1,accept1,...,xi,accepti)\displaystyle ( x_{1} ,accept_{1} ,...,x_{i} ,accept_{i})$ ，其中 $xi～q(x)\displaystyle x_{i} \sim q( x)$ 是由 $q\displaystyle q$ 產生的，而且每個pair $(xi,accepti)\displaystyle ( x_{i} ,accept_{i})$ 都是相互獨立的，那么我們可以證明， $p(x∣accept)\displaystyle p( x|accept)$ 就等于我們想要的分布 $p(x)\displaystyle p( x)$ :

$p(x∣accept)=p(accept∣x)q(x)p(accpet)=p~(x)Mq(x)q(x)∫p~(x)Mq(x)q(x)dx=p~(x)∫p~(x)dx=p(x)/Z∫p(x)/Zdx=p(x)\begin{aligned} p( x|accept) & =\frac{p( accept|x) q( x)}{p( accpet)}\\ & =\frac{\frac{\tilde{p}( x)}{Mq( x)} q( x)}{\int \frac{\tilde{p}( x)}{Mq( x)} q( x) dx}\\ & =\frac{\tilde{p}( x)}{\int \tilde{p}( x) dx}\\ & =\frac{p( x) /Z}{\int p( x) /Zdx}\\ & =p( x) \end{aligned}$

通過拒絕采樣，我們可以知道，對于任意的分布，我們可以設計一種接受-拒絕的概率，從而使得所有接受的點都是該分布的樣本。然而他的問題在于，他對propose distribution q的要求很高，試想下，如果q與真實p的差距過大，那么幾乎每個樣本都會被拒絕掉，效率很低。那如何設計更好的 $q\displaystyle q$ 呢？一個辦法是設計一個 $t(x′∣x)\displaystyle t( x'|x)$ ，他從上一次接受的樣本作為條件，然后經過轉換propose一個新的樣本，如果我們能夠約束這個轉換不會改變目標分布的，那么我們就能快速的找到 $p(x)\displaystyle p( x)$ 的樣本。這也正是MCMC的思想。

IMCMC

然而MCMC的方法少說也有成百上千種，如何統一起來呢？我們從另一個角度來考慮MCMC。首先MCMC最基本可以從離散馬爾科夫鏈說起，A是一個轉移矩陣， $π\displaystyle \pi$ 是一個向量，如果滿足下述公式

$πA=π\pi A=\pi$

則稱 $A\displaystyle A$ 的平穩分布是 $π\displaystyle \pi$ 。類似的在連續的情況下：

$∫t(x′∣x)p(x)dx=p(x′)\int t( x'|x) p( x) dx=p( x')$

我們希望找到一個轉換概率分布 $t(x′∣x)\displaystyle t( x'|x)$ ，且該轉換不會改變來自 $p(x)\displaystyle p( x)$ 的樣本的分布。事實上這個t不一定是隨機的，當他是確定的映射時，我們有 $t(x′∣x)=δ(x′?f(x))\displaystyle t( x'|x) =\delta ( x'-f( x))$ ,相當于 $x′=f(x)\displaystyle x'=f( x)$ 于是
$∫δ(x′?f(x))p(x)dx=p(x′)\int \delta ( x'-f( x)) p( x) dx=p( x')$

這意味著我們要找到一個映射函數，使得 $x′=f(x)\displaystyle x'=f( x)$ 并且他們的分布一樣，又因為根據分布變換公式

$px(x)=px′(x)px(x)=px′(f(x))∣df(x)dx∣=px(f(x))∣df(x)dx∣=px′(x)=px(f?1(x))∣df?1(x)dx∣p_{x}( x) =p_{x'}( x)\\ p_{x}( x) =p_{x'}( f( x))\left| \frac{df( x)}{dx}\right| =p_{x}( f( x))\left| \frac{df( x)}{dx}\right| \\ =p_{x'}( x) =p_{x}\left( f^{-1}( x)\right)\left| \frac{df^{-1}( x)}{dx}\right|$
（其實從上述公式中我們就能大致的能看出， $f(x)=f^{-1}(x)$ 這個條件的成立，這就是所謂的Involutive function，IMCMC就是圍繞如何設計這個f函數來做的）
不過這樣的f不好找，所以一般我們可以用以下這個：

$t(x′∣x)=δ(x′?f(x))min?{1,p(f(x))p(x)∣?f?x∣}?Paccept++δ(x′?x)(1?min?{1,p(f(x))p(x)∣?f?x∣})?Preject,\begin{aligned} t(x'|x)=\delta (x'-f(x))\underbrace{\min\biggl\{1,\frac{p(f(x))}{p(x)}\biggl|\frac{\partial f}{\partial x}\biggl|\biggr\}}_{P_{\text{accept}}} +\\ +\delta (x'-x)\underbrace{\biggl( 1-\min\biggl\{1,\frac{p(f(x))}{p(x)}\biggl|\frac{\partial f}{\partial x}\biggl|\biggr\}\biggr)}_{P_{\text{reject}}} , \end{aligned}$

我們一般可以通過拒絕采樣的方式來近似這個轉換概率(這部分我不是很理解他跟拒絕采樣的關系，沒看懂)。這意味著，當 $x′=f(x)\displaystyle x'=f( x)$ 時， $t(x′∣x)=min?{1,p(f(x))p(x)∣?f?x∣}\displaystyle t(x'|x)=\min\biggl\{1,\frac{p(f(x))}{p(x)}\biggl|\frac{\partial f}{\partial x}\biggl|\biggr\}$ ，而當 $x′=x\displaystyle x'=x$ 時， $t(x′∣x)=(1?min?{1,p(f(x))p(x)∣?f?x∣})\displaystyle t(x'|x)=\biggl( 1-\min\biggl\{1,\frac{p(f(x))}{p(x)}\biggl|\frac{\partial f}{\partial x}\biggl|\biggr\}\biggr)$ ，顯然，如果 $x～p(x)\displaystyle x\sim p( x)$ ，且 $x′=f(x)\displaystyle x'=f( x)$ ，我們寫的稍微簡單一點，那么 $px′(x′)=t(x′∣x)p(x)=p(f(x))∣?f?x∣=px(f?1(x))∣df?1(x)dx∣\displaystyle p_{x'}( x') =t( x'|x) p( x) =p(f(x))\biggl|\frac{\partial f}{\partial x}\biggl| =p_{x}\left( f^{-1}( x)\right)\left| \frac{df^{-1}( x)}{dx}\right|$ 。事實上，這等價于認為， $f\displaystyle f$ 滿足 $f(x)=f?1(x)\displaystyle f( x) =f^{-1}( x)$ ，這樣的函數稱為involutive function。然而，這個分布只是在兩個點之間反復橫跳，我們可以引入一個輔助變量 $v\displaystyle v$ ，我們隨便抽取 $v\displaystyle v$ 的樣本，從而得到不同的 $x′\displaystyle x'$ ，這樣就不會反復橫跳了。可以證明，引入一個輔助變量后，只要滿足 $f(x,v)=f?1(x,u)\displaystyle f( x,v) =f^{-1}( x,u)$ ，那么x的平穩分布仍然是目標分布。

$t (x^{'}, v^{'} ∣ x, v) p (x, v) = t (x, v ∣ x^{'}, v^{'}) p (x^{'}, v^{'}) .$

可以證明，t的邊緣分布仍然服從detailed balance
$t^(x′∣x)=∫dvdv′t(x′,v′∣x,v)p(v∣x)t^(x′∣x)p(x)=t^(x∣x′)p(x′).\begin{aligned} \hat{t} (x'|x)=\int dvdv't(x',v'|x,v)p(v|x) \end{aligned} \begin{aligned} \hat{t} (x'|x)p(x)=\hat{t} (x|x')p(x'). \end{aligned}$

通過引入輔助變量，我們可以設計出iMCMC的算法，只需保證f是involutive function。

MH

一個特例是 $f(x,v)=v,x\displaystyle f( x,v) =v,x$ ，他將 $x,v\displaystyle x,v$ 的位置交換了一下，這時候iMCMC等價于MH算法。顯然他的接收概率為

$P=min?{1,p(f(x,v))p(x)q(v∣x)}=min?{1,p(v)q(x∣v)p(x)q(v∣x)}P=\operatorname{min} \{1,\frac{p(f(x,v))}{p(x)q(v|x)} \}=\operatorname{min} \{1,\frac{p(v)q(x|v)}{p(x)q(v|x)} \}$

HMC

另外一個經典的例子是HMC算法，他通過構造一個復雜的函數 $f\displaystyle f$ ，從而使得接受率非常高，幾乎不會拒絕（實際做的時候就直接接受，不判斷拒不拒絕），但是 $q\displaystyle q$ 又很簡單，實際上啟發了很多基于神經網絡的MCMC就是在搞這個 $f\displaystyle f$ 。

HMC依賴于用leap-frog算子 $L\displaystyle L$ 來求解Hamiltonian dynamics的數值積分。對于 $L:[x(t),v(t)]→[x(t+ε),v(t+ε)]L:[x(t),v(t)]\rightarrow [x(t+\varepsilon ),v(t+\varepsilon )]$ ，

$v(t+ε/2)=v(t)?ε2?x(?log?p(x(t)))x(t+ε)=x(t)+ε?v(?log?p(v(t+ε/2)))v(t+ε)=v(t+ε/2)?ε2?x(?log?p(x(t+ε)))\left. \begin{array}{ l } v(t+\varepsilon /2)=v(t)-\frac{\varepsilon }{2} \nabla _{x} (-\operatorname{log} p(x(t)))\\ x(t+\varepsilon )=x(t)+\varepsilon \nabla _{v} (-\operatorname{log} p(v(t+\varepsilon /2)))\\ v(t+\varepsilon )=v(t+\varepsilon /2)-\frac{\varepsilon }{2} \nabla _{x} (-\operatorname{log} p(x(t+\varepsilon ))) \end{array}\right.$

后 $F\displaystyle F$ 則是一個簡單的將輔助變量取負號的函數 $F:[x,v]=[x,?v]\displaystyle F:[ x,v] =[ x,-v]$ ，可以論文中給出了證明 $FLk=(FLk)?1\displaystyle FL^{k} =\left( FL^{k}\right)^{-1}$ ，而且他的jacobian矩陣的行列式也是為1。里面還有很多算法，這里就不寫了，有興趣自己去看。

參考資料

Neklyudov, Kirill, et al. “Involutive MCMC: a Unifying Framework.” arXiv preprint arXiv:2006.16653 (2020).

How does the proof of Rejection Sampling make sense?

(ML 17.13) Proof of rejection sampling (part 1)

Murphy, Kevin P. Machine learning: a probabilistic perspective. MIT press, 2012.

Bishop, Christopher M. Pattern recognition and machine learning. springer, 2006.

總結

以上是生活随笔為你收集整理的MCMC算法大统一: Involutive MCMC的全部內容，希望文章能夠幫你解決所遇到的問題。

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

上一篇： MCMC算法
下一篇： python文件操作实验总结,[干货分享