Hamiltonian Monte Carlo: How Momentum Turns Random Walks into Efficient Sampling

Most Bayesian models are easy to write down and hard to sample from. You can usually evaluate an unnormalised posterior density, or its gradient, but you cannot draw independent samples from it directly. Markov chain Monte Carlo (MCMC) solves this by constructing a chain whose long-run distribution is the target. The question is how efficiently that chain explores.

Hamiltonian Monte Carlo (HMC) is one of the most effective answers for continuous, differentiable models. Instead of proposing a small random step and hoping it lands somewhere useful, HMC borrows an idea from physics: give the current point momentum, let it glide along the contours of the target distribution, and only then decide whether to accept the move. The result is a sampler that can travel far in one step while still keeping a high acceptance rate.

The problem with random walks

Consider a standard random-walk Metropolis sampler. From the current state, it proposes a new state by adding Gaussian noise, then accepts or rejects that proposal based on the density ratio. This works, but it has a structural weakness: to keep the acceptance rate reasonable, the step size must stay small, and small steps mean the chain diffuses slowly.

In high dimensions this gets worse. A random walk needs on the order of the square of the distance it must travel to cover that distance, because successive steps partly cancel each other out. If the target is strongly correlated, for example a narrow ridge in parameter space, an isotropic proposal will either step off the ridge and be rejected, or crawl along it very slowly.

HMC avoids this by using gradient information to propose coherent, long-distance moves.

A physical analogy that is more than an analogy

Imagine a particle sliding without friction over a landscape. Valleys correspond to regions of high posterior probability, and hills correspond to regions of low probability. If you give the particle a push, it will roll down into a valley, climb partway up the other side, turn around, and keep moving in a way that respects the shape of the landscape.

HMC makes this literal. Let

$q$

denote the position, which is the parameter vector we want to sample, and let

$p$

denote an auxiliary momentum variable of the same dimension. Define a potential energy from the target density

$\pi(q)$

as

$$U(q) = -\log \pi(q)$$

so that high probability means low potential energy. Define a kinetic energy as

$$K(p) = \frac{1}{2} p^{T} M^{-1} p$$

where

$M$

is a symmetric positive definite mass matrix, often taken as the identity or as a diagonal estimate of the posterior scales. The total energy, or Hamiltonian, is

$$H(q, p) = U(q) + K(p)$$

The dynamics are given by Hamilton's equations:

$$\frac{dq}{dt} = \frac{\partial H}{\partial p} = M^{-1} p, \quad \frac{dp}{dt} = -\frac{\partial H}{\partial q} = -\nabla U(q) = \nabla \log \pi(q)$$

Two properties matter. First, the dynamics conserve total energy, so a trajectory stays on a level set of

$H$

Second, the dynamics preserve volume in

$(q, p)$

space and are reversible in time. Those properties are what allow HMC to propose distant states and still correct exactly for numerical error with a Metropolis step.

The HMC algorithm

One HMC iteration works as follows:

1. Sample fresh momentum from its Gaussian distribution, independently of the current position. With the kinetic energy above, this is

$$p \sim \mathcal{N}(0, M)$$

1. Simulate Hamiltonian dynamics for a fixed trajectory length, starting from the current

$(q, p)$

pair. This produces a proposed pair

$(q^{}, p^{})$

1. Accept the proposal with probability

$$\min\left(1, \exp\left(H(q, p) - H(q^{}, p^{})\right)\right)$$

and otherwise stay at the current position. Finally, discard the momentum: the marginal distribution of

$q$

under the joint distribution is exactly the target, because the joint density is proportional to

$\exp(-H(q, p))$

and integrating out a Gaussian momentum leaves

$\pi(q)$

If the dynamics were simulated exactly, energy would be conserved and every proposal would be accepted. In practice we must discretise, and the Metropolis correction removes the bias introduced by that discretisation.

Why the leapfrog integrator matters

HMC almost always uses the leapfrog integrator, not a generic numerical solver, because leapfrog preserves the two properties that make the correction work: it is volume-preserving and time-reversible. A single leapfrog step of size

$\epsilon$

is:

$$p_{t + \epsilon/2} = p_{t} - \frac{\epsilon}{2} \nabla U(q_{t})$$

$$q_{t + \epsilon} = q_{t} + \epsilon M^{-1} p_{t + \epsilon/2}$$

$$p_{t + \epsilon} = p_{t + \epsilon/2} - \frac{\epsilon}{2} \nabla U(q_{t + \epsilon})$$

Repeat this for

$L$

steps to get a trajectory of total length

$L \epsilon$

Leapfrog is symplectic, meaning it approximately preserves the geometric structure of Hamiltonian dynamics, so energy error stays bounded over long trajectories instead of drifting steadily. A Runge-Kutta method, by contrast, is generally not volume-preserving, and using it would break the simple Metropolis correction.

The step size controls accuracy: too large, and energy error grows and proposals are rejected; too small, and computation is wasted on many tiny steps. The trajectory length controls how far the chain travels: too short, and HMC degenerates back toward a random walk; too long, and the trajectory may loop back near its starting point, wasting gradient evaluations. A practical rule is to randomise the trajectory length slightly on each iteration, which avoids periodic behaviour where a fixed length repeatedly returns to almost the same place.

A concrete example: a correlated Gaussian

The advantage of HMC is easiest to see on a strongly correlated two-dimensional Gaussian. A random-walk proposal aligned with the coordinate axes must take tiny steps to stay on the narrow ridge. HMC, by contrast, follows the gradient, so its trajectory curves along the ridge itself.

A minimal leapfrog implementation in NumPy makes the structure clear:

import numpy as np

def leapfrog(q, p, grad_log_pi, epsilon, n_steps):
    # grad_log_pi returns the gradient of log pi at q
    p = p + 0.5 * epsilon * grad_log_pi(q)
    for i in range(n_steps):
        q = q + epsilon * p
        if i != n_steps - 1:
            p = p + epsilon * grad_log_pi(q)
    p = p + 0.5 * epsilon * grad_log_pi(q)
    # Negate momentum for reversibility; harmless for a symmetric proposal
    return q, -p

def hmc_step(q, log_pi, grad_log_pi, epsilon, n_steps):
    p0 = np.random.randn(*q.shape)
    q_new, p_new = leapfrog(q.copy(), p0, grad_log_pi, epsilon, n_steps)
    log_accept = (
        log_pi(q_new) - 0.5 * np.dot(p_new, p_new)
        - (log_pi(q) - 0.5 * np.dot(p0, p0))
    )
    if np.log(np.random.rand()) < log_accept:
        return q_new, True
    return q, False

Here the mass matrix is the identity, so the kinetic energy is simply half the squared norm of the momentum. In a real model, you would replace the identity with a diagonal mass matrix estimated from warmup samples, which rescales coordinates so that a single step size works reasonably well in every direction.

For a Bayesian logistic regression, the same code pattern applies with the log posterior and its gradient in place of the toy density. The gradient is the expensive part, but it is also the reason each iteration can move much further than a gradient-free proposal.

Tuning, warmup, and the No-U-Turn Sampler

In practice, nobody hand-tunes

$\epsilon$

and

$L$

for every model. Modern implementations split sampling into a warmup phase and a sampling phase. During warmup, the step size is adapted to hit a target acceptance rate, often around 0.8, and the mass matrix is estimated from the samples seen so far. Those adapted values are then frozen for the sampling phase, so the chain is a valid Markov chain with a fixed transition kernel.

The No-U-Turn Sampler (NUTS), used by default in Stan and PyMC, goes further: it automatically chooses a trajectory length by simulating both forward and backward in time until the trajectory starts to turn back on itself, the "U-turn" condition. This removes the most fragile tuning choice in HMC, at the cost of a more complex algorithm and a variable number of gradient evaluations per iteration. If you use a probabilistic programming library, you are almost certainly using NUTS rather than hand-tuned HMC, but the underlying dynamics are the same.

Diagnostics: what to watch

HMC gives you diagnostic signals that a random-walk sampler does not. The most important is the divergent transition. A divergence occurs when the leapfrog integrator encounters a region of very high curvature, for example the neck of a funnel-shaped posterior, and the energy error becomes large. Divergences are not just rejected proposals: they indicate that the sampler cannot reliably explore that region, and posterior estimates from a chain with many divergences should not be trusted. The usual remedies are a smaller step size, a higher target acceptance rate during warmup, or a reparameterisation of the model.

Beyond divergences, the usual MCMC checks still apply: split-R-hat near 1 across independent chains, effective sample size large relative to the number of draws, and trace plots without long flat stretches or trends. HMC typically achieves a much higher effective sample size per iteration than random-walk Metropolis on smooth, continuous problems, which is why the extra gradient cost is usually worth paying.

Limitations and common pitfalls

HMC is not a universal sampler. It requires a differentiable log density, so it does not apply directly to discrete parameters. A common workaround is to marginalise discrete parameters out when possible, or to use a different sampler for the discrete part in a Gibbs-style scheme.

Multimodal targets are another weakness. Hamiltonian trajectories conserve energy, so a particle starting in one mode with a typical momentum draw will rarely gain enough energy to cross a low-density valley into another mode. Tempering, or a sampler designed for mode jumping, is a better fit there.

Three pitfalls recur in practice. First, forgetting that the gradient must be correct: a wrong sign, a missing Jacobian adjustment after a transformation, or a gradient computed at the wrong point will silently produce a wrong stationary distribution, because the Metropolis correction only corrects integrator error, not gradient error. Second, adapting the step size or mass matrix throughout sampling rather than only during warmup, which breaks the Markov property. Third, interpreting a high acceptance rate as success on its own: a chain with tiny steps accepts almost everything and explores almost nothing. Acceptance rate, effective sample size, and divergences must be read together.

Takeaway

Random-walk Metropolis explores by diffusion, which is slow in high dimensions and on correlated targets. Hamiltonian Monte Carlo explores by dynamics: momentum carries the chain coherently along the contours of the posterior, the leapfrog integrator preserves the geometry that makes the correction exact, and a Metropolis step removes the remaining numerical bias. That combination is why HMC, and its automatic cousin NUTS, became the default engine for modern Bayesian computation on continuous models. It asks for more per iteration, a gradient and a careful warmup, and in return it gives back far more per iteration in effective exploration.

---

Connect with me on LinkedIn: linkedin.com/in/umeshtharukaofficial