Hamiltonian Monte Carlo

Definition

Hamiltonian Monte Carlo (HMC) is a Markov Chain Monte Carlo (MCMC) method that uses gradient information to make large, directed moves in the posterior, suppressing the random-walk behaviour that makes simple Metropolis-Hastings (MH) slow in high dimensions. It augments the target distribution over parameters $q$ with fictitious momentum variables $p$ , simulates Hamiltonian dynamics jointly over $(q, p)$ via the leapfrog integrator, then accepts or rejects the trajectory endpoint using the Metropolis criterion.

Key Ideas

Augmented state: introduce $p \sim \mathcal{N}(0, M)$ independently of $q$ . Define the Hamiltonian $H(q,p) = U(q) + K(p)$ where $U(q) = -\log\pi(q)$ is the potential energy and $K(p) = p^T M^{-1} p / 2$ is kinetic energy. The joint density $\exp(-H) = \pi(q) \times \mathcal{N}(p; 0, M)$ , so the marginal over $p$ is the target.
One HMC step: (1) resample $p \sim \mathcal{N}(0, M)$ ; (2) simulate $L$ leapfrog steps from $(q, p)$ with stepsize $\varepsilon$ to produce proposal $(q^*, p^*)$ ; (3) accept with probability $\min(1, \exp(H(q,p) - H(q^*,p^*)))$ . Negate $p^*$ before presenting the proposal to restore detailed balance.
Three key properties of Hamiltonian dynamics:
- Reversibility: negating momentum reverses the trajectory, making the proposal time-reversible.
- Volume preservation (Liouville's theorem): the leapfrog map is symplectic — it exactly preserves phase-space volume, so the Jacobian is 1 and no correction is needed in the acceptance probability.
- Approximate energy conservation: $H$ changes only due to leapfrog discretization error. Because $H \approx \text{const}$ , proposals are distant in $q$ but have high acceptance probability.

How It Works

Leapfrog Integrator

The Störmer-Verlet / leapfrog scheme interleaves half-steps for momentum with full steps for position, for i = 1, ..., L:

$p(t + \varepsilon/2) = p(t) - (\varepsilon/2)\,\nabla U(q(t))$
$q(t + \varepsilon) = q(t) + \varepsilon M^{-1} p(t + \varepsilon/2)$
$p(t + \varepsilon) = p(t + \varepsilon/2) - (\varepsilon/2)\,\nabla U(q(t + \varepsilon))$

This is the unique first-order symplectic integrator for separable Hamiltonians; it exactly preserves a "shadow" Hamiltonian $\tilde{H} \approx H + O(\varepsilon^2)$ , ensuring acceptance rates remain high even with moderately large stepsizes.

Tuning

Two hyperparameters: stepsize $\varepsilon$ and number of leapfrog steps $L$ (trajectory length $\varepsilon L$ ). Optimal acceptance rate ≈ 65% (Roberts-Gelman-Gilks 1997, analogous to the 0.234 rule for random-walk Metropolis in $d$ dimensions). Trajectory length should be long enough for q* to be nearly independent of q, but not so long that the trajectory loops back. The No-U-Turn Sampler (NUTS, Hoffman-Gelman 2011) automates both: it grows the trajectory until the simulation turns around (U-turns), eliminating manual tuning. NUTS is the default in Stan.

Mass Matrix

$M$ acts as a preconditioning matrix. Setting $M = \Sigma^{-1}$ (inverse posterior covariance) decorrelates parameters and equalises leapfrog stepsizes across dimensions. In practice $M$ is estimated from warmup iterations (diagonal approximation is standard).

Why It Matters

Scaling: the total arithmetic cost per independent sample scales as $O(d^2)$ for random-walk Metropolis versus $O(d^{5/4})$ for HMC — a $d^{3/4}$ improvement. (Equivalently: HMC needs $O(d^{1/4})$ gradient evaluations per sample and RWM needs $O(d)$ density evaluations, each evaluation itself costing $O(d)$ arithmetic.) In high-dimensional posteriors (modern hierarchical models, neural networks, latent Gaussian fields) this is a decisive advantage.

Probabilistic programming: HMC/NUTS is the sampler in Stan, PyMC, and Turing.jl. The ability to target arbitrary differentiable posteriors without hand-crafted proposal distributions enabled a generation of applied Bayesian models that would have been computationally infeasible with Gibbs or random-walk MH.

Limitations:

Requires continuous, differentiable $U(q) = -\log\pi(q)$ . Cannot directly handle discrete parameters or non-differentiable densities.
Ineffective for highly multimodal distributions where modes are separated by very low-probability regions; tempered transitions or parallel tempering must be combined with HMC for such targets.
The gradient $\nabla U$ must be computed at each leapfrog step, which is expensive when the likelihood involves large datasets (though stochastic gradient HMC variants exist).

Open Questions

Optimal trajectory length and stepsize adaptation in non-stationary warm-up phases.
Handling of mixed discrete-continuous posteriors (relaxation / marginalization strategies).
Validity of subsampled/stochastic gradient HMC (Welling-Teh 2011 Stochastic Gradient Langevin Dynamics (SGLD), Chen et al. 2014 Stochastic Gradient HMC (SGHMC)): detailed balance is broken; correctness requires vanishing step-size schedules or correction terms.
Riemannian HMC (Girolami-Calderhead 2011) uses a position-dependent mass matrix $M(q)$ equal to the Fisher information; dramatically improves geometry adaptation but requires costly tensor computations.