馬可夫鏈蒙地卡羅: MCMC

8 min readOct 26, 2023

<圖學玩家第030篇原創文>

Introduction

在理解Diffusion Model時，除了DreamFusion Part 1這篇的視角外，還能透過Score-Based Generative Model與Energy-Based Model的角度去理解。而馬可夫鏈蒙地卡羅 (MCMC) 在理解Energy-Based Model時，會是不可或缺的角色。

另外如貝氏推斷中求期望值的積分部分，也很常會需要用到MCMC這類蒙地卡羅的方法進行抽樣以進行估計。

筆者假設讀者對基本的蒙地卡羅方法已有一定的認識，在此循序漸進地介紹幾種非均勻採樣，以協助讀者理解MCMC。

非均勻採樣

在計算積分I = ∫ f(x) dx的時候，積分區域可能是任意形狀。此時若使用MC積分(蒙地卡羅法)中的平均值法去均勻採樣並估計I，就可能造成採樣到的f(x)接近於0，進而導致效能上的浪費。

因此有了一系列非均勻採樣的方法，本文主要介紹接受-拒絕抽樣，重要採樣以及馬可夫鏈蒙地卡羅 (MCMC)。

接受-拒絕抽樣

接受-拒絕抽樣(Accept-Reject Sampling)的精神在於，提出一個分布 q(x)，其與欲採樣的目標分布p(x)之關係為p(x) < M * q(x)。其中M為某非負常數，也是我們設定的參數。

接著進行如下算法:

for i = 1 to N:
    // 產生一個均匀分布的隨機變量u，和服從q(x)的隨機變量x
    if u < p(x)/ M*q(x)，
       那麼x就被選為一樣本（accept）
       i+=1
    else：reject

這個概念可以用下圖直觀的解釋，虛線的部份可以想成M * q(x)，藍線的部分為p(x)。

這個算法的限制在於，不一定總能找到p(x) < M * q(x)這樣的關係，以及如果M設太大的话，抽樣的速度會很慢。

重要採樣

重要採樣 (Importance Sampling)的核心觀念就是，在分布機率π(x)高的地方多採樣些點，分布機率π(x)小的地方少採樣些點。

首先我們知道函數f(x) 在概率分布π(x)下的期望為:

由於π(x)我們無法求得。因此我們需要借助一個方便採樣的概率分布 p(x) 。這樣一來當我們在 p(x) 上採樣 {x1, x2, …, xi} 後，f(x)的期望為:

接著我們做以下式子的轉換

這樣一來，當我們在p(x)採樣 {x1, x2, …, xi} 後，f(x)的期望即為:

其中π(xi) / p(xi)稱為重要性權重。由於p(x)是我們提出的proposal，而這樣做的原因就是，我們預期π(x)在某些區間的值會大一點，因此設計對應的p(x)在該處多採樣些; 與之相對的，權重就會減少些。

事實上，在Rendering Equation中的積分部分，也是採用Importance Sampling這種蒙地卡羅法進行估計

馬可夫鏈蒙地卡羅 (MCMC)

馬可夫鏈

如下圖，馬可夫鏈被定義為狀態空間中從一個狀態到另一個狀態轉換的隨機過程。而狀態轉移的過程可以用矩陣P來表示，Pij即表示從第i個狀態轉移到第j個狀態的概率。

其中很重要的馬可夫性質是，當下的狀態僅與前一刻的狀態有關，如下表示:

平穩分布

如果一個非周期的馬可夫鏈有狀態轉移矩陣P，且任何兩個狀態是連通(從任意一個狀態可以通過有限步到達其他的任意一狀態)的，則我們有以下結論:

我們假設馬可夫鏈中的每個狀態都有其對應機率分布，則π(j)指的是狀態j所對應的機率分布。上式的意思也就是說，無論當前的狀態i為何，轉換到狀態j的機率都會固定為π(j)。且π(j)也可以被表示為:

按照上面描述，整個轉移矩陣P就會長的像下圖

而下述π = [π(1), π(1) …, π(j),…π(K)]即稱為該馬可夫鏈的平穩分布。

Detailed Balance

Detailed Balance的定義如下，只要滿足Detailed Balance，馬可夫鏈就會有收斂到平穩分布的性質

其詳細證明可以參考這篇。

MCMC

那麼說這麼多究竟要幹嘛?下圖其實可以解釋我們要做的事。左方可以看出，我們在Importance Sampling所提出的Proposal Q(x)，無法達到在P(x)機率較高的地方進行較多的採樣的目的。

因此我們希望可以透過Adaptive的方式，讓Proposal Q(x' | x)逐步在P(x)各個地方進行採樣，而這就是為何我們引進馬可夫鏈。

讀者必須先改變一下對於前面馬可夫鏈介紹中離散狀態的觀念。在連續的分布P(x)上，每個x都可以被視為一個狀態。

如上圖右側，MCMC的步驟就是先隨機初始一個x0 when t = 0，並帶入
Q(x1 | x0)這個條件分布以得到x1 when t = t+1，並接續進行下去。

在連續分布下，Q(x’ | x)就是離散馬可夫鏈中的轉移矩陣。通常Q(x’ | x)如上圖右側取一常態分佈即可。

接著我們假設P(x)分布就是某個馬可夫鏈最終的穩定分布。按照前面章節所述，我們所用的Proposal Q(x’ | x)必須滿足Detailed Balance，然而大多數的情況下這並不可能

這時候MCMC就引進了一個接受率α(x’ | x)，令其滿足Detailed Balance。

而該接受率分布也很好得到:

把Q(x’| x) α(x’ | x) = T(x’| x) ，我們就可以得到Detailed Balance條件

這裡的接受率α(x’ | x)與前面章節所述的接受-拒絕抽樣很類似，這裡的思維是一個常見的條件分布 Q(x’ | x) 通過一定的接受-拒絕概率得到目標條件分布T(x’ | x) 。

MCMC Algorithm

初始化起始狀態x; 並設計條件機率分布Q(x’ | x)，通常採常態分佈即可。
設定目標分布P(x)
建構一條長度為 n 的馬可夫鏈 : 每次迭代遵照以下步驟 :
3–1. 針對當前狀態 x，從Q(x’ | x)得出提議狀態x'
3–2. 從均勻分布中取樣 u ~ U(0, 1)，如果
u < α(x' | x) = P(x’) * Q(x | x’) ，則轉移狀態至 x’，否則停留在 x
3–3. 迭代 n 次，得出的馬可夫鏈即為目標分布中的採樣結果

然而在第3-2中，由于α(x’ | x)可能非常的小，比如0.1，導致大部分的採樣值都被拒絕轉移，效率很低。有可能我們採樣了上百萬次馬可夫鏈都没有收斂，一直停某個狀態x。

MCMC Improvement

對於上述接受率過低的問題，M-H採樣有提出對應的解法。對於M-H採樣以及為解決多維度資訊的Gibbs採樣有興趣的讀者，可以再另外閱讀如這篇文章或Google一下相關資訊，本篇文章僅就MCMC的主要精神做介紹。

如果喜歡筆者的文章分享，可以幫忙追蹤圖學玩家，你們的閱讀與追蹤是筆者繼續為大家分享的動力~