Hamiltonian Monte Carlo (HMC) in Python: Overcoming the Metropolis Algorithm's Weaknesses with Leapfrog Integration

Introduction

https://yuhi-sa.github.io/en/posts/20260226_mcmc/1/ covered the Metropolis-Hastings algorithm and Gibbs sampling, noting practical caveats like “consecutive samples are correlated” and “the acceptance rate needs tuning.” But on distributions with high dimensionality and strong correlation between variables, this autocorrelation problem becomes severe, and a random-walk proposal struggles to explore efficiently. This article derives and implements Hamiltonian Monte Carlo (HMC), which constructs proposals from Hamiltonian dynamics using gradient information, and quantitatively measures how much it improves sampling efficiency over basic Metropolis.

The Metropolis Algorithm’s Weakness: The Random Walk

The Metropolis-Hastings algorithm from https://yuhi-sa.github.io/en/posts/20260226_mcmc/1/ proposes an isotropic random step from the current location and accepts or rejects it based on the target distribution. This “random walk” style exploration is highly inefficient on distributions where variables are strongly correlated (e.g., a 2D Gaussian with correlation 0.95). Because the high-probability region only exists along a narrow, elongated “valley,” most randomly-directed proposals point toward low-probability territory and get rejected — and even accepted ones only creep a small distance along the valley.

HMC: Efficient Proposals Using Gradients

HMC treats the log probability density \(\log p(\mathbf{x})\) as the negative of a “potential energy,” introduces an auxiliary “momentum” variable \(\mathbf{p}\) , and generates proposals by evolving the system according to Hamiltonian dynamics from physics.

Evolving \((\mathbf{x}, \mathbf{p})\) over time according to Hamilton’s equations of motion

conserves the energy \(H\) almost exactly (up to numerical integration error). Because the momentum accelerates in the direction of the gradient \(\nabla \log p(\mathbf{x})\) , the trajectory can travel far along high-density regions, sidestepping the inefficient random-walk exploration.

Leapfrog Integration

To numerically solve equation (2) while preserving energy conservation and reversibility, we use leapfrog integration:

Repeating this \(L\) times, we then apply the Metropolis criterion

to accept or reject. Since numerical integration error means \(H\) isn’t perfectly conserved, this Metropolis correction guarantees convergence to the correct target distribution.

Python Implementation

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 implements the iteration in equation (3); hmc’s acceptance test implements the Metropolis correction in equation (4).

Numerical Experiment: Comparison Against Metropolis on a Correlated Gaussian

Using a 2D Gaussian target \(\mathcal{N}(\mathbf{0}, \Sigma)\) with correlation \(\rho=0.95\) , \(\Sigma = \begin{pmatrix}1 & 0.95 \\ 0.95 & 1\end{pmatrix}\) , we drew 5000 samples each with Metropolis (from https://yuhi-sa.github.io/en/posts/20260226_mcmc/1/, step size tuned to roughly 38% acceptance) and HMC (20 leapfrog steps).

Acceptance Rates

Effective Sample Size (ESS)

HMC was roughly 45x more efficient than Metropolis — the same 5000 samples carry roughly 45x as much information, measured in independent-sample-equivalents, as the Metropolis chain. HMC’s accepted samples carry essentially no autocorrelation, behaving almost like independent draws.

Covariance Recovery Accuracy (after discarding 500 burn-in samples)

Even with 4500 post-burn-in samples, Metropolis underestimates the variance by roughly 15-20%. This reflects slow movement along the strongly correlated valley, meaning the chain hasn’t adequately explored the full posterior. HMC, using the same sample budget, recovers a covariance estimate that nearly matches the truth.

Why HMC/NUTS Is the Modern Default

This efficiency gap is exactly why modern probabilistic programming languages — Stan, PyMC, NumPyro — default to HMC (or its auto-tuning extension, NUTS: the No-U-Turn Sampler) as their sampler. When working with high-dimensional, highly-correlated posteriors in practice (e.g., hierarchical Bayesian models), basic Metropolis often can’t produce converged samples in a reasonable amount of time, making gradient-based HMC-family methods close to essential.

Related Articles

• Markov Chain Monte Carlo (MCMC) Fundamentals: Metropolis and Gibbs Sampling - This article follows up on the autocorrelation problem noted there, solving it with a gradient-based method.
• Bayesian Optimization: Theory and Python Implementation - Relevant to using MCMC-estimated posteriors for Bayesian optimization.
• Gaussian Process Regression (GPR) in Practice - Fully Bayesian hyperparameter inference often relies on efficient MCMC methods like HMC.
• Bayesian Linear Regression Fundamentals - A contrast with analytic (non-MCMC) Bayesian inference.
• SGD and Adam: Theory and Comparison - HMC’s use of gradients has a conceptual kinship with gradient-descent optimizers; comparing the two deepens intuition for both.

References

• 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.