Hamiltonian Monte Carlo
Markov Chain Monte Carlo
What you need to know first 6 concepts, 2 layers
The requisite-knowledge inventory for this page, bottom-up: the primitives at the base, combined upward until you reach what this page assumes. Skim the layers you already own; start wherever the ground gets unfamiliar.
- base
- Bayes’ rule
- Hamiltonian mechanicsconcept
- Lennard-Jones MD (velocity Verlet)
- Linear algebraconcept
- MCMC from scratch (Metropolis)
- L1
- ↳you are here
2 of these are concepts without a dedicated page yet — the grey chips. Following the linked ones first makes the rest land.
Two pages on this site already contain everything Hamiltonian Monte Carlo needs. The Lennard-Jones molecular dynamics page integrates Newton's equations with leapfrog; the MCMC-from-scratch page corrects biased proposals with Metropolis. HMC is what you get when you point the MD integrator at a probability distribution and let Metropolis clean up after it — and it turns a sampler that crawls into one that flies. Five ways in.
Plain Metropolis proposes its next sample by jiggling the current one a little: with isotropic noise. On a target with strong correlations — a long thin ridge — that is a disaster. Big steps mostly point off the ridge and get rejected; small steps stay on it but inch along. The chain diffuses, random-walk style, and the distance it covers grows like , not .
The demo makes it concrete: on a ridge, random walk takes 127 steps to produce one effectively-independent sample. HMC takes one. The whole method exists to replace blind jiggling with a proposal that knows which way the ridge runs.
Turn the target into a landscape. Define a potential energy , so high-probability regions are valleys and low-probability regions are hills. A sampler that just wanted to sit at the bottom would be an optimizer; we want to explore the valley, not collapse into it.
So give the point a mass and a random shove. Now it is a ball: kick it in a random direction (sample a momentum) and let it roll under gravity-on-this-landscape for a while. It coasts along the valley floor, trading potential for kinetic energy and back, and comes to rest far from where it started — but still somewhere probable, because the ball can't climb a hill it lacks the energy for. Distant proposals, all in the valley. That is the entire idea; everything else is making it exact.
"Let the ball roll" means integrating Hamilton's equations, , . You already know the right integrator: leapfrog, the velocity-Verlet scheme from molecular dynamics — half-kick the momentum, drift the position, half-kick again.
It is chosen for a reason beyond familiarity: leapfrog is symplectic and exactly reversible, so it preserves phase- space volume and very nearly preserves energy even over long trajectories. Those two properties are precisely what the Metropolis correction needs to stay valid. The same code that ran your argon droplet proposes your Monte Carlo moves.
Why does a long, distant jump get accepted when a random-walk jump of the same length would be rejected? Because Hamiltonian dynamics conserve energy. The Metropolis acceptance probability is , and along an exact trajectory doesn't change, so the exponent is zero and the proposal is accepted with probability one — no matter how far it traveled.
Leapfrog isn't exact, so drifts by a small amount, and the accept/reject step exists exactly to absorb that discretization error and keep the samples unbiased. Shrink and acceptance climbs toward 100%; grow it and you trade acceptance for reach. In the demo, watch the rejected (red) trajectories appear only when you push too hard.
Here is why it's exact, not just clever. Invent an auxiliary momentum and pair it with the position through the Hamiltonian . The joint distribution
factorizes: and are independent, and the -marginal is exactly the target you wanted. So sample the joint by alternating two moves that each leave it invariant — draw a fresh , then run a volume-preserving, energy-conserving Hamiltonian flow — and simply throw the momentum away at the end. What's left is a draw from . The physics was scaffolding; the marginal is the point.
Watch it fly along the ridge
The target is a Gaussian with correlation — the thin diagonal ridge in the demo. Step through HMC and watch each leapfrog trajectory arc the length of the ridge in one proposal; then flip to random walk and watch it jitter in place. Push until the red rejections appear. The sample clouds fill in at wildly different rates — that difference is the whole reason HMC exists.
The core, in code you can run
The entire method is the leapfrog integrator plus one Metropolis test. If the leapfrog looks familiar, it should — it is line-for-line the velocity-Verlet step from the molecular-dynamics page, now walking on instead of a Lennard-Jones potential.
import numpy as np
# Target pi(q) ∝ exp(-U(q)). U is the "potential energy"; grad U is all HMC needs.
def leapfrog(q, p, eps, L, gradU): # this IS velocity-Verlet from MD
p = p - 0.5*eps*gradU(q) # half kick
for i in range(L):
q = q + eps*p # drift (mass matrix M = I)
if i != L-1: p = p - eps*gradU(q) # full kick
p = p - 0.5*eps*gradU(q) # half kick
return q, -p # negate p to keep the proposal reversible
def hmc_step(q, eps, L, U, gradU, rng):
p0 = rng.standard_normal(q.shape) # fresh momentum ~ N(0, I)
qp, pp = leapfrog(q, p0, eps, L, gradU)
H0 = U(q) + 0.5*p0@p0 # energy before
H1 = U(qp) + 0.5*pp@pp # energy after the trajectory
if np.log(rng.random()) < H0 - H1: # Metropolis correction
return qp, True # accept the distant proposal
return q, False # reject, stay putRun it against random-walk Metropolis on the correlated ridge and the gap is not subtle. Both chains are unbiased — they find the same mean and variance — but the effective sample size tells the real story: how many truly independent draws you got out of 20,000.
# Correlated Gaussian target, rho = 0.95 (a thin diagonal ridge), N = 20000
mean var cov_01 accept autocorr τ ESS
HMC (0.00, 0.00) (0.95, 0.95) 0.90 0.96 1.0 20000 (100%)
RW-MH (0.04, 0.04) (1.02, 1.00) 0.96 0.61 127.0 157 (0.8%)
# both unbiased; HMC mixes ~127x better PER SAMPLE on the correlated ridge.
# (exact target moments: mean 0, var 1, cov 0.95)A factor of 127, on a problem with only two dimensions and mild curvature. In the hundred-dimensional posteriors that turn up in Bayesian inference, the gap is the difference between a usable sampler and a hopeless one — which is why HMC (and its auto-tuned descendant, the No-U-Turn Sampler) is the engine inside modern probabilistic programming.
Where this goes
Two honest limits worth knowing. HMC needs the gradient — fine when you have it (often by autodiff), useless when you don't. And it follows continuous valleys, so it struggles to cross the energy barrier between well-separated modes — the same failure the demo's single-mode Gaussian politely hides. The natural next steps: the No-U-Turn Sampler that picks the trajectory length for you, a mass matrix tuned to the target's geometry, and the leap to Riemannian HMC where the metric itself bends with the curvature.