Hamiltonian Monte Carlo

Suppose you want to know the most plausible value of 100 parameters given your data. Bayes's theorem tells you the shape of the answer — the posterior distribution — but computing it requires sampling from a probability density that lives in a 100-dimensional space.

The naive approach, random-walk Metropolis, shuffles its proposal one tiny step at a time. In low dimensions that works. In high dimensions the sampler spends nearly all its time in "corners" where the density is negligible, and the sliver of space that actually carries probability mass is so thin that random steps almost never land on it. Mixing becomes catastrophically slow.

Hamiltonian Monte Carlo (HMC), introduced by Simon Duane et al. in 1987 in the context of lattice field theory and popularized for statistics by Radford Neal and later the Stan team, sidesteps this by treating the negative log-posterior as a potential energy surface. It augments each parameter with a fictitious momentum, then simulates the physical trajectory of a particle rolling across that surface. Because the particle follows the gradient, it glides toward regions of high probability instead of stumbling around randomly — and it can leap across the full width of the posterior in a single proposal.

The mathematics is elegant: Hamiltonian dynamics conserve total energy and are time-reversible, so the Metropolis acceptance step still guarantees the sampler converges to the exact posterior. In practice, HMC accepts proposals at rates above 90 % even in hundreds of dimensions, while random walk would need thousands of tiny steps to travel the same distance.

The demo below draws samples from a banana-shaped posterior — a classic stress-test for samplers because it is highly curved and its width varies dramatically across dimensions. Run both algorithms and compare how many steps each needs to cover the distribution.

<div class="controls">
  <button id="btn-rw" type="button">{{btn_rw}}</button>
  <button id="btn-hmc" type="button">{{btn_hmc}}</button>
  <button id="btn-reset" type="button" class="ghost">{{btn_reset}}</button>
</div>
<div class="legend">
  <span class="dot rw"></span> {{legend_rw}} &nbsp;
  <span class="dot hmc"></span> {{legend_hmc}}
</div>
<canvas id="cv" width="480" height="380"></canvas>
<div id="stats" class="stats"></div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; margin: 0; background: #f8f9fb; }
.controls { display: flex; gap: .5rem; flex-wrap: wrap; margin-bottom: .5rem; }
button { font: 600 13px system-ui; padding: .4rem .85rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 7px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
button:disabled { opacity: .4; cursor: default; }
.legend { font-size: .82rem; color: #444; margin-bottom: .3rem; display: flex; align-items: center; gap: .4rem; }
.dot { display: inline-block; width: 10px; height: 10px; border-radius: 50%; }
.dot.rw  { background: #2196f3; }
.dot.hmc { background: #ff6b35; }
canvas { display: block; border: 1px solid #dde3ea; border-radius: 10px; background: #fff;
         max-width: 100%; }
.stats { margin-top: .45rem; font-size: .82rem; color: #333; min-height: 1.3em; }

// Code not found

Notice how random walk (blue dots) clusters near its starting point for many iterations before slowly diffusing outward, while HMC (orange dots) leaps to distant, high-probability regions immediately. The "banana" shape is deliberately curved so that a sampler must follow the ridge — HMC's gradient information lets it do exactly that; random walk has no such map.

HMC introduces a momentum variable $p$ for each parameter $\theta$, drawn from a Gaussian at the start of each iteration. The joint distribution over $(\theta, p)$ is:

$$H(\theta, p) = U(\theta) + K(p)$$

where $U(\theta) = -\log \pi(\theta \mid \text{data})$ is the potential energy and $K(p) = \tfrac{1}{2} p^\top M^{-1} p$ is the kinetic energy ($M$ is a mass matrix, usually diagonal). Hamilton's equations of motion are:

$$\frac{d\theta}{dt} = \frac{\partial H}{\partial p} = M^{-1} p$$

$$\frac{dp}{dt} = -\frac{\partial H}{\partial \theta} = -\nabla U(\theta)$$

Because these equations preserve total energy $H$ and volume (Liouville's theorem), they define a transition that is:

• Reversible — running the dynamics backward returns to the start.
• Volume-preserving — the Jacobian of the map equals 1.

Both properties are exactly what a Metropolis proposal needs to leave the target distribution invariant.

In practice the equations are integrated with the leapfrog algorithm, a symplectic integrator that keeps energy errors bounded rather than growing, at the cost of a small, correctable bias:

1. Half-step update of momentum: $p \leftarrow p - \tfrac{\varepsilon}{2} \nabla U(\theta)$
2. Full-step update of position: $\theta \leftarrow \theta + \varepsilon M^{-1} p$
3. Half-step update of momentum: $p \leftarrow p - \tfrac{\varepsilon}{2} \nabla U(\theta)$

Repeat for $L$ leapfrog steps, then accept/reject the endpoint $(\theta^*, p^*)$ with probability $\min\!\left(1,\, \exp\!\left(H(\theta,p) - H(\theta^*,p^*)\right)\right)$. Because $H$ is approximately conserved, acceptance rates are very high even for large $L$ — which means each proposal jumps far.

Why does it scale? Random-walk Metropolis needs $O(d^{2})$ steps to explore a d-dimensional Gaussian; HMC needs only $O(d^{5/4})$ — an exponential win in practice. The No-U-Turn Sampler (NUTS), developed by Hoffman and Gelman in 2011, eliminates the need to tune L by automatically stopping the leapfrog trajectory when it doubles back on itself, making HMC practical with no hand-tuning. NUTS is the default in Stan, PyMC, and NumPyro.

See also the related challenge of non-convex optimization, where gradient information similarly separates efficient algorithms from hopeless ones.

HMC is the computational backbone of modern Bayesian data analysis. Its ability to explore high-dimensional posteriors efficiently opens problems that were computationally intractable just a decade ago:

• Clinical trials and epidemiology: hierarchical models with hundreds of random effects — one per patient, school, or hospital — are routinely fit with NUTS in Stan. Estimating COVID-19 transmission rates across 50 US states simultaneously, each with its own dynamics, is a canonical example.
• Astrophysics and cosmology: fitting stellar population models or gravitational-wave parameters involves likelihoods with dozens of correlated parameters and complex, multimodal shapes. The LIGO collaboration uses HMC-based tools to estimate black-hole masses and spins.
• Natural language processing: Bayesian topic models (LDA variants) and neural network uncertainty quantification use HMC when exact posterior geometry matters.
• Drug development: pharmacokinetic/pharmacodynamic (PK/PD) models embed differential equations inside a Bayesian model; HMC is the only practical sampler for the resulting curved, high-dimensional posteriors.
• Ecology: species-distribution models with spatial random fields and occupancy data rely on HMC in INLA and Stan to quantify uncertainty across thousands of grid cells.

The unifying theme is hierarchical structure: whenever parameters share a prior and data is grouped, the posterior is strongly curved and correlated — exactly the setting where HMC's gradient guidance pays off most.

Compare the approach to Bayesian inference, which lays out the framework HMC makes computationally feasible at scale.

Hamiltonian Monte Carlo is a beautiful collision of physics and statistics. By treating probability as a landscape and momentum as a guide, it turns the brutal curse of dimensionality into a navigable slope — the sampler glides where random walk would stumble.

The key insight is not just algorithmic: it says that the geometry of the posterior matters. Problems that seemed like pure data-analysis questions — "fit this model to this data" — turn out to have a physical structure that can be exploited. Gradient information, energy conservation, and time-reversibility are not just mathematical conveniences; they are the reason modern Bayesian science can handle models with thousands of parameters.

NUTS completed the picture by removing the last manual knob, and today HMC in Stan or PyMC is the default choice whenever a posterior is complex enough that random walk would fail. The next time a Bayesian inference problem seems intractable, the answer might be to simulate a particle rolling downhill.