白話馬爾科夫鏈蒙特卡羅方法(MCMC)

前言

你清茶園不是人待的地方！
里面的個個都是人才，說話又好聽——就是我太菜了啥也聽不懂，這次期中還考的賊**爛，太讓人郁悶了。
最近課上講這個馬爾科夫鏈蒙特卡羅方法，我也學得一塌糊塗。這時我猛然想起了自己的博客園密碼（霧），來更個博客吧。

[Warning] 本人數學水平差勁，下文用詞不嚴謹、缺少部分證明，請酌情閱讀。若出鍋，歡迎指正。

啥是馬爾科夫鏈？

馬爾科夫鏈(Markov Chain)，簡單來說就是一個用來隨機游走的有向圖，每條邊(u, v)的邊權\(p_{uv}\)代表“當前在u，下一步走到v”的概率，顯然需要

\[p_{uv}\ge 0, \sum_{v}p_{uv}=1. \]

下文中我們假設這個有向圖是強連通的，即任取兩個點u和v，都存在從u到v的、邊權都大於0的路徑（當然從v到u的路徑也要存在）。

馬爾科夫鏈（也就是這個隨機游走過程）的美妙性質在於它收斂。怎么個收斂法呢？

設這個圖有n個點，令\(n\)維行向量\(\mathbf p(t)\)表示隨機游走了t步之后的概率分布（在時間t分別位於每個點的概率），\(\mathbf p(t)_i\) 就是第t步到點i的概率。初始狀態\(\mathbf p(0)\)是隨便欽定的。

再定義一個量\(\mathbf a(t)\)，名叫“長期平均概率分布(long-term average probability distribution)”，

\[\mathbf a(t) = \frac1t \sum_{i=0}^{t-1} \mathbf p(i). \]

顧名思義，就是把前\(t\)個時間點的概率分布取個平均。

定理（馬爾科夫鏈基本定理）

• 存在唯一概率分布向量\(\mathbf \pi\)使得\(\mathbf \pi P = \mathbf \pi\)。（這個P就是由\(p_{xy}\)構成的矩陣。）
• \(\lim_{t\rightarrow \infty} \mathbf a(t)\) 存在，而且就等於\(\mathbf \pi\)。

(證明待補充……放假再敲？）

這個定理就很令人開心了——不管欽定初始狀態\(\mathbf p(0)\)的你是歐皇還是非酋，只要游走足夠多步，\(\mathbf a(t)\)肯定會收斂到唯一的答案。

誒？為啥要定義這個a(t)呢？p(t)自己它不收斂么？
還真不收斂。考慮一下\(w_{12} = 1, w_{21}=1, p(0)=(1, 0)\)。在這個圖上游走就是在兩個點之間反復橫跳，p(t)是不會收斂的！但是a(t)就能收斂到\((\frac12, \frac12)\)。

啥是蒙特卡羅方法？

蒙特卡羅(Monte Carlo)是誰？不是誰，這是個賭場的名字 (=_=|||)。取這個名字大概只是因為……它是隨機的，賭場也是隨機的？不得不說這個洋氣名字實在太勸退了。

其實蒙特卡羅方法就是……抽樣估計。小學學的撒豆子求面積啊，蒲豐投針計算圓周率啊，都可以視作蒙特卡羅算法。

啥是馬爾科夫鏈蒙特卡羅方法(MCMC)？

I have a Markov Chain, I have a Monte Carlo, ah, Markov Chain Monte Carlo!

我們已經知道，馬爾科夫鏈是個好東西，它保證能收斂到某個概率分布。現在我們已知一個概率分布，想要構造出相應的馬爾科夫鏈。這有什么用呢？有些時候，概率分布長得比較復雜，直接根據它生成隨機變量非常困難，例如想在一個高維空間中的凸多邊形中隨機取一個點——這時候我們就可以先構造出這個概率分布對應的馬爾科夫鏈，然后在上面隨機游走來取點。

不同的馬爾科夫鏈可能收斂到相同的概率分布，但是收斂的速率有快有慢。在應用中，我們肯定希望馬爾科夫鏈收斂得快一點，少游走幾步就能獲得和想要的概率分布差不多的結果。下面是兩個比較好用的構造方法。

Metropolis-Hasting 算法

有的時候又叫Metropolis算法，反正是個一個名字挺長的算法。以下簡稱MH算法。

和它不好記的名字相反，MH算法描述起來非常簡單、非常符合直覺！

首先你要有一個無向連通圖。先不管這個圖是怎么構建的，很有可能是出題人送給你的（假如你在做相關的練習題的話）。

設概率分布是\(\mathbf p\)，點\(i\)的概率是\(p_i\)。假如你當前在點\(i\)：

• 先從與點\(i\)相鄰的所有點中，等概率隨機取一個點\(j\);
• 如果\(p_j \ge p_i\)，則立刻走到\(j\)；
• 如果\(p_j < p_i\)， 隨機一下，有\(\frac{p_j}{p_i}\)的概率走到\(j\)，否則待在原地不動。

這樣\(p_j\)更大的\(j\)更有可能被走到，很符合直覺吧！

把上面的文字策略用公式表示一下，馬爾科夫鏈中點\(i\)到點\(j\)的邊權\(p_{ij}\)就是

\[p_{ij} = \frac1r \min(1, \frac{p_j}{p_i}), \]

而剩下的概率都分給了\(p_{ii}\)，

\[p_{ii} = 1 - \sum_{j\ne i} p_{ij}. \]

那么這個馬爾科夫鏈到底能不能收斂到\(\mathbf p\)這個概率分布呢？

引理

令\(\mathbf \pi\)為一個概率分布，\(P\)為馬爾科夫鏈的邊權矩陣，如果對任意兩點\(x,y\)，均有

\[\pi_x p_{xy} = \pi_y p_{yx}, \]

則\(\mathbf \pi\)就是馬爾科夫鏈收斂到的概率分布。

證明：若\(\pi_x p_{xy} = \pi_y p_{yx}\)，那么枚舉\(y\)，對等式兩邊分別求和，就有\(\pi_x = \sum_y \pi_y p_{yx}\)，即\(\mathbf \pi = \mathbf \pi P\)。根據上面的馬爾科夫鏈基本定理可知收斂到\(\mathbf \pi\)。

把MH算法構造的馬爾科夫鏈代入這個引理，可知確實收斂到\(\mathbf p\)！

（課上講了一個應用MH算法破譯密碼的有趣例子，我來不及敲完了，可以看看這篇文章的開頭。）

Gibbs 采樣

Gibbs 采樣(Gibbs Sampling)是另一種MCMC，它使用的是一種特殊的無向圖：每個點對應\(d\)維空間中的一個格點\(\mathbf x = (x_1, x_2, \cdots, x_d), x_i\in\mathbb Z\)，如果\(\mathbf x\)和\(\mathbf y\)只有一個個坐標不同，則二者之間有連邊。這樣形成的就是一個類似\(d\)維網格（lattice）的結構，但每條（與坐標軸平行的）直線上的點都形成一個團。

（一個用來Gibbs取樣的無向圖，圖片來自Blum, Hopcroft, Kannan, "Foundations of Data Science"）

對於相鄰的兩點\(\mathbf x\)和\(\mathbf y\)，不失一般性，假設它們第一維坐標不同（\(x_1 \ne y_1\)）而剩下的坐標都相同。那么就令

\[p_{\mathbf{xy}} = \frac1d p(y_1|x_2, x_3, \cdots, x_d). \]

而剩余的概率留給\(p_{\mathbf{xx}}\),

\[p_{\mathbf{xx}} = 1-\sum_{\mathbf y \ne x}p_{\mathbf{xy}}. \]

驗證：

\[p(x_1|x_2, x_3, \cdots, x_d) p_{\mathbf{xy}} = p(y_1|x_2, x_3, \cdots, x_d) p_{\mathbf{yx}}, \]

故根據條件概率公式有

\[p(x_1) p_{\mathbf{xy}} = p(y_1) p_{\mathbf{yx}}, \]

所以概率分布確實收斂到\(\mathbf p.\)

ε-混合時間(ε-mixing time)

（我沒查到這個ε-mixing time怎么翻譯成中文……有大佬知道的話評論區告訴我一下qaq）

盡管我們知道馬爾科夫鏈早晚會收斂的，但它到底早收斂還是晚收斂，對我們很重要。為了衡量收斂的快慢，定義ε-混合時間(ε-mixing time)為\(t\)的最小值，滿足：對任意初始狀態\(\mathbf p(0)\),\(\|\mathbf a(t) - \mathbf \pi(t)\|_1 < \epsilon\)（這個范數就是曼哈頓距離那種，各維度距離直接相加的和）。為了保證算法的效率，我們想知道這個ε-混合時間的一個上界。

一個上界是由導率(conductance)確定的。啥是導率呢？顧名思義，就和物理上的電導率差不多。對於節點集合\(S\)，令\(\pi(S) = \sum_{x\in S} \pi_x\)。令\(\bar S\)為\(S\)的補集，若\(\pi(S) \le \pi(\bar S)\)，那么導率\(\Phi(S)\)就是

\[\Phi(S) = \frac{\sum_{(x, y)\in(S, \bar S)}\pi_xp_{xy}}{\pi(S)} = \sum_{x\in S}\frac{\pi_x}{\pi(S)}\sum_{y\in\bar S}p_{xy}, \]

即“已知當前在\(S\)中，下一步跳出\(S\)的概率”。

直覺上，如果\(\Phi(S)\)都很大，那么我們可以在圖上來去自如，收斂就比較快；反之，如果某個\(\Phi(S)\)很小，一旦進入\(S\)就被困住、很難逃脫，收斂就可能很慢。
因此，定義整個馬爾科夫鏈的導率為

\[\Phi = \min_{S\subset V, S\ne \emptyset} \Phi(S). \]

定理

任何無向圖隨機游走的ε-混合時間有上界

\[O\left(\frac{\ln(1/\pi_{min})}{\Phi^2\epsilon^3}\right). \]

（證明……敲不完……下次也不一定）

例子

一條\(n\)個節點順序組成的鏈，\(\mathbf \pi = (\frac1n, \frac1n, \cdots, \frac1n)\)。

可以構造馬爾科夫鏈\(p(1, 1) = p(n, n) = \frac12, p(i, i+1) = p(i + 1, i) = \frac12, \forall 1 \le i \le n-1\)。

最難逃脫的\(S\)就由前半條鏈（后半條也行）構成的集合，\(\pi(S) = \frac12\),

\[\Phi = \Phi(S)= \frac{\frac1n \cdot \frac12}{\frac12} = \frac1n, \]

ε-混合時間的上界就是

\[O\left(\frac{\ln n}{\left(\frac1n\right)^2 \epsilon^3}\right) = O\left(\frac{n^2\ln n}{\epsilon^3}\right). \]

---

寫作業去啦 未完待續_(:з」∠)_