這篇文章將介紹一個非常自然的MCMC方法，也就是unadjusted Langevin algorithm（ULA）（實在沒找到中文叫什麼）。 這個方法其實並不算是Metropolis-Hastings方法，下篇文章將會介紹他的改進版本，也就是Metropolis adjusted Langevin algorithm（MALA）。但是為了解釋一下為什麼我覺得這個演算法很自然，還是先來回顧一下Metropolis-Hastings algorithm。

Metropolis Hastings Algorithm：

給定初始分佈

\nu

，

X_0\sim\nu

，n=0

n=n+1

Y_n\sim p(X_{n-1},\cdot)

，

U_n\sim Unif(0,1)

，

U_n

獨立於

Y_n

\begin{cases} If\ U_n<\alpha(X_{n-1},Y_n)\implies X_n:=Y_n\qquad&(accept)\\ else\ X_n:=X_{n-1}&(reject) \end{cases}

goto step 2

這裡acceptance probability alpha 的最佳選擇是

\alpha(x,y)=\min\{1,\frac{\mu(dy)p(y,dx)}{\mu(dx)p(x,dy)}\}\qquad\qquad\qquad\qquad（1）

對於Metropolis Hastings algorithm， 一個很重要的問題就是該如何選擇proposal kernel p， 因為其實只要選擇好p，acceptance probability 就可以直接根據式子（1）得到。在前幾篇文章裡面我們看到兩種不同的proposal kernel， 一個是independence sampling， 也就是我們選取一個與我們所處位置無關的kernel p，

p(x,dy)=\mu_0(dy)

。

另一個是random walk Metropolis， 也就是我們選取

p(x,dy)\sim N(x,h\cdot 1_d)

。

那這兩種方法的思路其實有點像是，我們得到一個目標分佈

\mu

，然後我們想方法找到一個合適的proposal kernel然後透過accept 和reject 來調整，最後得到一個以

\mu

為invariant measure 的馬氏鏈（其實是一個滿足detailed balance的馬氏鏈，也就是一個關於目標測度reversible的MC，事實上我們也可以構建一個non-reversible的MC，使得目標測度是他的invariant measure，因為畢竟detailed balance condition是強於invariant的，這個我們之後會寫到）。

那現在換一種思路，我們反著來。假設我們有一個目標測度

\mu

，我們可以先構造一個連續的馬可夫過程，使得該過程的invariant measure是目標測度，然後我們在將這個連續的過程離散化成為馬氏鏈， 這個思想就是Langevin algorithm。

假設我們的目標分佈為

\mu(dx)\propto e^{-U(x)}dx,\qquad U(x)\in C^2(\mathbb R^d)

。

那我們現在來構造一個關於目標分佈reversible的馬可夫過程。我們先來定義一個Dirichlet form

\mathcal E(f,g):=\frac 1 2\int \triangledown f(x)\cdot\triangledown g(x)\mu(dx)\qquad\forall f,g\in C^\infty_0(\mathbb R^d)

然後我們得到他的generator是

Lf=\frac 12 \triangle f-\frac 1 2\triangledown U\cdot\triangledown f

（這裡就不寫L的domain了，因為只是看一下這個馬可夫過程長什麼樣子）。根據SDE和martingale problem的關係，我們得到關於

\mu

reversible的過程應該是下面這個sde的解

dX_t=-\frac 1 2 \triangledown U(X_t)dt+dB_t\qquad\qquad (2)

接下來的問題就是這個sde是否存在一個non-explosive的解，要得到這樣的一個解我們需要一些關於U的額外條件，例如如果U滿足

\liminf_{|x|\to\infty}\frac{x\cdot \triangledown U(x)}{|x^2|}>-\infty

那麼式子（2）就存在一個non-explosive的解並且它的invariant measure就是

\mu

。並且在很多實際情況中，這個解的transition kernel滿足

W^2(\nu p_t,\mu)\le e^{-Kt}W^2(\nu,\mu),\nu\ arbitrary\ initial\ distribution.

也就是他會exponential 趨近於equilibrium。

但是呢，其實這個並沒啥用，因為這是一個連續的過程，計算機並不能模擬一個連續的過程，所以我們必須把這個連續的馬可夫過程離散化成為一個馬氏鏈。

這裡我們可以用Euler-Maruyama schema來離散化這個，也就是說

X_{n+1}=X_n-\frac h 2\triangledown U(X_n)+\underbrace{B_{(n+1)h}-B_{nh}}_{=\sqrt h Z},\qquad Z\sim N(0,I_d)

。

也就是

p(x,dy)\propto e^{-\frac{1}{2h}|y-x+\frac h 2 \triangledown U(x)|^2}

。

我們這裡可以看到和random walk Metropolis 的區別，也就是我們這裡還添加了第二項關於U的gradient的資訊，當然這樣可以增加我們的精確性，但其實也會造成困難，因為我們每一步都要計算U的gradient，在高維的情況下會是特別expensive的。

那我們現在就得到

Unadjusted Langevin algorithm:

X_0\sim \nu,\nu\ initial\ distribution

for n=1，2，。。。

X_{n}\sim N(X_{n-1}-\frac h 2 \triangledown U(x),h\cdot I_d)

我們可以看到在這個演算法裡面，我們並沒有使用accept和reject，所以這個演算法其實算不上是Metropolis-Hastings 演算法。

因為我們其實是離散化了一個以目標測度

\mu

為invariant measure的過程，所以這個方法有個硬傷（systematic error），也就是透過這個方法得到的馬氏鏈，它的invariant measure並不是

\mu

，而是另一個measure，所以真的是個硬傷啊。我們來看個簡單的例子。

假設我們想模擬的是高斯分佈，也就是

\mu\sim N(0,I_d)

，那我們得到

U(x)=\frac 1 2 |x|^2

。那我們將得到

X_{n+1}=(1-\frac h 2)X_n+\sqrt h Z

我們來求一下這個過程的invariant measure， 我們假設他的invariant measure是

N(0,v\cdot I_d)

，那我們將得到

X_{n+1}\sim N(0,[(1-\frac h 2)^2v+h]\cdot I_d)

所以我們有

(1-\frac h 2)^2v+h=v\implies v=1+\frac{h}{4-h}

。

也就是其實我們模擬的measure是

\mu_h:=N(0,(1+\frac{h}{4+h}))\neq N(0,I_d)=\mu

。但是呢，其實也沒那麼差，因為如果我們讓h趨近於0，我們還是能得到的，現在的問題是，我們該讓h以多快的速度趨近於0呢，那我們再來算一下他和目標測度之間的距離：

W^2(\mu,\mu_h)\approx \frac{h^2d}{64}

，

也就是我們應該選取

h\sim O(\frac {1}{\sqrt d} )

。

在random walk metropolis我們得到h應該是

O(\frac 1 d)

，也就是說這個方法即使沒有得到我們想要的measure，但是他趨近於我們想要的measure的速度還是要遠遠快於random walk metropolis的。

事實上我們還能證明很多其他這個演算法的優點的，但是就不在這裡繼續寫了，反正就是這個方法還是要好於random walk Metropolis的，就更不用說independence sampling了。

這篇文章我們介紹了unadjusted Langevin algorithm（ULA）， 他不是一個Metropolis-Hastings algorithm， 因為我們在這裡並沒有使用accept和reject。在選取步長的時候我們可以考慮

h\sim O(\frac{1}{\sqrt d})

， d是dimension，所以它是快於我們先前看到的兩個。但是他的硬傷就是模擬的並不是我們想要的measure，也就是存在systematic error。在下一篇文章我們將介紹Metropolis adjusted Langevin algorithm（MALA），也就是我們將在ULA的基礎上加上accept和reject以去除掉systematic error。