Hamiltonian Monte Carlo（HMC）の理論とPython実装：リープフロッグ積分でMetropolis法の弱点を克服する

はじめに

https://yuhi-sa.github.io/posts/20260226_mcmc/1/ ではメトロポリス・ヘイスティングス法とギブスサンプリングを解説し、「連続するサンプルは相関を持つ」「受容率の調整が必要」という実用上の注意点に触れました。しかし、次元が高く変数間の相関が強い分布では、この自己相関の問題が深刻化し、ランダムウォーク型の提案分布では効率的な探索が困難になります。本記事では、**勾配情報を使ってハミルトン力学に基づく提案分布を構成するHamiltonian Monte Carlo（HMC）**を導出・実装し、基本のメトロポリス法と比較してどれだけサンプリング効率が改善するかを定量的に検証します。

メトロポリス法の弱点：ランダムウォーク

https://yuhi-sa.github.io/posts/20260226_mcmc/1/ のメトロポリス・ヘイスティングス法は、現在地から等方的なランダムな一歩を提案し、目標分布に基づいて受理・棄却します。この「酔歩（ランダムウォーク）」的な探索は、変数間の相関が強い分布（例えば相関係数0.95の2次元ガウス分布）では非常に非効率です。細長い「谷」に沿ってしか高確率領域が存在しないため、ランダムな方向への提案のほとんどが低確率領域に向かい棄却されるか、受理されても谷に沿ってわずかしか進めません。

HMC：勾配を使った効率的な提案

HMCは、対数確率密度 \(\log p(\mathbf{x})\) を「ポテンシャルエネルギー」の負値とみなし、補助的な「運動量」変数 \(\mathbf{p}\) を導入して、物理学のハミルトン力学に従って提案を生成します。

ハミルトン力学の運動方程式

に従って \((\mathbf{x}, \mathbf{p})\) を時間発展させると、エネルギー \(H\) はほぼ保存されます（数値積分の誤差を除く）。勾配 \(\nabla \log p(\mathbf{x})\) の方向に運動量が加速されるため、確率密度の高い領域に沿って大きく移動でき、ランダムウォークのような非効率な探索を回避できます。

リープフロッグ積分

式(2)を数値的に解くには、エネルギー保存性・可逆性を保つリープフロッグ積分を使います。

これを \(L\) ステップ繰り返し、最後にメトロポリス基準

で受理・棄却を判定します。数値積分の誤差により \(H\) が完全には保存されないため、このメトロポリス補正で正しい目標分布への収束を保証します。

Python実装

import numpy as np

def leapfrog(x, p, grad_log_prob, step_size, n_steps):
    p = p + 0.5 * step_size * grad_log_prob(x)
    for _ in range(n_steps - 1):
        x = x + step_size * p
        p = p + step_size * grad_log_prob(x)
    x = x + step_size * p
    p = p + 0.5 * step_size * grad_log_prob(x)
    return x, p

def hmc(log_prob, grad_log_prob, n_samples, step_size, n_leapfrog, x0, rng):
    samples = np.zeros((n_samples, len(x0)))
    x = x0.copy()
    n_accept = 0
    for i in range(n_samples):
        p0 = rng.standard_normal(len(x0))
        x_new, p_new = leapfrog(x.copy(), p0.copy(), grad_log_prob, step_size, n_leapfrog)
        current_H = -log_prob(x) + 0.5 * p0 @ p0
        proposed_H = -log_prob(x_new) + 0.5 * p_new @ p_new
        if np.log(rng.uniform()) < current_H - proposed_H:
            x = x_new
            n_accept += 1
        samples[i] = x
    return samples, n_accept / n_samples

leapfrog が式(3)の反復、hmc の受理判定が式(4)のメトロポリス補正に対応します。

数値実験：相関ガウス分布でMetropolisとの比較

目標分布として、相関係数 \(\rho=0.95\) の2次元ガウス分布 \(\mathcal{N}(\mathbf{0}, \Sigma)\) 、\(\Sigma = \begin{pmatrix}1 & 0.95 \\ 0.95 & 1\end{pmatrix}\) を使い、https://yuhi-sa.github.io/posts/20260226_mcmc/1/ のメトロポリス法（ステップサイズを調整し受容率約38%）とHMC（リープフロッグ20ステップ）でそれぞれ5000サンプル生成しました。

受容率

有効サンプルサイズ（ESS）

HMCはメトロポリス法よりも約45倍効率的（同じ5000サンプルから、独立サンプル換算でメトロポリスの45倍の情報量を得られる）という結果になりました。HMCの受理サンプルはほぼ自己相関を持たず、実質的に独立サンプルに近い性質を持っています。

共分散復元精度（バーンイン500サンプル除去後）

メトロポリス法は4500サンプル（バーンイン後）を使っても、分散を約15〜20%過小評価しています。これは、相関の強い谷に沿った移動が遅く、事後分布の全域を十分に探索できていないことを示しています。一方HMCは、同じサンプル数で真の共分散にほぼ一致する推定を得ました。

なぜHMC/NUTSが現代の標準なのか

この効率差が、Stan・PyMC・NumPyroといった現代の確率的プログラミング言語がHMC（の自動チューニング拡張であるNUTS: No-U-Turn Sampler）をデフォルトサンプラーに採用している理由です。実務で高次元・高相関の事後分布（階層ベイズモデルなど）を扱う場合、基本のメトロポリス法では現実的な時間で収束したサンプルを得るのが困難なことが多く、勾配情報を使うHMC系の手法がほぼ必須になります。

関連記事

• マルコフ連鎖モンテカルロ（MCMC）の基礎：メトロポリス法とギブスサンプリング - 本記事はこの基礎記事で触れられた自己相関の問題を、勾配ベースの手法で解決する発展編です。
• ベイズ最適化の基礎とPython実装 - MCMCで推定した事後分布をベイズ最適化に活用する文脈で関連します。
• ガウス過程回帰（GPR）の実践 - ハイパーパラメータの完全ベイズ推定にはHMCのような効率的なMCMC手法が使われます。
• ベイズ線形回帰の基礎 - MCMCを使わない解析的なベイズ推定との対比になります。
• SGD・Adamの理論と比較 - HMCの勾配活用はSGDの勾配降下法と発想が近く、対比すると理解が深まります。

参考文献

• Neal, R. M. (2011). MCMC using Hamiltonian dynamics. In Handbook of Markov Chain Monte Carlo. Chapman & Hall/CRC.
• Hoffman, M. D., & Gelman, A. (2014). The No-U-Turn Sampler: Adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15(1), 1593-1623.
• Betancourt, M. (2017). A conceptual introduction to Hamiltonian Monte Carlo. arXiv:1701.02434.