趣旨

ハミルトニアンモンテカルロ法(HMC)を理解する(執筆当時ではまだ理解できてないけど)。

注意

物理ど素人なので完全には理解できてない。情報の正確性には注意。

理解してないなら「理解する」とか書くなよ

うるせーよ

では解説する。

HMCとは

Metropolis-Hastingsに代表されるこれまでのMCMCでは、あらかじめ与えられた初期値と初期提案分布からスタートし、棄却したりしなかったりしながらサンプリングを行うわけだが、提案分布が目的の分布と大きく外れていれば酷い出来になる。これを軽減するためにランダムウォークMetropolis-Hastingsがあるらしいのだが、あんまり探索的な(提案分布から推察すれば生起確率が低そうな)サンプリングを行うと高い確率で棄却されてしまう。要するに、HMCは強化学習における「探索と活用のトレードオフ」に似た問題を抱えている。これをハミルトン力学の観点から解決したのがHMC、らしい。つまり、ある程度探索的なサンプリングを確保しながら、同時に採択率を高く保てる方法がHMCなのである。

1. 初等物理量

ハミルトン「力学」とあるとおり、HMCの起源は物理学、より厳密には解析力学にある。よって、HMCの理解にはある程度の物理的知識が必要であるため、一旦物理のおさらいをしたい。ここで、時間を $\tau$ で表現することにする。

1-1. 位置、速度、加速度

空間上を動く物体を考える。この物体の時刻 $\tau$ における位置をベクトル $\boldsymbol{r} (\tau)$ で表現する。また、同時刻のおける物体の速度 $\boldsymbol{v} (\tau)$ は、速度が瞬間的な位置の変化を表現していることに注意すれば、

$\begin{align} \boldsymbol{v} (\tau) = \frac{d \boldsymbol{r} (\tau)}{d \tau} \end{align}$

と表現できるだろう。また、加速度 $\boldsymbol{a} (\tau)$ も同様に考えれば、

$\begin{align} \boldsymbol{a} (\tau) = \frac{d \boldsymbol{v} (\tau)}{d \tau} = \frac{d^2 \boldsymbol{r} (\tau)}{d \tau^2} \end{align}$

と書ける。

1-2. 質量、運動量、力

物体の質量を $m$ としたとき、時刻 $\tau$ における運動量 $\boldsymbol{p} (\tau)$ を

$\begin{align} \boldsymbol{p} (\tau) = m \boldsymbol{v} (\tau) \end{align}$

と定義する。これは誤解を恐れずに言うなら「物体の止めにくさ」を意味する。

例え高速で飛来する物体でも、それが野球ボールなら簡単に止められるだろう。一方で、ゆっくり動いていても、それが船ならまず人力では止められない。この辺りの違いを数値的に表現したのが運動量である。

また、運動量の変化を力 $\boldsymbol{F} (\tau)$ という。

$\begin{align} \boldsymbol{F} (\tau) = \frac{d \boldsymbol{p} (\tau)}{d \tau} \end{align}$

この表現と上の定義から、運動方程式と呼ばれる式が成立する。

$\begin{align} \boldsymbol{F} (\tau) = m \frac{d^2 \boldsymbol{r} (\tau)}{d \tau^2} = m \boldsymbol{a} (\tau) \end{align}$

この式は高校物理でもよく見るだろう。速度が光速に近いと成り立たないらしいが、そんなことは知らん。

1-3. 仕事

仕事とは、物体をある位置からある位置に動かすのに要するエネルギーのことである。具体的には仕事 $W$ は力の積分によって表される: 物体の移動経路を $C$ で表現すると、

$\begin{align} W = \int_C \boldsymbol{F} d \boldsymbol{r} \end{align}$

と表現される。ここで、積分は $\boldsymbol{r}$ について行われることに注意: 直感的には時間 $\tau$ を用いても良いような気がするが、仕事に置いて時間はさほど重要ではない。例えどれだけ時間をかけたとしても、結果として物体がほとんど動かなかったら、仕事としては小さい。極端な話、動かなかったら仕事は $0$ である。シビアやね。

上の式にあるように、仕事は経路 $C$ に依存するため、始点と終点以外の部分も考慮される: 両端点が同じでも、最短経路で動かした場合と、大回りしてたどり着いた場合では仕事の値が違う可能性がある。ここで、仕事が経路 $C$ に依存しない時、力 $\boldsymbol{F} (\tau)$ を保存力と呼ぶことにしよう。 $\boldsymbol{F} (\tau)$ が保存力であるという仮定のもとでは、先ほどの式は端点 $\boldsymbol{r}_s, \boldsymbol{r}_g$ のみを考慮した積分の式

$\begin{align} W = \int_{\boldsymbol{r}_s}^{\boldsymbol{r}_g} \boldsymbol{F} d \boldsymbol{r} \end{align}$

で表される。一般に仕事を評価する時、それは地点 $\boldsymbol{r} (0)$ から地点 $\boldsymbol{r} (T)$ まで、つまり時刻 $0$ から $T$ までの仕事を評価することがよくある。その場合は

$\begin{align} W(0, T) &= \int_{\boldsymbol{r} (0)}^{\boldsymbol{r} (T)} \boldsymbol{F} d \boldsymbol{r} \\\ &= \int_0^T \boldsymbol{F} (\tau) \frac{d \boldsymbol{r} (\tau)}{d \tau} d \tau \\\ &= \int_0^T \boldsymbol{F} (\tau) \boldsymbol{v} (\tau) d \tau \\\ \end{align}$

のように書ける。

2. 力学的エネルギー

運動する物体にはエネルギーを持っている。物体が存在する空間には摩擦や空気抵抗などが存在しないと仮定すると、物体が持つエネルギーはポテンシャルエネルギーと運動エネルギーの2種類に大別される。

2-1. ポテンシャルエネルギー

ポテンシャルエネルギーとは、物体が特定の位置に存在することによって潜在的に持つエネルギーで、位置エネルギーとも言う。保存力 $\boldsymbol{F}$ に対して、基準点 $\boldsymbol{r}_0$ を定めれば、以下のようなポテンシャルエネルギーを定義できる。

$\begin{align} U(\boldsymbol{r}) = -\int_{\boldsymbol{r}_0}^{\boldsymbol{r}} \boldsymbol{F} d \boldsymbol{r}\\\ \end{align}$

ここでもやはり、基準点を $\boldsymbol{r} (0)$ 、そして現在の位置を $\boldsymbol{r} (T)$ としてみると、

$\begin{align} U(\boldsymbol{r} (T)) = -\int_{\boldsymbol{r} (0)}^{\boldsymbol{r} (T)} \boldsymbol{F} d \boldsymbol{r} = -W(0, T) \\\ \end{align}$

と書ける。

例を挙げよう。今、鉛直上向きに $x$ 軸が伸びる座標系を考え、またこの系では鉛直下向きに質量に比例した力 $mg$ が与えられるものとする(要するに、高校力学と同じ環境)。この時、基準点を $0$ としたときの位置 $h$ にある物体のポテンシャルエネルギーは、力の向きに注意すると、

$\begin{align} U(h) = -\int_0^h (-mg) d x = mgh\\\ \end{align}$

と書ける。お馴染みの位置エネルギーである。

別の例として、 $x$ 軸だけが伸びる座標系において、基準点を $0$ としたとき位置 $h$ にある物体に力 $kh$ が与えられるものとする。この場合、

$\begin{align} U(h) = -\int_0^h kx d x = \frac{1}{2} k h^2\\\ \end{align}$

のようにポテンシャルエネルギーが与えられ、これはばねの弾性エネルギーと等価である。ここでいうポテンシャルエネルギーとは、高校物理における位置エネルギーより広い範囲のエネルギーであると言える。

2-2. 運動エネルギー

運動エネルギーは、以下の式で与えられるエネルギーである。

$\begin{align} K(\boldsymbol{p}) = \frac{1}{2} m {||\boldsymbol{v}||}^2 = \frac{1}{2m} {||\boldsymbol{p}||}^2 \\\ \end{align}$

さて、この運動エネルギーについての考察を得るため、両辺を $\tau$ で微分してみる。

$\begin{align} \frac{d K(\boldsymbol{p})}{d \tau} = m \boldsymbol{v} \cdot \frac{d \boldsymbol{v}}{d \tau} = m \boldsymbol{a} \cdot \frac{d \boldsymbol{r}}{d \tau} = \boldsymbol{F} (\tau) \boldsymbol{v} (\tau)\\\ \end{align}$

さらにこれを区間 $(0, T)$ で $\tau$ について積分してみると、

$\begin{align} K(\boldsymbol{p} (T)) - K(\boldsymbol{p} (0)) = \int_0^T \boldsymbol{F} (\tau) \boldsymbol{v} (\tau) d \tau = W(0, T) \\\ \end{align}$

が得られる。

2-3 エネルギー保存則

$\boldsymbol{F}$ が保存力である場合のポテンシャルエネルギー $U(\boldsymbol{r} (T))$ について、

$\begin{align} U(\boldsymbol{r} (T)) = -\int_{\boldsymbol{r} (0)}^{\boldsymbol{r} (T)} \boldsymbol{F} d \boldsymbol{r} = -W(0, T) \\\ \end{align}$

と書けることは上で確認した。ここで、この記法に従えば、 $U(\boldsymbol{r} (0)) = 0$ であることは容易に確認できるため、あえて

$\begin{align} U(\boldsymbol{r} (T)) - U(\boldsymbol{r} (0)) = -W(0, T) \\\ \end{align}$

と記載しておく。さて運動エネルギーの方はというと、

$\begin{align} K(\boldsymbol{p} (T)) - K(\boldsymbol{p} (0)) = W(0, T) \\\ \end{align}$

が成立しているのだった。2つの式の右辺を見てやると、

$\begin{align} U(\boldsymbol{r} (T)) - U(\boldsymbol{r} (0)) = -\lbrace K(\boldsymbol{p} (T)) - K(\boldsymbol{p} (0)) \rbrace \\\ \end{align}$

$\begin{align} \therefore U(\boldsymbol{r} (T)) + K(\boldsymbol{p} (T)) = U(\boldsymbol{r} (0)) + K(\boldsymbol{p} (0)) \\\ \end{align}$

が成立することがわかる。これをエネルギー保存則という。この法則は、ポテンシャルエネルギーと運動エネルギーの和は時間推移によって変化しないことを主張している。

2-4. ハミルトニアン、運動方程式

ここで、ハミルトニアン $H(\tau)$ を以下のように定義する。

$\begin{align} H(\tau) = U(\boldsymbol{r} (\tau)) + K(\boldsymbol{p} (\tau)) \end{align}$

このとき、 $H(\tau)$ は時間変化によらずに一定であると言える。

ここで、「物体」は空間上の曲面を転がっているとする:「物体」が存在する空間は位置に関わらず質量と加速度 $m, g$ に比例した力 $mg$ が与えられ、かつ「物体」の基準点からの高さ $h(\tau)$ は物体の(超)平面上の位置 $\boldsymbol{\theta} (\tau)$ に依存するものとする。この時のポテンシャルエネルギーは、

$\begin{align} U(h(\tau)) &= -\int_0^{h(\tau)} (-mg) d x = h(\tau) = mgh(\boldsymbol{\theta} (\tau)) \\\ \end{align}$

と書ける。ここで、 $U(h(\tau)) = U(\boldsymbol{\theta} (\tau))$ と表現しよう。さて、この時のハミルトニアンについて何が言えるだろうか。

曲面を転がる物体は、時間の変化によらずに、常にハミルトニアンが一定である。すなわち、

$\begin{align} \frac{d H(\tau)}{d \tau} &= 0 \\\ \end{align}$

$\begin{align} \frac{d U(\boldsymbol{\theta} (\tau))}{d \tau} + \frac{d K(\boldsymbol{p} (\tau))}{d \tau} &= 0 \end{align}$

これを変形する。

$\begin{align} \frac{d K(\boldsymbol{p} (\tau))}{d \tau} &= - \frac{d U(\boldsymbol{\theta} (\tau))}{d \tau} \\\ \frac{d K(\boldsymbol{p} (\tau))}{d \boldsymbol{p} (\tau)} \frac{d \boldsymbol{p} (\tau)}{d \tau} &= - \frac{d U(\boldsymbol{\theta} (\tau))}{d \boldsymbol{\theta} (\tau)} \frac{d \boldsymbol{\theta} (\tau)}{d \tau} \\\ \end{align}$

で、なんか知らんが(!?)ここから以下のような方程式が言えるらしい(というか、下の方程式を満たすとき、上の方程式を満たすような条件が見つかるのだろう)。

$\begin{align} \frac{d \boldsymbol{p} (\tau)}{d \tau} &= - \frac{d U(\boldsymbol{\theta} (\tau))}{d \boldsymbol{\theta} (\tau)} \\\ \frac{d \boldsymbol{\theta} (\tau)}{d \tau} &= \frac{d K(\boldsymbol{p} (\tau))}{d \boldsymbol{p} (\tau)} \\\ \end{align}$

これをハミルトンの運動方程式と言い、物体の運動を予測する微分方程式である。

2-5. オイラー法とリープフロッグ法

上で定義したハミルトンの運動方程式を用いて、ハミルトニアンの値を変化させないまま運動量 $\boldsymbol{p}$ や位置 $\boldsymbol{\theta}$ を変化させる、つまり物体の移動を表現することを考える。ここで、計算を簡略化するためにいささか恣意的な仮定を置く: 物体の質量と重力加速度について $m=1, g=1$ とする。これにより、

$\begin{align} \frac{d \boldsymbol{p} (\tau)}{d \tau} &= - \frac{d h(\boldsymbol{\theta} (\tau))}{d \boldsymbol{\theta} (\tau)} \\\ \frac{d \boldsymbol{\theta} (\tau)}{d \tau} &= \boldsymbol{p} (\tau) \\\ \end{align}$

のように簡潔に書ける。

さて、当然解析的には解けるはずがないのだが、数値的には解けそうである。オイラー法と呼ばれる方法は、上に対する素直な解法を提供してくれている。

$\begin{align} \boldsymbol{p} (\tau+ \varepsilon) &= p(\tau) + \varepsilon \frac{d \boldsymbol{p} (\tau)}{d \tau} = p(\tau) - \varepsilon \frac{d h(\boldsymbol{\theta} (\tau))}{d \boldsymbol{\theta}} \\\ \boldsymbol{\theta} (\tau+ \varepsilon) &= \boldsymbol{\theta} (\tau) + \varepsilon \frac{d \boldsymbol{\theta} (\tau)}{d \tau} = \boldsymbol{\theta} (\tau) + \varepsilon \boldsymbol{p} (\tau) \\\ \end{align}$

$\cdots$ まあ要するに最急降下法なのだが、あんまり精度が良くないらしい。なので、以下のような更新式によって精度を高めているらしい。

$\begin{align} \boldsymbol{p} (\tau+\frac{ \varepsilon}{2}) &= p(\tau) - \frac{\varepsilon}{2} \frac{d h(\boldsymbol{\theta} (\tau))}{d \boldsymbol{\theta}} \\\ \boldsymbol{\theta} (\tau+ \varepsilon) &= \boldsymbol{\theta} (\tau) + \varepsilon \boldsymbol{p} \left( \tau + \frac{ \varepsilon}{2} \right) \\\ \boldsymbol{p} (\tau+ \varepsilon) &= \boldsymbol{p} \left( \tau+\frac{ \varepsilon}{2} \right) - \frac{\varepsilon}{2} \frac{d h(\boldsymbol{\theta} (\tau+ \varepsilon))}{d \boldsymbol{\theta}} \\\ \end{align}$

これをリープフロッグ法と呼ぶ。おそらく位置変数 $\boldsymbol{\theta}$ と運動量変数 $\boldsymbol{p}$ を時刻 $\varepsilon/2$ だけずらしながら交互に更新する様をカエルに見立てたのだろう。

オイラー法でもリープフロッグ法でも、結局積分を数値積分で近似しているだけなので、誤差が生じうる。これのハンドリングは後述。

2-6. 位相空間

物理的な空間内にある物体の動きを考える際、時間 $\tau$ について考えることは重要であるが、この後統計学へ応用する際においてはさほど重要な概念ではない。この視点のもとでのエネルギーは、 $m=1, g=1$ に注意すると、

$\begin{align} U(\boldsymbol{\theta}) &= h (\boldsymbol{\theta}) \\\ K(\boldsymbol{p}) &= \frac{1}{2} {|| \boldsymbol{p} ||}^2 \end{align}$

のように表現できる。このときのハミルトニアンは

$\begin{align} H(\boldsymbol{\theta}, \boldsymbol{p}) = U(\boldsymbol{\theta}) + K(\boldsymbol{p}) = h (\boldsymbol{\theta}) + \frac{1}{2} {|| \boldsymbol{p} ||}^2 \end{align}$

のように位置と運動量の関数として表現される。

位置と運動量を座標とする空間を位相空間(phase space)と言う。

数学における位相空間(topological space)とはまた違う。いい加減にしろよ。

で、ハミルトニアンを固定した時、 $(\boldsymbol{\theta}, \boldsymbol{p})$ は特定の軌跡(これを特に等高線と呼ぶようだ)を描く。実空間内の物体は、常にこの等高線に沿って位置と運動量を変化させることになる。

さて、位相空間では2つの重要な性質がある。

可逆: 運動している物体に対して、任意の時間で運動を停止させ、これまでに受けてきた力を逆向きに全く同じ大きさで(つまり $-1$ 倍の運動量を)与えると、物体はこれまでたどってきた軌跡を正確に辿る
体積保存: 位相空間上の領域(点の集合)は、時間発展によって(各点が好き勝手動いて)形を変えることはあっても、その面積(体積とよく表記される)は不変である

可逆については直感的だろう。体積保存の性質は、リウビルの定理(Liouville's Theorem)と呼ばれる定理がその成立を保証してくれている。なぜ成立するかは他のサイトに譲るとして、上の2つの性質によって以下のような重要な式が導かれる。

$\begin{align} \text{Pr} (\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*} | \boldsymbol{\theta}, \boldsymbol{p}) = \text{Pr} (\boldsymbol{\theta}, \boldsymbol{p} | \boldsymbol{\theta}^{*}, \boldsymbol{p}^{*}) \end{align}$

これは、地点 $\boldsymbol{\theta}, \boldsymbol{p}$ から地点 $\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*}$ に遷移する確率は、 $\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*}$ から $\boldsymbol{\theta}, \boldsymbol{p}$ に遷移する確率と等しいことを意味している。

可逆性だけではこの式は導けない。なぜなら、物体に運動量 $\boldsymbol{p}$ を与える確率と、 $-\boldsymbol{p}$ を与える確率は等しいとは言えないからだ。リウビルの定理がその辺りを説き伏せてくれるのだろう。どこにも説明がない。マジでキレそう。

3. MCMCに求められる条件

ハミルトニアンの話を一旦脇に置いて、MCMCのおさらいをしよう。今我々が提案するマルコフ連鎖では、状態 $z$ の分布が目的の分布に収束する必要がある。

3-1. エルゴード性

まずは一意の分布に収束することについて考えよう。この収束先の分布のことを定常分布という。より具体的には、以下の数式を満たす分布 $p^* (z)$ のことを指す。

$\begin{align} p^* (z) = \int p(z | z') p^* (z') dz' \end{align}$

特定の状態にハマると抜け出せないようなマルコフ連鎖などは、この条件を満たさない。定常分布を持つための必要十分条件は、以下の3つの性質を全て満たすことである。

既約的: 任意の状態から他の任意の状態へ到達可能
非周期的: 周期性を持たない
正再帰的: 任意の状態は限りなく何度も訪問される

これら3つの性質を満たすとき、そのマルコフ連鎖はエルゴード性を持つという。

3-2. 詳細釣り合い条件

仮にマルコフ連鎖がエルゴード的であったとして、定常分布が目的の分布と全く違っては意味がない。ある確率分布 $p^* (z)$ が定常分布であるための必要十分条件は分かっていないが、十分条件は詳細釣り合い条件として知られる。

$\begin{align} p^* (z) p(z' | z) = p^* (z') p(z | z') \end{align}$

4. HMC

さて、ハミルトン力学の考え方をMCMCに導入することを考える。

4-1. ハミルトニアンの導入

今、目的の分布である事後分布 $\text{Pr}(\boldsymbol{\theta} | \boldsymbol{x})$ と、それとは独立の標準正規分布 $\text{Pr}(\boldsymbol{p})$ が存在するとする。同時分布は

$\begin{align} \text{Pr}(\boldsymbol{\theta}, \boldsymbol{p} | \boldsymbol{x}) = \text{Pr}(\boldsymbol{\theta} | \boldsymbol{x}) \text{Pr}(\boldsymbol{p}) \end{align}$

のように表現される。これの対数を取ると、

$\begin{align} \log \text{Pr}(\boldsymbol{\theta}, \boldsymbol{p} | \boldsymbol{x}) &= \log \text{Pr}(\boldsymbol{\theta} | \boldsymbol{x}) + \log \text{Pr}(\boldsymbol{p}) \\\ \end{align}$

と書ける。右辺の各項について考察を深めよう。

$\cdots$ と言っておいてアレだが、右辺第一項 $\log \text{Pr}(\boldsymbol{\theta} | \boldsymbol{x})$ については仮定を置くだけである。ハミルトン運動方程式を導いた際に、物体は高さが位置に依存した関数 $h(\boldsymbol{\theta})$ で表現されるような曲面上を動くと書いた。ここでは、

$\begin{align} h(\boldsymbol{\theta}) = - \log \text{Pr}(\boldsymbol{\theta} | \boldsymbol{x}) \end{align}$

と置く。

右辺第二項 $\log \text{Pr}(\boldsymbol{p})$ については具体的に計算してみる。 $\text{Pr}(\boldsymbol{p})$ が標準正規分布であったことに注意すれば、

$\begin{align} \log \text{Pr}(\boldsymbol{p}) &= \log \left\lbrace \frac{1}{\sqrt{2 \pi}} \exp \left( -\frac{1}{2} {||\boldsymbol{p}||}^2 \right) \right\rbrace \\\ &= -\frac{1}{2} {||\boldsymbol{p}||}^2 + \text{Const.} \\\ \end{align}$

と書ける。これらを全部代入すると、各項が2章6節で定義したエネルギーと同じ形をしているため、

$\begin{align} \log \text{Pr}(\boldsymbol{\theta}, \boldsymbol{p} | \boldsymbol{x}) &= -h(\boldsymbol{\theta}) -\frac{1}{2} {||\boldsymbol{p}||}^2 + \text{Const.} \\\ &= -U(\boldsymbol{\theta}) -K(\boldsymbol{p}) + \text{Const.} \\\ &= -H(\boldsymbol{\theta}, \boldsymbol{p}) + \text{Const.} \end{align}$

のように、突然ハミルトニアンが出現する。要するに、

$\begin{align} \text{Pr}(\boldsymbol{\theta}, \boldsymbol{p} | \boldsymbol{x}) = \frac{1}{Z_H} \exp (-H(\boldsymbol{\theta}, \boldsymbol{p})) \end{align}$

のように書ける。ここで、定数項の影響は正規化定数 $1 / Z_H$ で吸収した。

4-2. 運動の解釈

曲面上の物体は、重力に従って高さ $h(\boldsymbol{\theta})$ が小さくなる方向に運動することは想像に難くないだろう。前節によれば、高さは

$\begin{align} h(\boldsymbol{\theta}) = - \log \text{Pr}(\boldsymbol{\theta} | \boldsymbol{x}) \end{align}$

のように定義されているので、高さ $h(\boldsymbol{\theta})$ が小さくなるということは、事後分布 $\text{Pr}(\boldsymbol{\theta} | \boldsymbol{x})$ を最大化することに等しい。

ただし実際は物体は放置されるわけではない。前節で、標準正規分布に従って $\boldsymbol{p}$ がサンプリングされることを仮定した。 $\boldsymbol{p}$ とは運動量のことであるから、物体はただ重力に従って曲面上を転げ落ちるのではなく、ランダムな力を常に加えられながら転がることを意味する。

定期的に物体を止めて、ランダムな方向と力で物体を弾くとイメージすれば良い。

このランダムネスによって、有効なサンプリングが期待される。

4-3. 詳細釣り合い条件

ここで提案されているHMCは、要するに以下の処理を交互に行うマルコフ連鎖であると言える。

運動量 $\boldsymbol{p}$ の確率的な更新
2章4節で定義したハミルトン運動方程式に基づいた、力学系の更新

3章2節で、マルコフ連鎖は詳細釣り合い条件を満足することが求められると説明した。ハミルトン運動方程式に従って $\boldsymbol{\theta}, \boldsymbol{p}$ を変化させた場合、詳細釣り合い条件は満たされるのだろうか。ハミルトン力学における記法に則れば、詳細釣り合い条件は以下のように表現される。

$\begin{align} \text{Pr} (\boldsymbol{\theta}, \boldsymbol{p} | \boldsymbol{x}) \text{Pr} (\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*} | \boldsymbol{\theta}, \boldsymbol{p}) = \text{Pr} (\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*} | \boldsymbol{x}) \text{Pr} (\boldsymbol{\theta}, \boldsymbol{p} | \boldsymbol{\theta}^{*}, \boldsymbol{p}^{*}) \end{align}$

2章6節で、位相空間においては以下の性質を満たすと書いたことを思い出して欲しい。

さらにここでハミルトニアンを保存しながら遷移することから、

$\begin{align} \text{Pr}(\boldsymbol{\theta}, \boldsymbol{p} | \boldsymbol{x}) = \text{Pr}(\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*} | \boldsymbol{x}) = \frac{1}{Z_H} \exp (-H(\boldsymbol{\theta}, \boldsymbol{p})) \end{align}$

が成立する。以上の二つを組み合わせると詳細釣り合い条件が成立する。

4-4. 数値誤差と棄却

上の話は理想的な運動をシミュレートできればの話である。2章5節で紹介したようなハミルトンの運動方程式を解く方法は所詮は近似的なものであり、そこでは数値誤差、つまりハミルトニアンの変動が発生する。そのため、Metropolis-Hastingsアルゴリズムで登場した棄却サンプリングを適用する。

usapyoi.hatenablog.com

このアルゴリズムでは、分布 $q_k (\boldsymbol{z}^{(\tau)})$ からサンプリングしてきた $\boldsymbol{z}^{*}$ に対して、以下のような確率をもって採択しているのだった。

$\begin{align} A_k (\boldsymbol{z}^*, \boldsymbol{z}^{(\tau)}) = \min \left\lbrace 1, \frac{\tilde{p}(\boldsymbol{z}^*) q_k (\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)}{\tilde{p}(\boldsymbol{z}^{(\tau)}) q_k (\boldsymbol{z}^*| \boldsymbol{z}^{(\tau)})} \right\rbrace \end{align}$

今回の例では、目的の分布は $\text{Pr}(\boldsymbol{\theta} | \boldsymbol{x})$ であり、状態遷移分布は $\text{Pr} (\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*} | \boldsymbol{\theta}^{(\tau)}, \boldsymbol{p}^{(\tau)})$ のように表現される。後者の分布は対称的であることに注意すると、

$\begin{align} \frac{\tilde{p}(\boldsymbol{z}^*) q_k (\boldsymbol{z}^{(\tau)} | \boldsymbol{z}^*)}{\tilde{p}(\boldsymbol{z}^{(\tau)}) q_k (\boldsymbol{z}^*| \boldsymbol{z}^{(\tau)})} &\rightarrow \frac{\text{Pr}(\boldsymbol{\theta}^* | \boldsymbol{x}) \text{Pr} ( \boldsymbol{\theta}^{(\tau)}, \boldsymbol{p}^{(\tau)} | \boldsymbol{\theta}^{*}, \boldsymbol{p}^{*}) }{\text{Pr}(\boldsymbol{\theta}^{(\tau)} | \boldsymbol{x}) \text{Pr} (\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*} | \boldsymbol{\theta}^{(\tau)}, \boldsymbol{p}^{(\tau)})} \\\ &= \frac{\text{Pr}(\boldsymbol{\theta}^* | \boldsymbol{x}) }{\text{Pr}(\boldsymbol{\theta}^{(\tau)} | \boldsymbol{x})} \\\ &= \frac{ \text{Pr}(\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*} | \boldsymbol{x}) }{ \text{Pr}(\boldsymbol{\theta}^{(\tau)}, \boldsymbol{p}^{(\tau)} | \boldsymbol{x})} \\\ &= \exp \left( H(\boldsymbol{\theta}^{(\tau)}, \boldsymbol{p}^{(\tau)}) - H(\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*}) \right) \end{align}$

と書ける。ただし、途中で $\text{Pr}(\boldsymbol{p})$ を分母分子にかけている。

結局、数値誤差によるズレは、以下のような確率で採択するMHアルゴリズム(より正確にはMetropolisアルゴリズム)でカバーするのであった。

$\begin{align} r^{(\tau)} = \min \left\lbrace 1, \exp \left( H(\boldsymbol{\theta}^{(\tau)}, \boldsymbol{p}^{(\tau)}) - H(\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*}) \right) \right\rbrace \end{align}$

$\exp$ の中身は実際はほぼ $0$ であり、すなわち採択率はほぼ $1$ になる。HMCの効率の良さはここから来てるんやね。ちなみにここでハミルトン力学とMetropolisアルゴリズムが融合するので、この手法はHybrid MCと呼ばれたりもする。どっちにしろHMCである。

4-4. 結局

RMC法におけるアルゴリズムは、以下のように表現される。

初期値 $\boldsymbol{\theta}^{(0)}$ 、ハイパーパラメタ $\varepsilon$ を定める。
ポテンシャルエネルギー $h(\boldsymbol{\theta})$ を目的分布 $\text{Pr}(\boldsymbol{\theta} | \boldsymbol{x})$ の負の対数として定義する。
$\tau = 1, 2, \cdots, T$ において、以下のSTEP 4, 5, 6, 7を繰り返す。
$\boldsymbol{p}^{(\tau-1)}$ を標準正規分布 $\mathcal{N}(0,1)$ からサンプリングする。
後に示すリープフロッグ法を用いて $\boldsymbol{\theta}^{*}$ を計算する。
確率に従って採択する
- 採択された場合は、 $\boldsymbol{\theta}^{(\tau)} \leftarrow \boldsymbol{\theta}^{*}$ とする。
- 棄却された場合は、 $\boldsymbol{\theta}^{(\tau)} \leftarrow \boldsymbol{\theta}^{(\tau-1)}$ とする。つまりやり直し。
得られた $\boldsymbol{\theta}^{(\tau)}$ をサンプルとして保存する。

リープフロッグ法は以下の通り。

$\begin{align} \boldsymbol{p}^{(\tau-1/2)} &= \boldsymbol{p}^{(\tau-1)} - \frac{\varepsilon}{2} \frac{d h(\boldsymbol{\theta} ^{(\tau-1)})}{d \boldsymbol{\theta}} \\\ \boldsymbol{\theta}^{*} &= \boldsymbol{\theta}^{(\tau-1)} + \varepsilon \boldsymbol{p}^{( \tau - 1/2)} \\\ \boldsymbol{p}^{*} &= \boldsymbol{p}^{( \tau-1/2 )} - \frac{\varepsilon}{2} \frac{d h(\boldsymbol{\theta}^{*})}{d \boldsymbol{\theta}} \end{align}$

また採択率は以下の通りとなる。

$\begin{align} r^{(\tau)} = \min \left\lbrace 1, \exp \left( H(\boldsymbol{\theta}^{(\tau-1)}, \boldsymbol{p}^{(\tau-1)}) - H(\boldsymbol{\theta}^{*}, \boldsymbol{p}^{*}) \right) \right\rbrace \end{align}$

まとめ

HMCについてまとめた。まとめられてないけど。多分一生腹落ちすることはないと思う。実装してみたい。

おまけ

かえるぴょこぴょこ。これまで変分ベイズの記事が最大文字数(1.5万文字ぐらい)だったが、この記事はそれを大きく上回り2.1万字となりました。

変分ベイズも見てやってね。