Hamiltonian Monte Carlo

Hamiltonian Monte Carlo

Hamiltonian Monte Carlo (HMC) is a gradient-informed Metropolis-Hastings sampler whose proposal follows a continuous-time physical trajectory determined by the target $\pi$. Augment the state $x\in {\mathbb{R}}^d$ with an auxiliary momentum variable $p\in {\mathbb{R}}^d$ and define the Hamiltonian
$H(x,p)=U(x)+K(p)=-\log \pi (x)+\frac{1}{2}{p}^{\top }{M}^{-1}p,$
where $U(x)=-\log \pi (x)$ is the potential energy, $K(p)$ is the kinetic energy, and $M$ is a positive-definite mass matrix (default $M=I_d$). The joint target on $(x,p)$ is $\pi (x,p)\propto \exp (-H(x,p))=\pi (x)\cdot \mathcal{N}(p;0,M)$, so the $x$-marginal is exactly $\pi$.

One iteration:

1. Resample momentum. Draw $p\sim \mathcal{N}(0,M)$, independent of the current position $x$.
2. Simulate Hamiltonian dynamics. Integrate
   $\dot{x}={\nabla }_{p}H=M^{-1}p,\dot{p}=-{\nabla }_{x}H=\nabla \log \pi (x)$
   for time $T=L\epsilon$ using $L$ leapfrog steps of size $\epsilon$:
   $p\leftarrow p+\frac{\epsilon }{2}\nabla \log \pi (x),x\leftarrow x+\epsilon M^{-1}p,p\leftarrow p+\frac{\epsilon }{2}\nabla \log \pi (x).$
3. Metropolis accept/reject. Let $(x^{\prime },p^{\prime })$ be the leapfrog endpoint. Accept with probability
   $\alpha =\min (1,\exp (H(x,p)-H(x^{\prime },p^{\prime }))).$

Why it works. Continuous-time Hamiltonian dynamics exactly preserve $H$, so acceptance would be identically $1$ if the integrator were exact. The leapfrog integrator is a symplectic scheme — it preserves phase-space volume exactly and conserves $H$ to second order in $\epsilon$, so the acceptance probability is nearly $1$ even for large numbers of steps $L$. Symplecticity plus time-reversibility of leapfrog (with a momentum flip) is exactly what makes the Metropolis step satisfy detailed balance with respect to $\pi (x,p)$.

Mechanical interpretation. HMC is the position-dependent acceleration extension of MALA. MALA evaluates $\nabla \log \pi$ once at the starting position $x$ and holds it fixed for the flight time $\sigma$; HMC updates the acceleration $\nabla \log \pi (x(t))$ continuously along the trajectory, so the particle “skates” along contours of $H$ rather than falling along a single gradient direction. The particle can travel long distances before Metropolis-checking, so HMC takes much larger effective steps than MALA.

Cost and scaling. One HMC iteration costs $L$ gradient evaluations ($L$ typically $10$–$100$). In exchange, HMC mixes with step size $\epsilon \propto d^{-1/4}$ (versus $d^{-1/6}$ for MALA and $d^{-1/2}$ for random-walk Metropolis), and the optimal acceptance rate is roughly $0.65$. For smooth high-dimensional targets, HMC’s effective-sample-size-per-gradient is dramatically higher than MALA’s, which is why it is the default sampler in modern probabilistic programming systems like Stan, PyMC, and NumPyro. The No-U-Turn Sampler (NUTS) adaptively chooses $L$ to avoid hand-tuning the trajectory length.

Prerequisites

• 3.11.3 MALA — the constant-acceleration predecessor; HMC generalizes to position-dependent acceleration.
• 3.11.1 Metropolis-Hastings — HMC is an MH instance whose proposal is the leapfrog endpoint and whose acceptance ratio uses the Hamiltonian.
• 3.11 MCMC — detailed balance, ergodic theorem, effective sample size.
• 5.4 Reversibility — leapfrog with momentum flip is time-reversible, which yields detailed balance.

Intuition

Imagine a frictionless puck sliding across a surface whose height is $U(x)=-\log \pi (x)$. Low-posterior regions correspond to tall hills; high-posterior regions are valleys. Give the puck a random push and let it coast: it will speed up sliding into valleys and slow climbing hills, so it naturally spends more time in valleys (high-density regions). If we track its position periodically and record samples, the histogram of positions is $\pi$. This is exactly what HMC simulates, with kinetic energy $K(p)=\frac{1}{2}{p}^{\top }{M}^{-1}p$ playing the role of the puck’s kinetic energy and momentum resampling replacing occasional “shakes” that keep the puck exploring all valleys rather than settling in the first one.

The central insight that makes HMC fast is that the puck’s trajectory stays on a level set of $H$, and level sets stretch across the distribution — a single simulation can traverse a long distance without ever climbing a hill. Compare this to MALA, which takes one gradient step and then pauses for a Metropolis check: MALA’s effective step size is $\sigma \sim d^{-1/6}$, while HMC’s effective distance traveled per iteration scales with the trajectory length $L\epsilon$, which can be much larger. The price is that Hamiltonian dynamics have no closed-form solution for general $\pi$ and must be discretized; leapfrog is the numerical scheme that preserves both the structure needed for the Metropolis step (time-reversibility, volume preservation) and the quality needed for high acceptance ($H$ conserved to $\mathcal{O}({\epsilon }^2)$).

Properties

(Position-Dependent Acceleration) The acceleration $\nabla \log \pi (x(t))$ is updated at every leapfrog step, unlike MALA.

The MALA proposal uses a single gradient evaluation at the starting point $x$ and pretends the acceleration stays constant for time $\sigma$; this is accurate only for small $\sigma$. HMC evaluates $\nabla \log \pi$ at every leapfrog half-step, so the trajectory bends along the contours of $\pi$. The payoff is that HMC can take many small leapfrog steps and cover a long total distance with the acceleration always pointing in the right direction — the chain tracks the geometry of $\pi$ rather than fighting it.

(Symplectic Leapfrog Integrator) Leapfrog preserves phase-space volume exactly and conserves $H$ to second order.

The leapfrog update is the canonical symplectic integrator: the Jacobian of one leapfrog step as a map on $(x,p)$ has determinant exactly $1$, so the measure $dxdp$ is preserved. This volume-preservation is necessary for the Metropolis step to have a simple acceptance formula (no Jacobian correction). Additionally, the numerical error in $H$ is $\mathcal{O}({\epsilon }^2)$ per step and does not accumulate — leapfrog nearly-conserves a shadow Hamiltonian $\tilde{H}=H+\mathcal{O}({\epsilon }^2)$, so even over long trajectories ($L\epsilon =\mathcal{O}(1)$) the energy error stays bounded.

(Time-Reversibility with Momentum Flip) The leapfrog map is reversed by a momentum flip $p\mapsto -p$.

If $(x^{\prime },p^{\prime })={\text{Leapfrog}}_L(x,p)$, then ${\text{Leapfrog}}_L(x^{\prime },-p^{\prime })=(x,-p)$. Combined with symplecticity, this reversibility is exactly the symmetry needed so that the MH acceptance probability $\min (1,e^{H(x,p)-H(x^{\prime },p^{\prime })})$ satisfies detailed balance with respect to $\pi (x,p)$ — the proposal is effectively symmetric on the $(x,p)$ augmented space once the momentum flip is absorbed into the Metropolis step. (Most implementations omit the explicit momentum flip because momentum is resampled at the start of every iteration anyway.)

(Metropolis Acceptance) The accept/reject step uses the change in $H$, with acceptance probability $\alpha =\min (1,e^{-\Delta H})$.

If leapfrog were exact, $H$ would be conserved and $\alpha \equiv 1$. Discretization makes $\Delta H=H(x^{\prime },p^{\prime })-H(x,p)=\mathcal{O}({\epsilon }^2)$ per leapfrog step, but because the shadow-Hamiltonian error does not accumulate, the total $\Delta H$ over $L$ steps is also $\mathcal{O}({\epsilon }^2)$. This gives acceptance probabilities in the $0.6$–$0.9$ range for typical step sizes, compared to $\sim 0.25$ for random-walk Metropolis at the same per-iteration cost.

(Optimal Acceptance ≈ 0.65 and $\epsilon \propto d^{-1/4}$) HMC scales better with dimension than MALA.

For i.i.d. product targets $\pi (\theta )={\prod }_{i=1}^{d}f({\theta }_i)$, theory (Beskos, Pillai, Roberts, Sanz-Serna, Stuart 2013) shows the optimal step size scales as $\epsilon \propto d^{-1/4}$ with optimal acceptance rate $\approx 0.65$. This is strictly better than MALA’s $d^{-1/6}$ scaling with $0.57$ acceptance and dramatically better than RWM’s $d^{-1/2}$ scaling with $0.23$ acceptance. The practical consequence is that HMC remains efficient at dimensions ($d\sim 10^3$–$10^5$) where MALA and RWM become prohibitively slow.

(Momentum Resampling Ensures Ergodicity) Resetting $p\sim \mathcal{N}(0,M)$ at the start of every iteration is what makes the chain explore the full space.

A single Hamiltonian trajectory is confined to the level set $\{(x,p):H(x,p)=H_0\}$ — it cannot reach regions of $\pi$ at different energies. Resampling $p$ from its marginal Gaussian changes the Hamiltonian level and lets the chain move between energy shells; after many iterations, the full $(x,p)$ distribution is explored, and the $x$-marginal is $\pi$. Without momentum resampling, HMC is not ergodic.

(Mass Matrix Preconditions the Dynamics) A non-identity mass matrix $M$ rescales the proposal to match the local geometry of $\pi$.

The default $M=I_d$ is the wrong shape when coordinates of $\theta$ live on very different scales. A posterior covariance estimate $\hat{\Sigma }$ (from a preliminary run or a Laplace approximation at the MAP) used as $M={\hat{\Sigma }}^{-1}$ equalizes the momentum-induced steps across dimensions. Modern HMC adaptations (warmup phase in Stan, NumPyro) estimate $M$ automatically during the tuning phase by running several HMC iterations, computing the sample covariance of the positions, and using its inverse as the mass matrix for the remaining iterations.

(No-U-Turn Sampler (NUTS)) Automatically chooses trajectory length $L$ to avoid backtracking.

Choosing $L$ too small gives short HMC trajectories that behave like MALA; choosing $L$ too large wastes computation on trajectories that loop back to their starting point (called a “U-turn”). The No-U-Turn Sampler (Hoffman & Gelman 2014) detects U-turns on the fly by growing the trajectory doubling each iteration and stopping when the momentum would reverse. NUTS is the default in Stan and PyMC because it removes the most delicate HMC hyperparameter ($L$) from the user’s hands.

(Requires Smooth, Gradient-Accessible Targets) HMC breaks on discrete parameters, discontinuities, and heavy tails.

Like MALA, HMC needs $\nabla \log \pi$ to be defined and computable. For discrete parameters, HMC cannot be used directly — collapsed Gibbs or block updates must handle them while HMC handles the continuous parameters. For non-differentiable log-posteriors (e.g., Laplace priors with exact zero, hinge-loss-like terms), HMC behaves erratically at the kinks; reparameterization or smoothing is typically needed. Heavy-tailed targets make the Hamiltonian dynamics unstable at large $\Vert x\Vert$ because the gradient grows without bound, and leapfrog energy errors compound — Riemannian HMC with position-dependent metrics partially addresses this but is rarely used in practice.

Examples

HMC on a 2D Gaussian

Target $\pi (x)=\mathcal{N}(x;0,\Sigma )$ with $\Sigma =\left(\begin{array}{cc}1& 0.9\\ 0.9& 1\end{array}\right)$. The potential and its gradient are
$U(x)=\frac{1}{2}{x}^{\top }{\Sigma }^{-1}x,\nabla U(x)={\Sigma }^{-1}x.$
With mass matrix $M=I_2$ and step size $\epsilon =0.1$, run $L=20$ leapfrog steps per iteration. Each leapfrog step is:

• $p\leftarrow p-\frac{\epsilon }{2}{\Sigma }^{-1}x$
• $x\leftarrow x+\epsilon p$
• $p\leftarrow p-\frac{\epsilon }{2}{\Sigma }^{-1}x$

The trajectory traces an ellipse in $(x,p)$ space determined by $H=\frac{1}{2}{x}^{\top }{\Sigma }^{-1}x+\frac{1}{2}\Vert p{\Vert }^2$. Over $N=5,000$ iterations the effective sample size for each coordinate is typically around $4,000$ — in contrast, random-walk Metropolis gives $\sim 500$ and MALA $\sim 2,000$ on the same target. Setting $M={\Sigma }^{-1}$ (the ideal preconditioner) produces exact independent samples after a single iteration because the Hamiltonian becomes rotationally symmetric.

Why leapfrog and not Euler?

The Euler integrator
$x\leftarrow x+\epsilon p,p\leftarrow p+\epsilon \nabla \log \pi (x)$
is not symplectic: the Jacobian determinant is $1-{\epsilon }^2(\text{stuff})$, so volume is not preserved. Worse, energy errors accumulate linearly: $\Delta H$ grows like $L{\epsilon }^2$, so long trajectories are rejected with probability approaching $1$. Leapfrog
$p\leftarrow p+\frac{\epsilon }{2}\nabla \log \pi (x),x\leftarrow x+\epsilon p,p\leftarrow p+\frac{\epsilon }{2}\nabla \log \pi (x)$
is symplectic (volume-preserving) and keeps the energy error bounded near the shadow Hamiltonian. Running $L=20$ Euler steps vs. $L=20$ leapfrog steps on the 2D Gaussian above gives acceptance rates near $0$ for Euler and near $0.95$ for leapfrog with the same step size. Symplectic integration is not a minor implementation detail — it is what makes HMC work.

Hierarchical model with funnel geometry

The Bayesian hierarchical model
$\tau \sim \text{HalfNormal}(0,1),{\theta }_i\mid \tau \stackrel{\text{iid}}{\sim }\mathcal{N}(0,{\tau }^2),y_i\mid {\theta }_i\sim \mathcal{N}({\theta }_i,{\sigma }_y^2)$
produces a “Neal’s funnel” posterior: as $\tau \to 0$, the ${\theta }_i$ concentrate near $0$; as $\tau$ grows, the ${\theta }_i$ spread out. HMC with a fixed mass matrix gets stuck near the funnel neck because the local curvature changes dramatically with $\tau$. The standard fix is a non-centered reparameterization: ${\theta }_i=\tau \cdot {\tilde{\theta }}_i$ with ${\tilde{\theta }}_i\sim \mathcal{N}(0,1)$ a priori. This decouples $\tau$ from $\tilde{\theta }$ and makes the posterior approximately isotropic, at which point vanilla HMC mixes well. Reparameterization is a common HMC workflow step: the sampler exposes geometry problems that are fixed not by tuning HMC but by rewriting the model.

Derivation

Goal. Show that the HMC proposal satisfies detailed balance with respect to the augmented target $\pi (x,p)=\pi (x)\cdot \mathcal{N}(p;0,M)$, hence its $x$-marginal is $\pi$.

Step 1 — Target augmentation. Define the joint target on $(x,p)\in {\mathbb{R}}^{2d}$ by
$\pi (x,p)\propto \exp (-H(x,p)),H(x,p)=-\log \pi (x)+\frac{1}{2}{p}^{\top }{M}^{-1}p.$
Because $H$ separates additively, $x$ and $p$ are independent under $\pi (x,p)$: $p\sim \mathcal{N}(0,M)$ marginally, and the $x$-marginal is exactly $\pi$. Any sampler that preserves $\pi (x,p)$ as its stationary distribution thus provides samples from $\pi$ by projection.

Step 2 — Momentum resampling preserves $\pi (x,p)$. At the start of each iteration, replace the current momentum with a fresh draw $p\sim \mathcal{N}(0,M)$, independent of $x$. Because $\pi (x,p)$ factors and the $p$-conditional is exactly $\mathcal{N}(0,M)$, this Gibbs-style update leaves $\pi (x,p)$ invariant. (In MH language, this is a full-conditional step with acceptance probability $1$.)

Step 3 — Exact Hamiltonian dynamics preserve $H$ and hence $\pi (x,p)$. Hamilton’s equations
$\dot{x}={\nabla }_{p}H=M^{-1}p,\dot{p}=-{\nabla }_{x}H=\nabla \log \pi (x)$
conserve $H$ (classical mechanics: $\frac{dH}{dt}={\nabla }_{x}H\cdot \dot{x}+{\nabla }_{p}H\cdot \dot{p}={\nabla }_{x}H\cdot M^{-1}p-{\nabla }_{p}H\cdot {\nabla }_{x}H=0$ by substitution) and are volume-preserving by Liouville’s theorem. Under an exact Hamiltonian flow, if $(x_0,p_0)\sim \pi (x,p)$ then $(x_T,p_T)\sim \pi (x,p)$ as well — the dynamics leave the joint target invariant.

Step 4 — Leapfrog discretization is volume-preserving and time-reversible. The leapfrog map ${\Phi }_{\epsilon }:(x,p)\mapsto (x^{\prime },p^{\prime })$ given by
$p\leftarrow p+\frac{\epsilon }{2}\nabla \log \pi (x),x^{\prime }=x+\epsilon M^{-1}p,p^{\prime }=p+\frac{\epsilon }{2}\nabla \log \pi (x^{\prime })$
has Jacobian determinant exactly $1$ (direct computation: each half-step has a triangular Jacobian with unit diagonal, so the determinant is $1$). It is also time-reversible in the sense ${\Phi }_{\epsilon }^{-1}=F\circ {\Phi }_{\epsilon }\circ F$ where $F:(x,p)\mapsto (x,-p)$ flips momentum. These two properties — symplecticity and reversibility — are exactly what is needed for the Metropolis step to satisfy detailed balance.

Step 5 — Metropolis correction gives exact detailed balance. Let $(x^{\prime },p^{\prime })$ be the endpoint of $L$ leapfrog steps, and accept with probability $\alpha =\min (1,e^{H(x,p)-H(x^{\prime },p^{\prime })})$. The joint kernel is
$P((x,p),(x^{\prime },p^{\prime }))=Q((x,p),(x^{\prime },p^{\prime }))\alpha ((x,p),(x^{\prime },p^{\prime })),$
where $Q$ is the deterministic leapfrog transition. Because leapfrog is volume-preserving, the Jacobian factor in the density cancels; because it is reversible with a momentum flip, the forward and backward deterministic densities are linked by
$Q((x,p),(x^{\prime },p^{\prime }))=Q((x^{\prime },-p^{\prime }),(x,-p)).$
Combined with the symmetry $\mathcal{N}(p;0,M)=\mathcal{N}(-p;0,M)$, the acceptance rule makes $\pi (x,p)P((x,p),(x^{\prime },p^{\prime }))=\pi (x^{\prime },p^{\prime })P((x^{\prime },p^{\prime }),(x,p))$ — detailed balance on the augmented state space. Marginalizing over $p$ recovers samples from $\pi$.

Why leapfrog’s energy error stays small. Although ${\Phi }_{\epsilon }$ does not exactly preserve $H$, symplectic integrators preserve a nearby shadow Hamiltonian $\tilde{H}=H+\mathcal{O}({\epsilon }^2)$ exactly. Over $L$ steps the true $H$ error stays bounded by $|\tilde{H}-H|=\mathcal{O}({\epsilon }^2)$ rather than accumulating to $\mathcal{O}(L{\epsilon }^2)$, which is the specific reason HMC can afford long trajectories ($L\epsilon =\mathcal{O}(1)$) with acceptance near $1$.

See Also

• 3.11.3 MALA — constant-acceleration predecessor; HMC generalizes to position-dependent acceleration.
• 3.11.1 Metropolis-Hastings — HMC is an MH sampler on the augmented $(x,p)$ state space.
• 3.11.2 Gibbs Sampling — alternative when full conditionals are tractable or when parameters are discrete.
• 3.11 MCMC — parent index with tuning, diagnostics, effective sample size.
• 3.9 Laplace Approximation — source of the mass matrix $M$ for preconditioning.
• 4.6 Brownian Motion — contrast: MALA’s SDE has a Brownian driving term; HMC’s ODE is deterministic between momentum resets.