何となくHamiltonian Monte Carloの挙動を実際に見てみたかったので実装してみた
Hamiltonian Monte Carloについては
arxiv.org
がわかりやすかった*1
Markov Chain Monte Carlo
Hamiltonian Monte Carlo、略してHMCはMarkov Chain Monte Carlo(MCMC)に分子動力学法の知識を援用したものだけど、そもそもMCMCとは何かについて軽く触れておく
MCMCの前2ワードのMarkov chainとはインデックス集合が自然数であるMarkov processのことなのだけど、厳密さを抜きに話すと確率過程のうち各についての分布がに従うようなものを指す
このMarkov chainの各の分布をと書くことにすると、は適当な条件のもとで定常分布に収束することが知られている
そこで、その条件を満たしつつ所望の確率分布に収束するようなMarkov chainをうまく構成することができれば、そのMarkov chainを十分なステップ更新することで近似的にからサンプリングを行うことができ、この値を用いてMonte Carloを行う手法をMarkov Chain Monte Carloと呼ぶ
これだけでは何だかサンプリングするだけなのになぜこんな遠回りを、、、と思うかもしれないが、MCMCにはベイズ推定などと親和性が高いアルゴリズム的な理由がある
ベイズ推定では確率モデルと事前分布に対してが観測されたとすると、事後分布は
p(w) &= \frac{1}{Z_n(\beta)} \phi(w) \prod_{i=1}^n p(X_i | w)^\beta = \frac{1}{Z_n(\beta)} \exp(-\beta H(w))
\end{align*}
ベイズ推定では一般的に、正規化項を解析的に求めることが難しく、それゆえ事後分布から直接サンプリングを試みることは出来ないわけだが、MCMCでは正規化項が計算できなくとも、さえ計算できれば事後分布から近似的にサンプリングを行うことができる
Metropolis Algorithm
上ではベイズ推定を例に出したが、より一般にの形をしている分布に対し、正規化項を計算せずにサンプリングを行えることをMCMCの中でも一番有名なMetropolis法を用いて見ていく
適当な初期状態からはじめ、ステップ後にが得られている時、まず提案分布によって次の遷移先の候補をサンプリングする
ここでは基本的には何でも良いが、対称的、すなわちである必要がある*3
こうして得られた遷移先の候補に対して、
P := \min \{1, \exp(-\beta \Delta H)\} \quad (\text{ただし}\Delta H := H(w') - H(w))
\end{align*}
上で適当な条件のもとMarkov chainはに収束すると述べたが、その具体的な条件とは雑に言って
が満たされれば良く、MCMCでは特にが定常分布であることの十分条件である詳細釣り合い条件
p(w_b | w_a)\pi(w_a) = p(w_a | w_b)\pi(w_b)
\end{align*}
下のgifはMetropolis法を用いて混合正規分布
&0.3 \times \mathcal{N}((-1,-1)^\top, 1_2) + 0.7 \times \mathcal{N}((2,2)^\top, 1_2) \\
&= 0.3 \times \frac{1}{2\pi}e^{-\frac{(x_1 + 1)^2 + (x_2 + 1)^2}{2}} + 0.7 \times \frac{1}{2\pi}e^{-\frac{(x_1 - 2)^2 + (x_2 - 2)^2}{2}}
\end{align*}
大きい丸は遷移先の候補で、それが採用されると小さい丸になるようにプロットしている
このgifを見ると、結構な割合で候補が棄却されており何となく効率が悪そうに見えるのと、混合正規分布の2つの山の谷にあたる部分を飛び越えるのに苦労している、すなわち各ステップの相関が高いことがわかる
Hamiltonian Monte Carlo
HMCもMetropolis法の一種であり、同様の手続きを踏むことになるが、名前の通りHamiltonianの概念が登場する
HMCでは座標と運動量に対して
\mathcal{H}(w,p) = \frac{1}{2}p^\top p + \beta H(w)
\end{align*}
1. 状態から次の遷移を行う場合、まずを新たに標準正規分布からサンプリングし置き換える
2. 1で得られたを正準方程式に従って秒後まで更新する:
\frac{dw_i}{dt} &= \frac{\partial H}{\partial p_i} \\
\frac{dp_i}{dt} &= -\frac{\partial H}{\partial w_i}\end{align*}
3. 2で得られた遷移先の候補を、確率
P := \min\{1, \exp(-\Delta \mathcal{H})\} \quad (\text{ただし}\Delta \mathcal{H} := \mathcal{H}(w^*, p^*) - \mathcal{H}(w, p))
\end{align*}
この手続きを繰り返して得られるはに従い、の形からとは独立であるから、特にはに従う*7
HMCがを定常分布に持つことは詳細釣り合い条件が満たされることを確かめることによって示されるが、Metropolis法のように直接詳細釣り合い条件を確認するわけではない
というのも、上のステップ1は運動量のみを更新するため定常分布には影響を与えず、ゆえに2,3の手続きが詳細釣り合い条件を満たすことを確認すれば良く、これは簡単な計算によって示される
Leapfrog Method
上の正準方程式の数値解法としては、leapfrog法と呼ばれる手法を用いることが推奨されている:
p_i(t + \epsilon/2) &= p_i(t) - \frac{\epsilon \beta}{2} \frac{\partial H}{\partial w_i}(w(t)) \\
w_i(t + \epsilon) &= w_i(t) + \epsilon p_i(t + \epsilon/2) \\
p_i(t + \epsilon) &= p_i(t + \epsilon/2) - \frac{\epsilon \beta}{2} \frac{\partial H}{\partial w_i}(w(t + \epsilon))\end{align*}
Langevin Equation
HMCでは本来、正準方程式を時刻に依存しない、時間分だけ進めるわけだが、ここでleapfrog法を各ステップごとに1イテレーションだけ実行した場合、Langevin Monte Carlo (LMC)と呼ばれる、これまた別のMCMCと等価になることが知られている
\frac{d w}{dt} = -\beta \nabla H(w) + \frac{d B}{dt}
\end{align*}
w_i(t+\epsilon) = w_i(t) - \frac{\epsilon^2\beta}{2}\frac{\partial H}{\partial w_i}(w(t)) + \epsilon p_i(t)\end{align*}
となり、が正規分布からサンプリングされることから上の確率微分方程式のEuler-Maruyama近似になっていることがわかる*10
Langevin方程式は粘性が大きく加速度による影響を無視できるとした上での、粒子がポテンシャルとホワイトノイズとして定式化される他の粒子との衝突を受ける時の運動を記述したものであり、少なくとも連続的な確率微分方程式の場合はポテンシャルが適当な正則性を満たせば*11でに収束することが知られている
下のgifはMetropolis法と同じ混合正規分布を、leapfrog法を1イテレーションのみ更新するHMCでサンプリングした様子
ノイズが乗った勾配降下法のような挙動をしており、また、leapfrog法はHamiltonianを保存するため棄却率が大幅に減ったためMetropolis法よりは効率が良さそうだが、やはりランダムウォークのような動きをするため点同士の相関が高い
最後に、各ステップにleapfrog法を10イテレーション回した場合の様子*12
推進力、または突破力のようなものが加わり、各ステップの点の位置が大きく離れるようになった
*1:そもそもはHybrid Monte Carloという名前で提案された手法だけど、ここで引用した論文の著者はあえてHamiltonian Monte Carloと呼んでいる
*2:適当なから始めて、各ステップの確率で、の確率でに移動するようなもの
*3:何でも良いと言ってもMarkov chainが収束するための条件は満たしている必要があり、また対称的でなくても遷移候補の棄却確率を少し煩雑にすることで対応できる
*4:がに従うときもまたに従うこと
*5:が任意の状態について成り立つこと
*6:本来は2つ目と3つ目の条件も確認するべきなんだろうけれど、基本的に無視されていると思われる
*7:は正規化項
*8:なので天体の軌道計算などにも使われるみたいな話を何かの授業で聞いたことある気がする
*10:正確にはを新しいステップ幅にとる必要があるので、時間軸がスケール変換される気がする
*11:例えばポテンシャルの裾が重すぎると粒子が遠くに行きがちになり収束しない
*12:このgifではleapfrog法による途中経過もプロットしている