リープフロッグアルゴリズム

HMC法で用いられる分子動力学法であるリープフロッグ法を解説する。

馬跳法

HMC法では逆向きの時間発展が可能であるような時間発展アルゴリズム

\begin{aligned} (\bm y, \bm q) &= \operatorname{MD}(\bm x, \bm p) \\ (\bm x, -\bm p) &= \operatorname{MD}(\bm y, -\bm q) \end{aligned}

によって新たな提案点を得る。このアルゴリズムの主要なものにリープフロッグ積分法 (Leapfrog integration; 馬跳法) という数値積分法がある。変数 $\bm x, \bm p$ が「馬跳び」のように相手を飛び越すような挙動を示すことからこの名がついている。

アルゴリズム

以下、ハミルトニアンが次式で表されているとする。

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

リープフロッグ法は次のようなアルゴリズムによりハミルトニアンを時間発展させる。

Algorithm: Leapfrog algorithm.
$\begin{darray}{rl} 1. & \bm p_0 \gets \bm p \\ 2. & \bm x_0 \gets \bm x + \frac{1}{2} \Delta t \left.\frac{\partial \cal H}{\partial \bm p}\right|_{\bm p_0, \bm x} \\ 3. & \text{for } t \in \{1, 2, \dots, \tau-1\}: \\ 4. & \qquad \bm p_t \gets \bm p_{t-1} - \Delta t \left.\frac{\partial \cal H}{\partial \bm x} \right|_{\bm p_{t-1}, \bm x_{t-1}} \\ 5. & \qquad \bm x_t \gets \bm x_{t-1} + \Delta t \left.\frac{\partial \cal H}{\partial \bm x} \right|_{\bm p_{t}, \bm x_{t-1}} \\ 6. & \bm p_T \gets \bm p_{T-1} - \Delta t \left.\frac{\partial \cal H}{\partial \bm x} \right|_{\bm p_{T-1}, \bm x_{T-1}} \\ 7. & \bm x_T \gets \bm x_{T-1} + \frac{1}{2} \Delta t \left.\frac{\partial \cal H}{\partial \bm x} \right|_{\bm p_{T}, \bm x_{T-1}} \\ 8. & \text{return} (\bm p_T, \bm x_T) \\ \end{darray}$

ここで、

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

であることから

\begin{aligned} \frac{\partial H}{\partial \bm x} ={}& \frac{d S}{d \bm x} \\ \frac{\partial H}{\partial \bm p} ={}& \bm p \end{aligned}

が得られる。これを用いれば次のようなアルゴリズムになる。実用上はこちらを用いれば良い。

Algorithm: Leapfrog algorithm.
$\begin{darray}{rl} 1. & \bm p_0 \gets \bm p \\ 2. & \bm x_0 \gets \bm x + \frac{1}{2} \Delta t \bm p_0 \\ 3. & \text{for } t \in \{1, 2, \dots, \tau-1\}: \\ 4. & \qquad \bm p_t \gets \bm p_{t-1} - \Delta t \left.\frac{dS}{d \bm x} \right|_{\bm x_{t-1}} \\ 5. & \qquad \bm x_t \gets \bm x_{t-1} + \Delta t \bm p_{t} \\ 6. & \bm p_T \gets \bm p_{T-1} - \Delta t \left.\frac{dS}{d \bm x} \right|_{\bm x_{T-1}} \\ 7. & \bm x_T \gets \bm x_{T-1} + \frac{1}{2} \Delta t \bm p_{t} \\ 8. & \text{return} (\bm p_T, \bm x_T) \\ \end{darray}$

ハミルトンの運動方程式

時刻 $t$ における更新式に注目しよう。

\begin{aligned} \bm p_t ={}& \bm p_{t-1} - \Delta t \frac{\partial \cal H}{\partial \bm x} \\ \bm x_t ={}& \bm x_{t-1} + \Delta t \frac{\partial \cal H}{\partial \bm p} \\ \end{aligned}

これを変形すると次のようになる。

\begin{aligned} \frac{\bm p_t - \bm p_{t-1}}{\Delta t} ={}& {-}\frac{\partial H}{\partial \bm x} \\ \frac{\bm x_t - \bm x_{t-1}}{\Delta t} ={}& \frac{\partial H}{\partial \bm p} \\ \end{aligned}

さらに $\bm p_t = \bm p(t)$ および $\bm p_{t-1} = \bm p(t-\Delta t)$ と割り当てると、左辺はいずれも時間 $\Delta t$ における平均変化量になる。

\begin{aligned} \frac{\bm p(t) - \bm p(t-\Delta t)}{\Delta t} ={}& {-}\frac{\partial H}{\partial \bm x} \\ \frac{\bm x(t) - \bm x(t-\Delta t)}{\Delta t} ={}& \frac{\partial H}{\partial \bm p} \\ \end{aligned}

そして時間 $\Delta t$ を限りなく $+0$ に近づけていくと、次のハミルトンの運動方程式に帰着する。

\begin{aligned} \frac{d \bm p}{d t} &= -\frac{\partial H}{\partial \bm x} \\ \frac{d \bm x}{d t} &= \frac{\partial H}{\partial \bm p} \\ \end{aligned}

ハミルトニアン保存則

以上のようなことが成り立つ場合、ハミルトニアンを時間微分するとハミルトニアン保存則が成立する。

\begin{aligned} \frac{d H}{d t} ={}& \frac{d S(\bm x)}{dt} + \frac{1}{2} \frac{d}{dt}\|\bm p\|_2^2 \\ ={}& \frac{d \bm x}{d t}^\mathsf{T} \frac{d S(\bm x)}{d \bm x} + \frac{d \bm p}{dt}^\mathsf{T} \bm p \\ ={}& \frac{\partial H}{\partial \bm p}^\mathsf{T} \frac{d S(\bm x)}{d \bm x} - \frac{\partial H}{\partial \bm x}^\mathsf{T} \bm p \\ &\left|\small\quad\begin{aligned} \frac{\partial H}{\partial \bm x} &= \frac{d S}{d \bm x} \\ \frac{\partial H}{\partial \bm p} &= \bm p \\ \end{aligned}\right. \\ ={}& \bm p^\mathsf{T} \frac{d S(\bm x)}{d \bm x} - \frac{d S(\bm x)}{d \bm x}^\mathsf{T} \bm p \\ ={}& 0 \end{aligned}

すなわちハミルトニアンが時間依存しない。それゆえ、初期と終期のハミルトニアンが変化しないことになる。

\lim_{\Delta t \to +0} \Delta H = 0

したがって受理確率が $1$ になり、提案された $\bm y$ が必ず受理されることになる。

\begin{aligned} A(\bm y | \bm x) ={}& \exp(-\Delta H) \\ \to{}& \exp(0) \\ ={}& 1 \quad (\Delta t \to +0) \end{aligned}

もちろん $\Delta t$ を厳密に $+0$ に一致させることは不可能だし、もしもそのようなことをすると、そもそもこのアルゴリズムでは時間発展が行なわれなくなる。実際には、 $\Delta t$ は受理確率を $1$ に近づけるか、時間発展を速く行なうかのトレードオフを決めるパラメータとなる。