梅特羅波利斯修正朗之萬演算法

在計算統計學中，梅特羅波利斯修正朗之萬演算法（，簡稱MALA）或朗之萬蒙地卡羅方法（，簡稱LMC）是一種馬可夫鏈蒙地卡羅（MCMC）方法，用於從難以直接抽樣的機率分佈中獲取——即隨機觀測值的序列。顧名思義，MALA結合了兩種機制，用以生成具有目標機率分佈作為的隨機漫步狀態：

• 利用（過阻尼）朗之萬動力學提出新狀態，該動力學使用目標機率密度函數梯度的計算值；
• 再透過梅特羅波利斯-黑斯廷斯算法接受或拒絕這些提案，該演算法使用目標機率密度的計算值（而非其梯度）。

非正式地說，朗之萬動力學以梯度流的方式驅動隨機漫步朝向高機率區域，而梅特羅波利斯-黑斯廷斯算法接受/拒絕機制則改善了此隨機漫步的混合與收斂特性。MALA最初由於1994年提出（儘管智慧蒙地卡羅方法早在1978年便已問世），其特性則由與及共同進行了詳細探討。此後已衍生出許多變體與改進版本，例如吉羅拉米（Girolami）與卡爾德海德（Calderhead）於2011年提出的流形變體。該方法相當於僅採用單一離散時間步長的（混合蒙地卡羅）演算法。

細節
設 \pi 為定義於 \mathbb{R}^{d} 上的機率密度函數，我們希望從中抽取一組獨立同分布的樣本。我們考慮過阻尼的朗之萬伊藤過程：

:\dot{X} = \nabla \log \pi(X) + \sqrt{2} \dot{W},

其驅動力來自標準布朗運動 W 的時間導數。（注意此擴散過程的另一種常用正規化形式為

:\dot{X} = \frac{1}{2} \nabla \log \pi(X) + \dot{W}

其產生的動力學與上述相同。）當 t \to \infty 時，此 X(t) 的機率分佈 \rho(t) 趨近於一個穩態分佈，該分佈在擴散作用下亦保持不變，我們將其記為 \rho_\infty。事實上，可知 \rho_\infty = \pi。

朗之萬擴散的近似樣本路徑可透過多種離散時間方法生成。其中最簡單的方法之一是採用固定時間步長 \tau > 0 的歐拉-丸山法。我們設定 X_0 := x_0，然後透過以下公式遞迴定義真實解 X(k \tau) 的近似值 X_k：

:X_{k + 1} := X_k + \tau \nabla \log \pi(X_k) + \sqrt{2 \tau} \xi_k,

其中每個 \xi_{k} 皆為從 \mathbb{R}^{d} 上的多變量常態分布獨立抽取的樣本，其均值為0，且共變異數矩陣等於 d \times d 單位矩陣。注意X_{k + 1} 服從常態分佈，其期望值為 X_k + \tau \nabla \log \pi(X_k)，且共變異數等於 2 \tau 倍的 d \times d 單位矩陣。

相較於用於模擬朗之萬擴散的歐拉-丸山法，該方法總是根據以下更新規則更新 X_k：

:X_{k + 1} := X_k + \tau \nabla \log \pi(X_k) + \sqrt{2 \tau} \xi_k,

MALA則加入了額外的一步。我們將上述更新規則視為定義了一個針對新狀態的「提案」 \tilde{X}_{k + 1}，

:\tilde{X}_{k + 1} := X_k + \tau \nabla \log \pi(X_k) + \sqrt{2 \tau} \xi_k

此提案將根據梅特羅波利斯-黑斯廷斯演算法予以接受或拒絕：設定

:\alpha := \min \left\{ 1 , \frac{\pi(\tilde{X}_{k + 1}) q(X_{k}\mid\tilde{X}_{k + 1})}{\pi({X}_{k}) q(\tilde{X}_{k + 1}\mid X_k)} \right\},

其中

:q(x'\mid x) \propto \exp \left( - \frac{1}{4 \tau} \| x' - x - \tau \nabla \log \pi(x) \|_2^2 \right)

是從 x 到 x' 的遷移機率密度（注意一般而言 q(x'\mid x) \neq q(x\mid x')）。 設 u 取自區間 [0, 1] 上的連續型均勻分布。若 u \leq \alpha，則接受該提案，並設定 X_{k + 1} := \tilde{X}_{k + 1}；否則，拒絕該提案，並設定 X_{k + 1} := X_k。

朗之萬擴散與梅特羅波利斯-黑斯廷斯演算法結合後的動力學，滿足了存在唯一、不變且穩態的分布 \rho_{\infty} = \pi 所必需的細緻平衡條件。相較於原始的梅特羅波利斯-黑斯廷斯演算法，MALA的優勢在於它通常會建議向 \pi 機率較高的區域移動，而這些移動被接受的可能性也更高。另一方面，當 \pi 呈現強烈的各向異性（即在某些方向上的變化遠比其他方向更快）時，必須取 0 才能正確捕捉朗之萬動力學；使用一正定矩陣 A \in \mathbb{R}^{d \times d} 可透過以下方式生成移動建議，從而緩解此問題：
:\tilde{X}_{k + 1} := X_k + \tau A \nabla \log \pi(X_k) + \sqrt{2 \tau A} \xi_k,
使得 \tilde{X}_{k + 1} 的期望為 X_k + \tau A \nabla \log \pi(X_k)，且共變異數為 2 \tau A。

對於特定類別的目標分佈，可以證明此演算法的最佳接受率為 0.574；若實際情況發現有顯著差異，則應相應調整 \tau。

參考資料