Leapfrog Integration in Hamiltonian Monte Carlo

Why the Leapfrog Map Matters

In Hamiltonian Monte Carlo, the posterior density determines a potential energy, while an auxiliary momentum variable supplies kinetic energy. The resulting Hamiltonian equations describe a trajectory through parameter space. Leapfrog integration approximates that trajectory using a finite sequence of position and momentum updates. Its importance is not simply that it approximates a differential equation: its reversibility and volume preservation make it suitable for constructing valid HMC proposals, with a Metropolis correction for the energy error.

For an NLME posterior, the position vector \(q\) can contain unconstrained coordinates representing population parameters, random effects, or other unknowns. The potential energy is the negative log posterior, up to an additive constant. The momentum \(p\) is an auxiliary variable, not a proportion. This tutorial focuses on the numerical map that advances \((q,p)\); it assumes that the potential and its gradient can be evaluated at the positions encountered.

Deriving the Updates From Hamilton’s Equations

Hamilton’s equations for this Hamiltonian are $$ \frac{dq}{dt}=M^{-1}p,\qquad \frac{dp}{dt}=-\nabla U(q). $$ The position derivative depends on momentum, and the momentum derivative depends on position. Leapfrog is obtained by splitting these equations into two simpler systems and composing their exact solutions.

First consider a momentum-only substep: hold \(q\) fixed and evolve \(dp/dt=-\nabla U(q)\) for time \(a\). Since the right-hand side is constant during this substep, its exact update is \(p\leftarrow p-a\nabla U(q)\). Next consider a position-only substep: hold \(p\) fixed and evolve \(dq/dt=M^{-1}p\) for time \(a\). Its exact update is \(q\leftarrow q+aM^{-1}p\). Applying a momentum substep of duration \(\epsilon/2\), a position substep of duration \(\epsilon\), and another momentum substep of duration \(\epsilon/2\) gives the stated leapfrog equations.

The order matters. The first half momentum update uses the gradient at the old position \(q\); the last half update uses the gradient at the new position \(q'\). Using the old gradient twice, or taking a full momentum update before and after the position update, gives a different numerical method. For a fixed \(M\), the equations also show precisely where parameter scaling enters: position changes in the direction \(M^{-1}p_{1/2}\).

Reversibility and Volume Preservation

A numerical integrator need not reproduce the exact Hamiltonian trajectory. For HMC, two structural properties of the leapfrog map are especially useful: it can be reversed by reversing momentum and integrating forward again, and it preserves volume in position-momentum space. These properties hold when the updates are well-defined, \(M\) is fixed and positive definite, and the potential is sufficiently differentiable for the map and its derivatives to exist.

To see reversibility, define \(K_a\) as a momentum update of duration \(a\), \(K_a(q,p)=(q,p-a\nabla U(q))\), and \(D_a\) as a position update of duration \(a\), \(D_a(q,p)=(q+aM^{-1}p,p)\). One leapfrog step is the composition $$ \Phi_\epsilon=K_{\epsilon/2}\circ D_\epsilon\circ K_{\epsilon/2}. $$ Each substep is inverted by changing the sign of its duration: \(K_a^{-1}=K_{-a}\) and \(D_a^{-1}=D_{-a}\). Consequently, the inverse step is \(K_{-\epsilon/2}\circ D_{-\epsilon}\circ K_{-\epsilon/2}\).

Let \(R(q,p)=(q,-p)\) be momentum reversal. Applying \(R\) before and after a kick changes its duration from \(a\) to \(-a\); the same is true for a drift. Therefore \(R\Phi_\epsilon R=\Phi_\epsilon^{-1}\). In practical terms, if a leapfrog step takes \((q,p)\) to \((q',p')\), applying the same step to \((q',-p')\) returns \((q,-p)\). The momentum flip is essential to this statement; the map is not generally its own inverse without it.

Volume preservation follows by examining the Jacobian matrices of the substeps. If \(U\) is twice differentiable, the kick has the block-triangular Jacobian $$ J_K= \begin{pmatrix} I&0\\ -a\nabla^2U(q)&I \end{pmatrix}, \qquad \det(J_K)=1. $$ The drift has Jacobian $$ J_D= \begin{pmatrix} I&aM^{-1}\\ 0&I \end{pmatrix}, \qquad \det(J_D)=1. $$ A composition has determinant equal to the product of the component determinants, so a leapfrog step has Jacobian determinant one. It therefore preserves local volume in \((q,p)\)-space. This does not mean that the map leaves every distance or shape unchanged; it means that infinitesimal volumes are not expanded or compressed.

Energy Error and the Metropolis Correction

Exact Hamiltonian dynamics conserve \(H(q,p)\), but a finite leapfrog step generally does not. After \(L\) steps, write the numerical endpoint as \((q_L,p_L)\) and define the energy error as $$ \Delta H=H(q_L,p_L)-H(q_0,p_0). $$ A positive value means the proposed endpoint has higher Hamiltonian energy; a negative value means it has lower energy. The sign alone does not determine whether the numerical integration is accurate: the size of the error matters.

HMC uses a Metropolis acceptance probability \(\min(1,\exp(-\Delta H))\) for this proposal, with the usual momentum reversal convention. This correction compensates for the fact that leapfrog does not conserve the exact Hamiltonian. Reversibility and volume preservation are what allow the proposal and its reverse to be compared without an extra volume factor. The correction does not make a poor numerical trajectory useful: large energy errors can make proposals likely to be rejected.

Diagnostics, Alternatives, and Limitations

A practical check is to examine energy errors across representative starting states and trajectories. Errors that are unexpectedly large, highly variable, or sensitive to small changes in step size can indicate that the discretization is not resolving the local curvature of the potential. A second check is a numerical reversibility test: take a step, reverse momentum, and step again. The recovered state should agree with the start up to numerical rounding. Failure beyond that scale can point to inconsistent update order, inaccurate gradient evaluations, or numerical instability. These checks do not establish that the posterior is well explored; they assess the integration map itself.

Leapfrog is not the only way to approximate differential equations, but its structural properties distinguish it from many generic integrators. A straightforward Euler update is simpler, yet generally lacks the same combination of reversibility and volume preservation. Higher-order Runge–Kutta methods can provide accurate solutions to ordinary differential equations, but ordinary implementations do not automatically preserve these HMC properties. Exact Hamiltonian solutions are preferable when available, but are uncommon for complex NLME posteriors. A random-walk Metropolis proposal avoids integrating Hamilton’s equations, but it does not use their gradient information to construct trajectories.

The guarantees described here have limits. They presume a fixed mass matrix and a well-defined gradient at the evaluated positions. Nondifferentiable points, discontinuities, invalid parameter regions, or unreliable gradient calculations can undermine the numerical map. Even with a smooth potential and correct implementation, leapfrog only approximates the trajectory, and finite precision can introduce small discrepancies. A valid Metropolis correction addresses the energy change of a correctly constructed proposal; it cannot repair an incorrectly implemented update or an invalid gradient.

A common mistake is to say, “Leapfrog conserves energy, so its endpoint must have the same Hamiltonian as its start.” That confuses an exact property of Hamiltonian flow with a property of its numerical approximation. For a finite leapfrog step, \(\Delta H\) is generally nonzero; the method’s useful structural guarantees are reversibility and volume preservation, while the Metropolis correction accounts for energy error. Another mistake is to omit the final half momentum update or evaluate its gradient at the old position. Either change breaks the leapfrog map defined here and can invalidate the reversibility argument.