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

分子動力学法を用いたMH法であるハミルトニアン・モンテカルロ法を導入する。

「ハミルトニアン」を用いたM-H法導出分子動力学法の条件受理確率はハミルトニアンの差分の関数になるアルゴリズムちなみに

「ハミルトニアン」を用いたM-H法

ハミルトニアン・モンテカルロ法は、M-H法の1つ。ギブスサンプリングと同様に提案分布 $Q(\bm y | \bm x)$ を工夫したものであるが、こちらは提案分布が必ず $1$ になるものではなく、できるだけ $1$ に近づけるようなダイナミクスを用いている。その工夫が $\bm y \sim Q(\bm y | \bm x)$ という操作に分子動力学法による時間発展を加えるというものであり、「エネルギー保存則」を近似的に利用することで提案分布を1に近づけている。

…という説明をされると難しそうに聞こえるが、1つ1つ読み解いていけばそれほどややこしいアルゴリズムではないとわかるはずだ。

導出

確率分布 $P(\bm x)$ が以下の形をしているとする。

\begin{aligned} P(\bm x) = \frac{1}{Z} \exp(-S(\bm x)) \end{aligned}

ここで $S(\bm x)$ は「作用」と呼ばれる関数である。ここに「共役運動量」と呼ばれる新たな変数 $\bm p$ を導入し、次のようなハミルトニアンを定義する。

\begin{aligned} H(\bm x, \bm p) \coloneqq S(\bm x) + \frac{1}{2} \|\bm p\|_2^2 \end{aligned}

ただし共役運動量 $\bm p$ は多変量標準正規分布からサンプルされたものとする。

\begin{aligned} \bm p \sim {\cal N}_N (\bm p | \bm 0, \bm I) \end{aligned}

このハミルトニアンを「分子動力学法」と呼ばれるアルゴリズムにより時間発展させ、時間発展後の $\bm x$ を、提案 $\bm y$ として扱うことにする。

H(\bm x, \bm p) \to H(\bm y, \bm q)

分子動力学法の条件

分子動力学法による時間発展 $H(\bm x, \bm p) \to H(\bm y, \bm q)$ に伴う変化およびその出力を次のように表すことにする。

\begin{alignedat}2 (\bm y, \bm q) \gets{}& \operatorname{MD}(\bm x, \bm p) \end{alignedat}

「分子動力学法」の具体的なアルゴリズムはここには書かないが、以下の2つの要件は重要なので記しておく。

共役運動量 $\bm p$ を逆向きにすることで、逆向きの時間発展が可能でなくてはならない:
$\begin{aligned} (\bm y, \bm q) &= \operatorname{MD}(\bm x, \bm p) \\ (\bm x, -\bm p) &= \operatorname{MD}(\bm y, -\bm q) \end{aligned}$
時間発展前後におけるハミルトニアンの差分 $\Delta H = H(\bm y, \bm q) - H(\bm x, \bm p)$ はできるだけ小さいほうが好ましい:
$| \Delta H | \ll 1$

受理確率はハミルトニアンの差分の関数になる

以上の提案アルゴリズムは、提案分布を次のように設定することに相当する。

\begin{aligned} Q(\bm y | \bm x) &\coloneqq \left\{\begin{aligned} & {\cal N}_N (\bm p | \bm 0, \bm I) && \text{if } (\bm y, \bm q) = \operatorname{MD}(\bm x, \bm p) \\ & 0 && \text{otherwise} \end{aligned}\right. \\ Q(\bm x | \bm y) &\coloneqq \left\{\begin{aligned} & {\cal N}_N (-\bm q | \bm 0, \bm I) && \text{if } (\bm x, -\bm p) = \operatorname{MD}(\bm y, -\bm q) \\ & 0 && \text{otherwise} \end{aligned}\right. \end{aligned}

これを用いて受理確率を計算すると、その値はハミルトニアンの差分 $\Delta H = H(\bm y, \bm q) - H(\bm x, \bm p)$ の値に依存することがわかる。

\begin{aligned} A(\bm y | \bm x) ={}& \frac{ P(\bm y) Q(\bm x | \bm y) }{ P(\bm x) Q(\bm y | \bm x) } \\ &\left|\small\quad\begin{aligned} Q(\bm y | \bm x) &= {\cal N}_N (\bm p | \bm 0, \bm I) = \sqrt\frac{1}{(2\pi)^N} \exp \left( -\frac{1}{2} \|\bm p\|_2^2 \right) \\ Q(\bm x | \bm y) &= {\cal N}_N (-\bm q | \bm 0, \bm I) = \sqrt\frac{1}{(2\pi)^N} \exp \left( -\frac{1}{2} \|\bm q\|_2^2 \right) \\ \end{aligned}\right. \\ ={}& \frac{ P(\bm y) \exp \left( -\dfrac{1}{2} \|\bm p\|_2^2 \right) }{ P(\bm x) \exp \left( -\dfrac{1}{2} \|\bm q\|_2^2 \right) } \\ &\left|\small\quad\begin{aligned} P(\bm y) ={}& \frac{1}{Z_2} \exp(-S(\bm y)) \\ P(\bm x) ={}& \frac{1}{Z_2} \exp(-S(\bm x)) \\ \end{aligned}\right. \\ ={}& \frac{ \exp(-S(\bm y)) \exp \left( -\dfrac{1}{2} \|\bm p\|_2^2 \right) }{ \exp(-S(\bm x)) \exp \left( -\dfrac{1}{2} \|\bm q\|_2^2 \right) } \\ ={}& \frac{ \exp\left(-S(\bm y) - \dfrac{1}{2} \|\bm p\|_2^2 \right) }{ \exp\left(-S(\bm x) - \dfrac{1}{2} \|\bm q\|_2^2 \right) } \\ ={}& \frac{ \exp(-H (\bm y, \bm q)) }{ \exp(-H (\bm x, \bm p)) } \\ ={}& \exp( -(H (\bm y, \bm q) - H (\bm x, \bm p)) ) \\ ={}& \exp( -\Delta H ) \end{aligned}

差分 $\Delta H$ が $0$ に近いほど好ましいのは、そのほうが受理確率が $1$ に近づくためである。

アルゴリズム

以上の考察から、次のようなアルゴリズムを得る。

Algorithm: Hamiltonian Monte Carlo.
$\begin{darray}{rl} 1. & X \gets \{\} \\ 2. & \text{for } t \in \{1, 2, \dots, n\} : \\ 3. & \qquad \bm p \sim \mathcal N_N (\bm p | \bm 0, \bm I) \\ 4. & \qquad (\bm y, \bm q) \gets \operatorname{MD}(\bm x, \bm p) \\ 5. & \qquad \Delta H \gets H (\bm y, \bm q) - H (\bm x, \bm p) \\ 6. & \qquad r \sim \mathcal U(r | 0, 1) \\ 7. & \qquad \text{if } \exp(-\Delta H) \gt r: \\ 8. & \qquad \qquad \bm x \gets \bm y \\ 9. & \qquad X \gets X \cup \{ \bm x \} \\ 10. & \text{return } X \end{darray}$

このようなM-H法の実装を特にハミルトニアン・モンテカルロ法 (Hamiltonian Monte Carlo; HMC法) という。

ちなみに

ハミルトニアン・モンテカルロ法 (HMC法) はハイブリッド・モンテカルロ法 (hybrid Monte Carlo; HMC) と呼ばれることもある。