ハミルトニアンモンテカルロ法

ハミルトニアン モンテカルロ法の導入の理由

HMC法(ハミルトニアンモンテカルロ法)はMH法(メトロポリスヘイスティング法)と比べて圧倒的に受容率の高いサンプリング手法です。
低次元であればMH法でも使えますが、10次元ぐらいになると、受容率の低さ故使い物になりません。
そこで、HMC法を導入します。
ハミルトニアンは物理学におけるエネルギーに対応する物理量です。
解析力学、量子力学で利用されますが、本記事では解析力学における性質を説明します。

ハミルトニアン

※以下で、ハミルトニアンについて説明しますが、厳密ではありませんので、詳しくは解析力学の書籍を参考にしてください。

坂を転げ落ちる物体を例にとります。

f:id:olj611:20210325045938p:plain:w500

時間を $\tau$ として、位置 $\theta(\tau)$ 、速度 $v(\tau)$ 、加速度 $a(\tau)$ は以下の関係式が成り立ちます。

$\begin{eqnarray} &&v(\tau)=\frac{d}{d\tau}\theta(\tau)\tag{1}\\ &&a(\tau)=\frac{d}{d\tau}v(\tau)=\frac{d^2}{d\tau^2}\theta(\tau)\tag{2}\\ \end{eqnarray}$

図の赤い矢印は速度を表しています。

物体の質量を $m$ として、運動量 $p(\tau)$ と力 $F(\tau)$ は以下のようになります。

$\begin{eqnarray} &&p(\tau)=mv(\tau)\tag{3}\\ &&F(\tau)=\frac{d}{d\tau}p(\tau)=\frac{d}{d\tau}mv(\tau)=ma(\tau)\tag{4}\\ \end{eqnarray}$

位置エネルギーはポテンシャルエネルギーとも呼ばれます。
ポテンシャルエネルギーは坂の一番上にいる時が最大で、坂の一番下にいる時が最小です。
ポテンシャルエネルギー $U(\tau)$ は力×移動距離で表されるので、以下のようになります。

$\begin{eqnarray} U(\tau)=F(\tau)\times h(\tau)=ma(\tau)h(\tau)=mgh(\tau)\tag{5} \end{eqnarray}$

ここで、 $g$ は重力加速度と呼ばれる定数なので、 $(5)$ の式変形で $g=a(\tau)$ としました。

物体の移動距離 $は\frac{1}{2}a(\tau)\tau^2$ で表されます。
運動エネルギー $K(\tau)$ も力×移動距離で表されるので、以下のようになります。

$\begin{eqnarray} K(\tau)=F(\tau)\times \frac{1}{2}a(\tau)\tau^2=ma(\tau)\frac{1}{2}a(\tau)\tau^2=\frac{1}{2m}p(\tau)^2\tag{6} \end{eqnarray}$

ポテンシャルエネルギー $P(\tau)$ と運動エネルギー $K(\tau)$ の和をハミルトニアン $H(\tau)$ といいます。

$\begin{eqnarray} H(\tau)=U(\tau)+K(\tau)\tag{7} \end{eqnarray}$

摩擦や空気抵抗などがない理想状態では、物体はハミルトニアンが一定になるように運動します。

あとで、位置 $\theta(\tau)$ を母数空間に対応させ、高さ $h(\tau)$ を負の対数確率に対応させる予定ですので、
自然界とは異なり、以下のように高さを位置で表現することにします。

$\begin{eqnarray} h(\tau)=h(\theta(\tau))\tag{8} \end{eqnarray}$

これにより、ポテンシャルエネルギー $U(\tau)$ は以下のようになります。

$\begin{eqnarray} U(\tau)=mgh(\theta(\tau))\tag{9} \end{eqnarray}$

このとき、ハミルトニアンは以下のようになります。

$\begin{eqnarray} H(\tau)&=&U(\tau)+K(\tau)\\ &=&mgh(\theta(\tau))+\frac{1}{2m}p(\tau)^2\tag{10} \end{eqnarray}$

物体の運動は、ハミルトンの正準方程式と呼ばれる以下の微分方程式によって定められます。

$\begin{eqnarray} &&\frac{d}{d\tau}p(\tau)=-\frac{d}{d\theta(\tau)}U(\theta(\tau))\tag{11}\\ &&\frac{d}{d\tau}\theta(\tau)=\frac{d}{dp(\tau)}K(p(\tau))\tag{12}\\ \end{eqnarray}$

このハミルトンの正準方程式 $(11),(12)$ が成り立つとき、ハミルトニアンは時間変化によらず、一定になります。
このことを確認してみます。

$\begin{eqnarray} \frac{d}{d\tau}H(\tau)&=&\frac{d}{d\tau}U(\tau)+\frac{d}{d\tau}K(\tau)\\ &=&\frac{d}{d\theta(\tau)}U(\theta(\tau))\frac{d}{d\tau}\theta(\tau)+\frac{d}{dp(\tau)}K(p(\tau))\frac{d}{d\tau}p(\tau)\\ &=&-\frac{d}{d\tau}p(\tau)\frac{d}{d\tau}\theta(\tau)+\frac{d}{d\tau}\theta(\tau)\frac{d}{d\tau}p(\tau)\\ &=&0\tag{13} \end{eqnarray}$

(13)の2行目から3行目の変形では $(11),(12)$ を代入しました。
(13)より、ハミルトンの正準方程式が成り立つとき、ハミルトニアンは時間変化によらず、一定になることが確認できました。

ここで、 $m=1,g=1$ として、ハミルトンの正準方程式を書き直します。

$\begin{eqnarray} &&\frac{d}{d\tau}p(\tau)=\frac{d}{d\theta(\tau)}h(\theta(\tau))=-h'(\theta(\tau))\tag{14}\\ &&\frac{d}{d\tau}\theta(\tau)=p(\tau)\tag{15}\\ \end{eqnarray}$

$(14),(15)$ を解析的に解くのではなく、繰り返して解くことを考えます。
オイラー法という手法は以下の方法です。

$\begin{eqnarray} &&p(\tau+1)=p(\tau)-\epsilon h'(\theta(\tau))\tag{16}\\ &&\theta(\tau+1)=\theta(\tau)-\epsilon p(\tau)\tag{17}\\ \end{eqnarray}$

このオイラー法は、ハミルトニアンの保存の精度が良くないことが知られているので
リープフロッグ法と呼ばれる手法が用いられます。

$\begin{eqnarray} &&p(\tau+\frac{1}{2})=p(\tau)-\frac{\epsilon}{2}h'(\theta(\tau))\tag{18}\\ &&\theta(\tau+1)=\theta(\tau)+\epsilon p(\tau+\frac{1}{2})\tag{19}\\ &&p(\tau+1)=p(\tau+\frac{1}{2})-\frac{\epsilon}{2}h'(\theta(\tau+1))\tag{20}\\ \end{eqnarray}$

リープフロッグ法は半単位時間離れた運動量 $p(\tau+\frac{1}{2})$ を使って更新する手法です。

ここまで、 $\tau$ は物理世界の実際の時間を表していました。
しかし、HMC法を考える時は確率密度関数内を1回遷移する際の時間になるため、 $\tau$ を消します。

すると、ポテンシャルエネルギーと運動エネルギーは以下のようになります。

$\begin{eqnarray} U(\theta)=h(\theta)\tag{21}\\ K(p)=\frac{1}{2}p^2\tag{22}\\ \end{eqnarray}$

ポテンシャルエネルギー $U$ が位置 $\theta$ の関数、運動エネルギー $K$ が運動量 $p$ の関数になっていることが分かります。
この時、ハミルトニアンは以下のように、位置 $\theta$ と運動量 $p$ の関数になります。

$\begin{eqnarray} H(\theta,p)=U(\theta)+K(p)=h(\theta)+\frac{1}{2}p^2\tag{23} \end{eqnarray}$

位置と運動量からなる空間を位相空間と言います。
位相空間には、ハミルトニアンの等高線を描くことができます。

f:id:olj611:20210325183117p:plain

青色の点はリープフロッグ法で更新した $(p,\theta)$ です。等高線に沿って動いていることが図から分かります。

位相空間には可逆と体積保存という2つの性質があります。
可逆とは、運動している物体を任意の時間で止め、逆向きの同じ強さの運動量を与えると今来た経路を正確に逆戻りする性質です。
体積保存とは、位相空間内の流れは、すべての点で流れ込む量と流れ出る量が釣り合っているという性質です。

HMC法

$\bf x$ はデータの集合、 $\theta$ はパラメータとします。
我々が欲しいのは、事後分布 $f(\theta|{\bf x})$ に従うサンプルです。

事後分布 $f(\theta|{\bf x})$ とそれとは独立な標準正規分布 $f(p)$ との同時分布 $f(\theta,p|{\bf x})$ を考えます。

$\begin{eqnarray} f(\theta,p|{\bf x})=f(\theta|{\bf x})f(p)\tag{24} \end{eqnarray}$

$f(\theta|{\bf x})$ と $f(p)$ は独立なので、 $f(\theta,p|{\bf x})$ からの乱数は $f(\theta|{\bf x})$ からの乱数と同じです。

$(24)$ を式変形します。

$\begin{eqnarray} f(\theta,p|{\bf x})&=&\exp(\ln(f(\theta,p|{\bf x})))\\ &=&\exp(\ln f(\theta|{\bf x})+\ln f(p))\\ &=&\exp\left(\ln f(\theta|{\bf x})+\ln \left(\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{1}{2}p^2\right)\right)\right)\\ &\propto&\exp\left(\ln f(\theta|{\bf x})+\ln \exp\left(-\frac{1}{2}p^2\right)\right)\\ &=&\exp\left(-h(\theta)-\frac{1}{2}p^2\right)\\ &=&\exp(-H(\theta,p))\tag{25}\\ \end{eqnarray}$

なんと、ハミルトニアンが出てきました！
$(25)$ の4行目から5行目の式変形で対数事後分布 $\ln f(\theta|{\bf x})$ を $-h(\theta)$ とおきました。
ということは、 $\theta$ を物理における位置とみなしたということです。

乱数発生のイメージは次のようなものです。
物体に標準正規分布に従う初速度を与え、一定時間 $L$ の経過後に止め、その位置を記憶します。
その位置から再び物体に標準正規分布に従う初速度を与え、一定時間 $L$ の経過後に止め、その位置を記憶します。
これを繰り返します。

初速度を与えた瞬間の物体の状態を $(\theta,p)$ とし、時間 $L$ 後の遷移先での状態を $(\theta',p')$ とします。
この時、遷移確率は、位相空間の可逆と体積保存という性質により

$\begin{eqnarray} f(\theta',p'|\theta,p)=f(\theta,p|\theta',p')\tag{26} \end{eqnarray}$

と等しくなります。
$(26)$ は $(\theta,p)$ から $(\theta',p')$ への遷移確率と $(\theta',p')$ から $(\theta,p)$ への遷移確率が等しいことを表しています。(下図参照)

f:id:olj611:20210325183139p:plain

また、ハミルトニアンは変化しないので

$\begin{eqnarray} f(\theta,p|{\bf x})=f(\theta',p'|{\bf x})\tag{27} \end{eqnarray}$

となります。

$(26),(27)$ より、

$\begin{eqnarray} f(\theta',p'|\theta,p)f(\theta',p'|{\bf x})=f(\theta,p|\theta',p')f(\theta,p|{\bf x})\tag{28} \end{eqnarray}$

であり、詳細釣り合い条件が成り立っています。

リープフロッグ法では、 $\epsilon$ を使って物体の移動をシミュレートしているので
数値計算上の誤差が生じます。これをMH法のように補正係数を使い、確率的補正します。

$\begin{eqnarray} r=\frac{f(\theta^{(a)},p^{(a)}|{\bf x})}{f(\theta^{(t)},p^{(t)}|{\bf x})}=\exp(H(\theta^{(t)},p^{(t)})-H(\theta^{(a)},p^{(a)}))\tag{29} \end{eqnarray}$

$(29)$ の $\theta^{(t)},p^{(t)}$ は $t$ 回目の点であり、 $\theta^{(a)},p^{(a)}$ はリープフロッグ法により遷移した $t+1$ 回目の候補点です。

HMC法アルゴリズム

1: 初期値 $\theta^{(1)},\epsilon,L,T,$ バーンイン期間を求めます。 $t\leftarrow1$ とします。
2: 独立な標準正規乱数 $p^{(t)}$ を発生します。
3: リープフロッグ法で遷移させ、候補点 $\theta^{(a)},p^{(a)}$ を入手します。
4: 確率 $\min(1,r)$ で受容 $\left(\theta^{(t+1)}=\theta^{(a)}\right)$ し、そうでなければその場に留まります $\left(\theta^{(t+1)}=\theta^{(t)}\right)$ 。
5: $T=t$ なら、終了します。
6: $t\leftarrow t+1$ として、2へ戻ります。