HMC:哈密頓蒙特卡洛方法

哈密頓蒙特卡洛方法

1. 簡介

哈密頓蒙特卡洛方法（Hamiltonian Monte Carlo）這一算法最早出現在由Dunne等人在1987年發表的一篇論文中，最初他們稱之為混合蒙特卡洛（Hybird Monte Carlo），是因為該方法將分子動力學理論和馬爾科夫鏈蒙特卡洛方法（Markov Chain Monte Carlo）混合在了一起。但Hamiltonian MC這個稱呼更能具體的描述所引入的Hamiltonian動力學，因而更廣為人知。HMC在統計問題上的應用始於（Neal,1996a）用於神經網絡模型的參數估計。在Neal看來，這是一個被統計學家低估的算法，時至今日，對於高維的參數分布估計問題，HMC仍然是最好的方法，如貝葉斯深度學習中使用HMC推斷網絡參數。
掌握HMC方法需要把握三個關鍵：

• Metropolis-Hasting算法
• 哈密頓動力學基礎
• 蛙跳法積分

2. 准備

2.1 Metropolis-Hasting算法

本文假設讀者擁有一定的相關統計知識，更詳細的馬爾可夫鏈蒙特卡洛內容等可以參考

1. 統計學習方法 李航
2. 隨機過程基礎
3. 博客：基礎馬爾可夫鏈概念，MCMC算法。

簡單來說，MH算法是MCMC算法的一種易於實現的改良版本。在回顧MH算法之前，這里對Markov Chain相關概念也進行簡要介紹。馬爾可夫鏈是在參數空間中通過序列地應用被稱作馬爾可夫轉移的隨機映射而產生一系列連續點集。一個任意的馬爾可夫鏈可以簡單地在參數空間中漫游，它本身在計算期望是沒有特殊作用的，然而當它漫游次數足夠，接近並保持為目標分布，這就意味着該鏈是通過目標分布生成的一個樣本集合。而且后續的漫游產生的樣本集合也任然保持在這個目標分布當中。這樣的樣本屬於目標分的典型集（Typical set，可以類似理解為大部分的樣本所在的集合中，如正態分布的\(\mu\pm2\sigma\)區間）。因此，理論上無論從參數空間的什么位置開始，相應的馬爾可夫鏈最終都會進入並且跨越目標分布的典型集。通過這些樣本集進行期望計算 \(E(\widehat{f})\)，以及估計整個參數空間的期望 \(E(f)\)。

圖 1 綠色的部分代表馬爾可夫轉移密度，它決定了新點產生的概率

MCMC方法定義了在給定目標分布的馬爾可夫轉移情況下量化典型集的一般策略。然而如何構建這樣的轉移本身就是一個不小的問題。幸運的是，對於任意給定的目標分布，都有各種自動構建合適的轉移的程序，其中最重要的就是MH算法。
MH算法分兩步：建議（proposal）與 修正（correction）。這里的建議是指給初始狀態施加一個任意的隨機擾動，而修正指的是拒絕任何偏離目標分布的典型集太遠的的建議。更正式的說，假設 \(Q(q^*|q)\)為定義的每一個建議的概率密度，接受給定的建議的的概率則是：

\[a(q^*|q)= \min(1, \frac{Q(q|q^*)\pi(q^*)}{Q(q^*|q)\pi(q)}) \]

\(\pi\)表示目標分布，\(q\)是感興趣的變量，一般來說，采用高斯分布作為建議，這一算法至今仍然被廣泛使用。即：

\[Q(q^*|q)=\mathcal{N}(q^*|q,\Sigma)) \]

采用高斯擾動的好處是建議機制在初始點和建議點是對稱等價的，因此它的接受率可以簡化成;

\[a(q^*|q)= \min(1, \frac{\pi(q^*)}{\pi(q)}) \]

隨機游走的MH算法實現簡單，非常符合直覺，建議點的轉移趨向於概率大的鄰域（若是從后續提到的正則分布中的能量的角度來說，建議點的轉移則是優先能量低的點）。這樣的 建議 然后修正 的組合程序優先選擇了落在高概率質量的點，並集中於典型集之內。但隨着目標分布的維度和復雜性的增加，這一算法的性能大打折扣。盡管可以通過調整建議分布的尺度來獲得更高的接受率，但會導致馬爾可夫鏈移動緩慢，甚至是失敗的探索，得到高度偏置的估計量。如果需要擴展到高維概率分布中去，那么就需要一個更好的方法來探索典型集，特別是需要利用好典型集本身的集合性質，HMC就是這樣的方法。

2.2 哈密頓動力學簡介

哈密頓動力學（Hamiltonian Dynamics）是拉格朗日力學的勒讓德變換。它可以從一個直觀且經典的物理解釋來簡單理解。拉格朗日運動方程是基於廣義坐標和廣義速度的，而哈密頓運動方程則是基於廣義坐標和廣義動量的。想像一個冰球在一個光滑表面上滑動，這個系統的能量狀態包括冰球的勢能 \(U\) 和動能 \(K\)。冰球的勢能 \(U(q)\) 與當前的位置 \(q\) 有關，它的動能 \(K(q)\) 與動量 \(p\) 有關。當冰球遇到上坡會減速，動能減少，勢能增加；當它遇到下坡會加速，動能增加，勢能減少。在非物理的MCMC應用中，HMC中的位置對應於感興趣的變量。勢能等於這些變量的負對數概率密度。每一個位置變量都會對應一個動量變量，動量由人為引入。這一物理先驗的引入，最大的意義在於在MCMC方法中(更確切的說MH算法)，通過動量項避免了采樣點完全的隨機游走。

2.2.1 哈密頓方程

哈密頓動力學作用於一個 \(d\)維位置向量 \(q\)和一個 \(d\) 維動量向量 \(p\) ，因此整個狀態空間是2維的。這個系統由一個成為哈密頓量的函數 \(H=(q,p)\)。
它的運動方程描述如下：

\[\begin{align} &\frac{dq_{i}}{dt}= \frac{∂H}{∂p_{i}} \tag{2.1}\\ &\frac{dp_{i}}{dt}=-\frac{∂H}{∂q_{i}} \notag \\ \end{align} \]

對於任何的時間間隔 \(s\),都定義了一個從任意時間狀態 \(t\) 到 \(t+s\) 的映射 \(T_{s}\)。作為替代，可以使用向量\(q\)和向量\(p\)構造一個新的二維狀態向量\(z=(p,q)\)，此時哈密頓方程可以寫為

\[\frac{dz}{dt}=J∇H(z) \]

\(∇H\) 是 \(H\) 的梯度（即 \([∇H]_{k}=∂H/∂z_{k}\)），並且

\[J=\left[\begin{array}{rl} 0_{d \times d} & I_{d \times d} \\ -I_{d \times d} & 0_{d \times d} \end{array}\right]\]

是一個由單位矩陣和零矩陣構建的 \(2d×2d\) 的矩陣。通常情況下，HMC當中可以將哈密頓方程寫成如下形式 $$H(q, p)=U(q)+K(p) \tag{2.2}$$ 這里的 \(U(q)\)被稱為勢能，並且被定義為我們想要采樣的參數 \(q\)對應分布的負對數概率密度加上一個常數。\(K(p)\)被稱作動能，通常被定義為

\[K(p)=p^{T}M^{-1}p/2 \tag{2.3} \]

這里 \(M\) 是一個對稱對角正定的“質量矩陣”，通常為單位矩陣的標量倍數。這種形式的的 \(K(p)\) 對應於一個均值為0，協方差矩陣為\(M\)的高斯分布的對數概率密度。 \(H\)，\(K\)在這種形勢下，方程\((2.1)\)可以寫成

\[\begin{align} \frac{d q_{i}}{d t} &=\left[M^{-1} p\right]_{i} \tag{2.4}\\ \frac{d p_{i}}{d t} &=-\frac{\partial U}{\partial q_{i}} \notag \end{align} \]

2.2.2 哈密頓動力學的性質

在構造MCMC更新時，哈密頓動力學的幾個性質是至關重要的，這些性質的在相關文獻中有詳細介紹：

1. 可逆性：保持期望的分布不變
2. 守恆性
3. 體積保持性：在進行MH更新時，不需要考慮積分空間體積變化，概率質量得以保證。
4. 辛性

2.3 離散哈密頓方程——蛙跳法（Leapfrog Method）

為了在計算機上模擬哈密頓動力學方程，它必須使用若干離散的小步長 \(\varepsilon\)的時間來求和近似。從時刻0開始計算 \(\varepsilon,2\varepsilon,3\varepsilon...\)。通常假設哈密頓量的形式為\((2.2)\)，雖然\(K(p)\)可以時任意形式，但為了簡單起見，假設動能形如\((2.3)\)。而且，\(M\)是由對角元素 \(m_{1},...,m_{d}\) 構成的對角矩陣，所以

\[K(p)=\sum^{d}_{i=1}\frac{p^2_{i}}{2m_{i}} \tag{2.5} \]

從方程 \((2.4)\)利用歐拉法可以很簡單的得到：

\[\begin{aligned} &p_{i}(t+\varepsilon)=p_{i}(t)+\varepsilon \frac{d p_{i}}{d t}(t)=p_{i}(t)-\varepsilon \frac{\partial U}{\partial q_{i}}(q(t)) \\ &q_{i}(t+\varepsilon)=q_{i}(t)+\varepsilon \frac{d q_{i}}{d t}(t)=q_{i}(t)+\varepsilon \frac{p_{i}(t)}{m_{i}} \end{aligned} \]

但是歐拉法容易發散，在分子動力學中的運動方程模擬主要還是使用蛙跳法，可以在精度和收斂方便做得更好，且蛙跳法積分在時間上是可逆的，這對於馬爾可夫鏈達到細致平衡條件更有利。

\[\begin{aligned} p_{i}(t+\varepsilon / 2) &=p_{i}(t)-(\varepsilon / 2) \frac{\partial U}{\partial q_{i}}(q(t)) \\ q_{i}(t+\varepsilon) &=q_{i}(t)+\varepsilon \frac{p_{i}(t+\varepsilon / 2)}{m_{i}} \\ p_{i}(t+\varepsilon) &=p_{i}(t+\varepsilon / 2)-(\varepsilon / 2) \frac{\partial U}{\partial q_{i}}(q(t+\varepsilon)) \end{aligned} \]

它先更新半個時間步長獲取動量新值，然后用動量的新值對位置變量進行一個完整的時間步長 \(\varepsilon\)的更新，最后再用新的位置變量對動量變量的新值進行另一個半步更新獲得完整的動量更新。這種更新方法每次只更新動能或者位置，但因為\((2.6)\)是剪切變換，因此蛙跳法准確的保留了積分空間體積。

圖 2 $H(q,p)=q^2/2+p^2/2$時，不同的積分算法和參數設置對精度的影響。模擬20個時間步

3. HMC算法

推薦閱讀：

• 博客：HMC采樣可視化，可互動，方便理解。
• 博客：Matlab實現，MH與HMC對比，抽樣。
• MCMC using Hamiltonian dynamics： 主要參考文獻（70%）
• A Conceptual Introduction to Hamiltonian Monte Carlo：主要參考文獻（30%）

3.1 概率和哈密頓量：正則分布

准備好以上知識之后，HMC算法便可以這樣理解：將目標分布的概率密度函數轉化為哈密頓動力學中對應的勢能函數，然后引入動量變量 \(p\)來配合位置變量 \(q\)（感興趣的參數）模擬該變量一定時間長度的運動學軌跡並到達該次模擬的終點（即建議點），然后對哈密頓動力學發現的建議點進行Metropolis-Hasting更新，之后在每個迭代中重新采樣動量，重復執行直至模擬出一個馬爾可夫鏈。
實際上，在任何MCMC方法中我們想要采樣的分布都可以通過統計力學中的吉布斯正則分布（Canonical Distribution）與能量函數聯系起來。給定某個物理系統的某個狀態 \(x\)的能量函數 \(E(x)\)，在統計力學中它的正則分布有概率且概率密度函數為： $$P(x)=\frac{1}{Z}\exp(\frac{-E(x)}{T}) \tag{3.1}$$ 這里的 \(T\)代表系統的溫度，\(Z\)代表歸一化常數。反過來看，如果我們對一些概率密度函數 \(P(x)\)感興趣，可以假設 \(T=1\)，\(E(x)=-\log{P(x)}-\log{Z}\)來得到一個正則分布，其中 \(Z\) 可以是任意值，因為后續的計算中會約去。
HMC中哈密頓量是位置 \(q\)和動量\(p\)聯合的能量函數，因此 $$P(q,p)=\frac{1}{Z}\exp(\frac{-H(q,p)}{T})$$如果 \(H(q,p)=U(q)+K(p)\), 那么聯合概率分布為：

\[P(q,p)=\frac{1}{Z}\exp(\frac{-U(q)}{T})\exp(\frac{-K(q)}{T}) \tag{3.2} \]

並且注意到 \(p\)和\(q\)是相互獨立的且均存在一個能量分別為 \(U(q)\), \(K(p)\)的獨立的正則分布。在貝葉斯統計學中，參數的后驗分布往往是關注的焦點。這些參數扮演着HMC中的位置變量 \(q\)。我們可以將后驗分布的正則分布 \((T=1)\)用勢能函數表達為： $$U(q)=-\log{[L(q|D)\pi(q)]} \tag{3.3}$$這里的 \(\pi(q)\)表示 \(q\)的先驗概率密度函數 \(L(q|D)\)表示在給定數據 \(D\)情況下的似然函數。
至此，我們對HMC有了整體的了解。HMC方法只能用於 \(R^d\)上的連續分布中取樣，這些分布的概率密度函數可以被評估且處處不為零。我們還必須能夠計算該概率密度函數的對數的偏導數。HMC從公式 \((3.2)\)定義的正則分布中抽取 \(q\)和\(p\)，其中 \(q\)由勢能函數\(U(q)\)指定。我們可以選擇動量 \(p\)的分布，正如之前在公式 \((2.3),(2.5)\)所描述一樣。

3.2 HMC的兩個步驟

HMC算法的每個迭代都有兩個步驟：第一步動量會改變，第二步可能同時改變位置和動量。在第一步中，動量變量是從給定的高斯分布中隨機抽取的，與位置變量的當前值無關。對於公式\((2.5)\)中每個動量變量 \(p_i\)它的方差為 \(m_i\)。在第二步中，利用哈密頓動力學模擬並建議一個新的位置變量，然后進行一次Metropolis更新。從當前狀態 \((q,p)\)開始，利用leapfrog方法(或其他可以保持體積的可逆方法)進行 \(L\) 次步長為 \(\varepsilon\)的動力學模擬。這兩個參數需要調節以獲得最好的模擬效果，且較為敏感，一旦設置不合理，不能保證MCMC方法所需要滿足的各種性質將不能獲得正確的結果。在模擬完該次動力學軌跡后，在建議位置將動量變量置負 (保證Metropolis建議的對稱性，實踐中是不需要的，因為二次動能 \(K(p)=K(-p)\))，然后給出該次模擬最終的建議狀態 \((q^*,p^*)\)。這一狀態被作為馬爾可夫鏈的下一個狀態的概率，即接受率為：

\[\min[1,\exp(-H(q^*,p^*)+H(q,p))] = \min[1,\exp(-U(q^*)+U(q)-K(q^*)+K(p))] \]

如果狀態不被接受，那下一個狀態與當前狀態一樣（拒絕的點也會參與蒙特卡羅積分，只是沒有發生轉移，但也是馬氏鏈的一個狀態）。
值得注意的是，第一步的動量重采樣對實施HMC是至關重要的，因為哈密頓量 \(H\)是不變或基本不變的，否則模擬的運動方程將只在同一概率測度的超表面（可以理解為能級）移動，並不能近似模擬出目標分布的典型集。通過不斷的改變每次模擬的動量，就像是在無摩擦光滑地面上給冰球隨便踹一腳，冰球就有機會遍歷該地面的不同勢能的位置。這種解釋不夠嚴謹但很直觀。另外第二步中的Metropolis更新在一定程度上補償了初始位置和建議位置的點因模擬而產生的能量差異，沒有被精確模擬的運動軌跡也有一定的概率接受，這也是HMC算法在超參數合理情況下接受率高的另一個原因。

4. Python簡單實例

在熟悉了HMC的理論之后，我們可以將其用於從已知概率密度復雜分布中抽取樣本。但很多情況下，我們面對的都是模型參數估計問題。在貝葉斯數據分析中，感興趣的往往是參數難以計算的后驗分布，應該這里給出一個簡單的例子來說明如何實現HMC並估計模型中的未知參數。
已知模型為 \(f(x)=ax^{b} + 450\)。 模型參數 \(q=(a,b)\), 其中 \(a=10\), \(b=0.2\)。假設我們現在只觀測到區間 \(x\in [0,4]\)區間上的若干帶噪音的數據 \(\mathcal{D}=(\boldsymbol{x_{obs}},\boldsymbol{y_{obs}})\)（圖3）, 並不清楚該非線性模型中參數的具體數據，希望通過這些觀測值去推理模型參數的大致的取值分布。

圖 3 觀測數據與真實數據

首先，我們需要構建公式\((3.3)\)所描繪的勢能函數將需要解決的問題與已知數據連接起來。觀測值與模型輸出之間的殘差為: $$\boldsymbol{\varepsilon}=\boldsymbol{y_{obs}}-f(\boldsymbol{x_{obs}},q)$$ 假設殘差\(\varepsilon\)為同方差高斯噪音,即 \(\boldsymbol{\varepsilon} \sim \mathcal{N}(0,I)\)，因為各觀測值相互獨立，模型參數 \(q\)的似然函數便可以寫為：

\[L(q|\boldsymbol{\varepsilon})=P(\boldsymbol{\varepsilon}|q)=\Pi^{n}_{i=0}P(\varepsilon_i|q)=\Pi^n_{i=0}\frac{1}{\sqrt{2\pi}}\exp(-\frac{\varepsilon_i^2}{2}) \]

兩邊取負對數：

\[U(q)=-\log{L(q|\boldsymbol{\varepsilon})}=\sum^n_{i=0}\frac{1}{2}\log{2\pi}+\frac{1}{2}\varepsilon_i^2 \]

因為這里使用的是極大似然估計，沒有關於參數的先驗信息（無信息先驗），該式即為本算例的勢能函數。利用構造的勢能函數帶入HMC流程中，模擬出目標參數的典型集即可求解模型的參數分布。圖4，圖5給出了最后的估計結果。算例代碼

圖 4 HMC采樣的結果。當積分步長和積分次數選擇得當，該系統運動方程將被精確模擬，因此接受率很高，且模擬的馬氏鏈可以很快接近並很好的跨越參數聯合分布的典型集。從圖中可以看出各狀態之間的相關性是比較低的，相鄰兩個狀態並不是只停留在很小的鄰域之內

圖 5 參數a,b的邊緣分布。因為HMC更好的模擬的參數的典型集樣本，所以對於參數的期望和方差分析也更為合理

總結

HMC算法毫無疑問是優秀的，它給高維復雜分布的蒙特卡洛積分提供了優秀的解決方案。但是HMC對於leapfrog積分方法中的\(L\)和\(\varepsilon\)太過敏感，這對於HMC在很多實際應用中的實施面臨困難，同時HMC也有可能面對陷入“局部”地區的困境。HMC依賴於勢能函數的一階導數，多次的積分運算模擬使得在單次蒙特卡洛模擬中耗時較多。但相比於MH采樣，整體是還是大大提升了MCMC方法在較復雜問題上的效率和精度，接受率大大提高，穿越並停留在目標分布的典型集的能力大大提升，因此也是優秀的不確定性量化方法。學者們還不遺余力的發展了諸如NUTS, SGHMC等一系列改進的算法。其中NUTS算法憑借其自適應的調節\(L\)和\(\varepsilon\)的能力和有效性在多個庫中作為默認采樣算法。