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は、同じサンプル数で真の共分散にほぼ一致する推定を得ました。

トレースプロットと自己相関の可視化

上記のESS・共分散の数値差が、実際のサンプル列としてどう現れるかを可視化します。HMC（ステップサイズ0.2、リープフロッグ20ステップ、受容率98.0%）とメトロポリス法（提案標準偏差0.6、受容率37.6%）でそれぞれ5000サンプルを生成し（バーンイン500除去後の4500サンプル、seed=20260715）、\(x_1\) 座標のトレースプロットと自己相関関数（ACF）を比較しました。

def metropolis(log_prob, n_samples, sigma, x0, rng):
    samples = np.zeros((n_samples, len(x0)))
    x = x0.copy()
    n_accept = 0
    for i in range(n_samples):
        x_new = x + sigma * rng.standard_normal(len(x0))
        if np.log(rng.uniform()) < log_prob(x_new) - log_prob(x):
            x = x_new
            n_accept += 1
        samples[i] = x
    return samples, n_accept / n_samples


def acf(x, max_lag):
    n = len(x)
    x = x - x.mean()
    var = np.dot(x, x) / n
    result = np.zeros(max_lag + 1)
    for k in range(max_lag + 1):
        result[k] = np.dot(x[: n - k], x[k:]) / ((n - k) * var)
    return result


def ess_geyer(x):
    """Geyerの初期正数列推定量によるESS（1+2*罰則付き自己相関和 = 積分自己相関時間）"""
    n = len(x)
    x = x - x.mean()
    var = np.dot(x, x) / n
    acf_full = np.correlate(x, x, mode="full")[n - 1 :]
    rho_k = (acf_full / (var * np.arange(n, 0, -1)))[1:]
    total, m = 0.0, (len(rho_k) // 2) * 2
    for k in range(0, m, 2):
        pair_sum = rho_k[k] + rho_k[k + 1]
        if pair_sum < 0:
            break
        total += pair_sum
    tau = max(1 + 2 * total, 1.0)
    return n / tau, tau


rng_hmc = np.random.default_rng(20260715)
samples_hmc, acc_hmc = hmc(log_prob, grad_log_prob, 5000, 0.2, 20, np.zeros(2), rng_hmc)
rng_mh = np.random.default_rng(20260715)
samples_mh, acc_mh = metropolis(log_prob, 5000, 0.6, np.zeros(2), rng_mh)

burn = 500
ess_hmc, tau_hmc = ess_geyer(samples_hmc[burn:, 0])
ess_mh, tau_mh = ess_geyer(samples_mh[burn:, 0])
print(f"HMC: accept={acc_hmc*100:.1f}%  tau_int={tau_hmc:.3f}  ESS={ess_hmc:.1f}")
print(f"MH : accept={acc_mh*100:.1f}%  tau_int={tau_mh:.3f}  ESS={ess_mh:.1f}")

HMC: accept=98.0%  tau_int=1.000  ESS=4500.0
MH : accept=37.6%  tau_int=73.692  ESS=61.1

トレースプロット（左図）では、HMC（青）が\(x_1\) 軸上を高速に往復しながら分布全域を埋めているのに対し、メトロポリス法（赤）は近傍にとどまる「粘着」した動きを示しています。自己相関関数（右図）ではこの違いがより明確で、HMCはラグ1で符号が反転し振動しながら急速にゼロへ収束するのに対し、メトロポリス法はラグ60でも自己相関0.2程度が残る緩やかな単調減衰を示します。この自己相関の振動（負の自己相関）は、リープフロッグ軌道が1回の提案で分布の谷を大きく横切ることの直接的な帰結であり、積分自己相関時間 \(\tau_{\mathrm{int}}\) がHMCで1.000（理論的下限に到達＝ほぼ独立サンプル）、メトロポリス法で73.7と、73.7倍の差になって現れています（この個別測定は、共分散復元精度の節で報告した45倍という数値と手法・サンプル系列は同一ですが、ESS推定量が異なるため厳密には一致しません。ともに数十倍規模の効率差を示す点で整合的です）。

発散する軌道：ステップサイズが大きすぎる場合

ここまではステップサイズ \(\epsilon\) を「適切に選ばれた値」として扱ってきましたが、実務でHMCを使う際に最初にぶつかる問題は、ステップサイズが大きすぎるとリープフロッグ積分がエネルギー保存性を失い、シミュレーションが「発散」することです。Stan・PyMCなどが警告する “divergent transitions” はまさにこの現象を指します（Betancourt, 2017）。

なぜ発散するのか：線形安定性解析

目標分布がガウス分布 \(\mathcal{N}(\mathbf{0}, \Sigma)\) の場合、リープフロッグ積分は精度行列 \(\Lambda = \Sigma^{-1}\) の固有ベクトル方向に分解でき、各固有値 \(\lambda_i\) の方向は角振動数 \(\omega_i = \sqrt{\lambda_i}\) の単振動として振る舞います。リープフロッグ（2次の陽的解法）はこの調和振動子に対して線形安定性の限界を持ち、次の条件を超えると数値誤差が指数的に増大します。

ここで \(\lambda_{\max}\) は精度行列 \(\Lambda\) の最大固有値です。本記事の相関ガウス分布（\(\rho=0.95\) ）では、精度行列の固有値は約 \(0.513\) と \(20.0\) なので、\(\epsilon_{\mathrm{crit}} = 2/\sqrt{20.0} \approx 0.447\) と予測されます。

数値実験：ステップサイズを増やしながらエネルギー誤差と受容率を追跡

\(L=20\) に固定し、ステップサイズを \(0.05\) から \(1.3\) まで変化させ、各ステップサイズで3000回の提案を生成して、平均エネルギー誤差 \(\overline{|\Delta H|}\) （提案前後のハミルトニアンの絶対差の平均）と受容率を記録しました（seed=20260715）。

def hmc_energy_error(log_prob, grad_log_prob, n_samples, step_size, n_leapfrog, x0, rng):
    x = x0.copy()
    n_accept = 0
    abs_dH = np.zeros(n_samples)
    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
        with np.errstate(over="ignore", invalid="ignore"):
            proposed_H = -log_prob(x_new) + 0.5 * p_new @ p_new
        if not np.isfinite(proposed_H):
            proposed_H = 1e12  # オーバーフロー時は「棄却確定」の大きな値で代用
        abs_dH[i] = abs(proposed_H - current_H)
        if np.log(rng.uniform()) < current_H - proposed_H:
            x, n_accept = x_new, n_accept + 1
    return n_accept / n_samples, abs_dH.mean()


eigvals = np.linalg.eigvalsh(Sigma_inv)
eps_crit = 2 / np.sqrt(eigvals.max())
print(f"精度行列の固有値: {eigvals}, 理論的発散限界 eps_crit={eps_crit:.3f}")

for eps in [0.05, 0.1, 0.2, 0.3, 0.4, 0.45, 0.5, 0.6, 0.7, 0.85, 1.0, 1.3]:
    rng = np.random.default_rng(20260715)
    acc, mean_dH = hmc_energy_error(log_prob, grad_log_prob, 3000, eps, 20, np.zeros(2), rng)
    print(f"eps={eps:5.2f}  accept={acc*100:6.2f}%  mean|dH|={mean_dH:.4g}")

精度行列の固有値: [ 0.51282051 20.        ], 理論的発散限界 eps_crit=0.447
eps= 0.05  accept= 99.73%  mean|dH|=0.007922
eps= 0.10  accept= 99.57%  mean|dH|=0.01297
eps= 0.20  accept= 98.17%  mean|dH|=0.03868
eps= 0.30  accept= 82.27%  mean|dH|=0.3823
eps= 0.40  accept= 83.20%  mean|dH|=0.3629
eps= 0.45  accept=  0.23%  mean|dH|=7.751e+04
eps= 0.50  accept=  0.00%  mean|dH|=3.325e+16
eps= 0.60  accept=  0.00%  mean|dH|=2.599e+27
eps= 0.70  accept=  0.00%  mean|dH|=5.279e+34
eps= 0.85  accept=  0.00%  mean|dH|=8.673e+42
eps= 1.00  accept=  0.00%  mean|dH|=2.279e+49
eps= 1.30  accept=  0.00%  mean|dH|=1.731e+59

理論的な発散限界 \(\epsilon_{\mathrm{crit}} \approx 0.447\) を境に、平均エネルギー誤差が0.36から77,506へ、5桁近くジャンプし、受容率は83.2%から0.23%へ崩壊しています。\(\epsilon=0.5\) 以降はエネルギー誤差が \(10^{16}\) を超え、受容率は完全に0%になります。これは、線形安定性解析で予測した通り、リープフロッグ積分が指数的に発散するモード（固有値20の「速い」方向）を含んでいるためです。実務上の教訓は次の2点です。

1. 受容率の急落は「ステップサイズ過大」の明確なシグナルであり、受容率を見ながらステップサイズを調整する（デュアル平均化などの適応則）ことは、発散を避けるための最低限の安全策になります。
2. エネルギー誤差そのものを監視する方がより早期に問題を検知できる：受容率が顕著に落ちる前（\(\epsilon=0.3\sim0.4\) ）でも、エネルギー誤差はすでに小さいステップサイズの40倍程度に増大しており、Stan/PyMCが「divergent transition」として個々の提案を検出・警告する設計（単純な受容率のモニタリングより感度が高い）の合理性を裏付けています。

ステップサイズとリープフロッグステップ数の感度：ESS/勾配評価という指標

ここまで見た受容率とESSは「1提案あたりの質」の指標ですが、HMCの計算コストはリープフロッグステップ数 \(L\) に比例します（1提案につき \(L+1\) 回の勾配評価）。したがって、サンプル数あたりのESS（本記事の冒頭で報告した100%やメトロポリスの2.20%）だけでは計算効率を正しく評価できません。ステップ数を増やせばESSは改善しやすくなりますが、その分1提案のコストも増えるため、真に比較すべき指標は**勾配評価1回あたりのESS（ESS/勾配評価数）**です。

数値実験：ステップサイズ × リープフロッグステップ数のグリッドサーチ

ステップサイズ \(\epsilon \in \{0.05, 0.10, \ldots, 0.40\}\) （発散が始まる\(0.447\) 未満に限定）とリープフロッグステップ数 \(L \in \{5, 10, 20, 30, 40, 60, 80, 120, 160, 200\}\) の全80通りの組み合わせで、それぞれ3000サンプル（バーンイン300、seed=20260715）を生成し、Geyerの初期正数列推定量でESSを計算、\(L+1\) 回の勾配評価で正規化しました。

step_sizes = [0.05, 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40]
n_leapfrogs = [5, 10, 20, 30, 40, 60, 80, 120, 160, 200]
n_samples, burn = 3000, 300

results = []
for eps in step_sizes:
    for L in n_leapfrogs:
        rng = np.random.default_rng(20260715)
        samples, acc = hmc(log_prob, grad_log_prob, n_samples, eps, L, np.zeros(2), rng)
        s = samples[burn:]
        ess_x0, _ = ess_geyer(s[:, 0])
        ess_x1, _ = ess_geyer(s[:, 1])
        ess_min = min(ess_x0, ess_x1)
        grad_evals = n_samples * (L + 1)
        results.append((eps, L, acc, ess_min, ess_min / grad_evals))

抜粋（全80通りのうち代表的な行、単位は ESS/勾配評価 \(\times 10^{-3}\) ）：

このヒートマップから2つの重要な事実が読み取れます。

（1）「ステップ数を増やせば増やすほど良い」は誤り。今回のグリッドでの最良点は \(\epsilon=0.25, L=10\) で、ESS/勾配評価は \(81.8\times10^{-3}\) でした。\(L=10\) 前後ですでにほぼ独立なサンプル（ESS \(\approx\) サンプル数）が得られており、\(L\) をそれ以上増やしても1提案あたりのESSはほとんど改善しない一方、勾配評価コストは線形に増え続けるため、ESS/勾配評価は単調に悪化します。

（2）\((\epsilon, L)\) の組み合わせによっては「受容率は高いのにESSが崩壊する」危険な組み合わせが存在する。例えば \(\epsilon=0.15, L=60\) では受容率99.3%と一見健全に見えますが、ESS(最小次元)はわずか22.1（同じ \(\epsilon\) で \(L=20\) や \(L=40\) なら2700.0＝サンプル数上限に達する）まで落ち込みます。これは、軌道長 \(L\epsilon\) が調和振動子の周期（またはその整数倍）に近づくと、リープフロッグ軌道がほぼ同じ位相に戻ってしまい、提案が事実上「その場に留まる」ことに起因する共鳴（periodicity）現象です。Neal (2011) はこの現象を理論的に指摘し、固定した \(L\) と \(\epsilon\) の組み合わせがこの共鳴に偶然一致すると、たとえ受容率が高くてもチェーンがエルゴード性を失いうることを示し、対策として**\(\epsilon\) または \(L\) を毎回の提案でランダムに選ぶこと**を提案しています。

この「\(L\) の最適値は目標分布の形状に依存し、かつ固定すると共鳴という思わぬ落とし穴がある」という事実こそが、次節で説明する**NUTS（No-U-Turn Sampler）**が \(L\) を自動決定する仕組みへと発展した動機です。

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

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

NUTSが解決する問題：リープフロッグステップ数の手動チューニング

前節で見た通り、リープフロッグステップ数 \(L\) の最適値は目標分布の形状（固有値の分布）に依存し、しかも固定すると共鳴による性能崩壊のリスクがあります。\(L\) が小さすぎればランダムウォークに近い非効率な探索に戻り（\(L=5\) でESS/勾配評価は最適点の20分の1以下）、大きすぎれば計算コストが無駄になるだけでなく、軌道が目標分布の周りを回りすぎて出発点付近に戻ってくる「Uターン」が発生し、獲得した勾配評価の多くが無駄になります。

Hoffman & Gelman (2014) が提案したNUTSは、この\(L\) の手動チューニングを自動化します。具体的には、現在地から前方・後方に軌道を倍々に伸ばしながら二分木を構築し、木のどちらかの端で運動量と位置ベクトルの内積が負に転じる（軌道が引き返し始める＝Uターンする）時点で木の伸長を打ち切ります。これにより、目標分布の局所的な形状に応じて\(L\) を毎回自動的に決定でき、実務家がステップサイズ以外のハイパーパラメータを手動で探索する必要がなくなります。

さらに近年の研究は、NUTSがまだ抱える限界（ステップサイズ \(\epsilon\) 自体は依然としてグローバルに固定される点や、階層モデルの「くびれ」構造で局所的に必要な分解能が変わる点）に取り組んでいます。Bou-Rabee, Carpenter & Marsden (2024) のGISTは、Gibbsサンプリングの枠組みでステップサイズや軌道長を位置・運動量に応じて局所適応させる一般的な理論的枠組みを提示し、NUTSやランダム化HMCをその特殊ケースとして統一的に説明しています。Modi (2024) のATLASは、局所的なヘッセ行列の低ランク近似からステップサイズを反復ごとに適応させ、Uターン条件の監視で軌道長も適応させることで、多重スケールの曲率を持つ事後分布での発散転移を大幅に削減しています。これらは、本記事で数値的に確認した「\(\epsilon\) と \(L\) の同時チューニングが難しく、固定値では危険な組み合わせが存在する」という問題に対する、2024年時点での最先端の解答です。

関連記事

• マルコフ連鎖モンテカルロ（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.
• Bou-Rabee, N., Carpenter, B., & Marsden, M. (2024). GIST: Gibbs self-tuning for locally adaptive Hamiltonian Monte Carlo. arXiv:2404.15253.
• Modi, C. (2024). ATLAS: Adapting trajectory lengths and step-size for Hamiltonian Monte Carlo. arXiv:2410.21587.