Alessio Borgi

Alessio Borgi

PhD in Graph Neural Networks and Generative AI 👨🏻‍💻🤖

Hamiltonian Dynamics: Phase Space, Liouville, and Why HMC Works

6 minute read

Published:

TL;DR: A Legendre transform swaps velocity for momentum and turns one second-order equation into two first-order ones on phase space. The resulting flow conserves energy and, by Liouville's theorem, conserves volume. HMC borrows exactly those two properties: volume preservation makes the proposal Jacobian equal to 1 and reversibility makes the proposal symmetric, so the Metropolis ratio reduces to \(e^{-\Delta H}\). A symplectic integrator such as leapfrog keeps both properties exactly; an Euler step keeps neither and the energy diverges.

From velocities to momenta

Start from the Lagrangian \(L(q,\dot q)\). Define the conjugate momentum and perform a Legendre transform in the velocity variable:

\[ p \;=\; \frac{\partial L}{\partial \dot q}, \qquad\quad H(q,p) \;=\; p\,\dot q \;-\; L(q,\dot q), \]

where \(\dot q\) on the right is eliminated in favour of \(p\). For \(L = \tfrac12 m\dot q^2 - V(q)\) we get \(p = m\dot q\), hence \(\dot q = p/m\) and

\[H = \frac{p^2}{m} - \left(\frac{p^2}{2m} - V(q)\right) = \frac{p^2}{2m} + V(q),\]

which is the total energy. The transform trades a function of position and velocity for a function of position and momentum; nothing is lost, because the map \(\dot q \mapsto p\) is invertible whenever \(L\) is convex in \(\dot q\).

Hamilton’s equations

Taking the differential of \(H = p\dot q - L\) and using Euler–Lagrange gives the pair

\[ \dot q \;=\; \frac{\partial H}{\partial p}, \qquad\qquad \dot p \;=\; -\,\frac{\partial H}{\partial q}. \]

The minus sign sits on the momentum equation only. For the example above, \(\dot q = p/m\) and \(\dot p = -V'(q)\), which is \(F = ma\) split into two first-order pieces.

Two consequences follow immediately. First, energy conservation:

\[\frac{dH}{dt} = \frac{\partial H}{\partial q}\dot q + \frac{\partial H}{\partial p}\dot p = \frac{\partial H}{\partial q}\frac{\partial H}{\partial p} - \frac{\partial H}{\partial p}\frac{\partial H}{\partial q} = 0\]

for a time-independent \(H\): the two terms cancel exactly. Second, Liouville’s theorem. The flow on phase space has velocity field \((\dot q,\dot p)\), whose divergence is

\[\frac{\partial \dot q}{\partial q} + \frac{\partial \dot p}{\partial p} = \frac{\partial^2 H}{\partial q\,\partial p} - \frac{\partial^2 H}{\partial p\,\partial q} = 0 .\]

A divergence-free flow transports volume without changing it. A blob of initial conditions may be stretched and folded beyond recognition, but its phase-space volume is invariant.

Why HMC needs exactly these two facts

To sample \(\pi(q) \propto e^{-U(q)}\) with \(U(q) = -\log \pi(q)\), HMC invents a momentum \(p \sim \mathcal N(0, M)\) and works with the joint density

\[ \pi(q,p) \;\propto\; e^{-H(q,p)}, \qquad H(q,p) = U(q) + \tfrac12 p^\top M^{-1} p . \]

The marginal in \(q\) is the target, so any correct sampler for the joint gives a correct sampler for \(\pi\). A proposal is generated by drawing a fresh \(p\) and simulating Hamiltonian dynamics for \(L\) steps, then negating the momentum. Metropolis–Hastings would in general require the proposal density ratio and a Jacobian determinant. Here both disappear:

  • Volume preservation \(\Rightarrow\) the Jacobian determinant of the proposal map is exactly 1, so no correction term.
  • Reversibility (the flow with the momentum flipped retraces the path) \(\Rightarrow\) the proposal is symmetric, so the density ratio cancels.

What survives is \(\min\!\left(1,\ e^{-H(q^{\ast},p^{\ast})+H(q,p)}\right)\). If the dynamics were simulated exactly, \(H\) would be conserved and every proposal accepted, however far it travelled. Discretisation error is the only source of rejection, which is why the integrator matters so much (Neal, 2011).

Leapfrog versus Euler, with numbers

Leapfrog takes a half kick, a full drift, then a half kick:

\[ p_{t+\frac{\epsilon}{2}} = p_t - \tfrac{\epsilon}{2}\nabla U(q_t), \quad q_{t+\epsilon} = q_t + \epsilon M^{-1} p_{t+\frac{\epsilon}{2}}, \quad p_{t+\epsilon} = p_{t+\frac{\epsilon}{2}} - \tfrac{\epsilon}{2}\nabla U(q_{t+\epsilon}). \]

Test both on the harmonic oscillator \(H = \tfrac12(q^2+p^2)\), so \(\nabla U(q)=q\) and \(M=1\). Explicit Euler gives \(q' = q+\epsilon p\), \(p' = p - \epsilon q\), whose Jacobian determinant is \(1+\epsilon^2 > 1\). Worse, the energy obeys \(q'^2+p'^2 = (1+\epsilon^2)(q^2+p^2)\) exactly, it multiplies by \(1+\epsilon^2\) every single step. With \(\epsilon = 0.1\) over 1000 steps the energy grows by a factor of \(1.01^{1000} \approx 2.1\times10^4\). Every proposal is rejected.

Composing the three leapfrog lines instead gives the linear map with matrix

\[A = \begin{pmatrix} 1-\tfrac{\epsilon^2}{2} & \epsilon \\[2pt] -\epsilon\left(1-\tfrac{\epsilon^2}{4}\right) & 1-\tfrac{\epsilon^2}{2}\end{pmatrix}, \qquad \det A = \left(1-\tfrac{\epsilon^2}{2}\right)^2 + \epsilon^2\left(1-\tfrac{\epsilon^2}{4}\right) = 1 .\]

The \(\epsilon^4\) terms cancel and the determinant is exactly 1, for any step size. Numerically, running \(\epsilon = 0.3\) from \((q,p)=(1,0)\) for 22 steps, the leapfrog radius stays between 0.989 and 1.000, while Euler’s grows to 2.58.

Phase-space trajectories of leapfrog and explicit Euler on a harmonic oscillatorTwo trajectories in the q-p plane starting at q equals 1, p equals 0, integrated for 22 steps with step size 0.3. The leapfrog trajectory stays on a closed loop with radius between 0.989 and 1.000. The explicit Euler trajectory spirals outwards, reaching radius 2.58 after 22 steps, so its energy grows without bound. start (1, 0) q p leapfrog, radius 0.989 to 1.000 explicit Euler, radius grows to 2.58 step size ε = 0.3, 22 steps, H = (q²+p²)/2 Euler energy multiplies by (1+ε²) each step; leapfrog Jacobian determinant is exactly 1.
Notice that leapfrog does not conserve energy exactly, its radius wobbles by about 1%, but the error oscillates instead of accumulating. Euler's error compounds multiplicatively, which is what destroys the acceptance rate.
Key Insight, bounded, not small: the virtue of a symplectic integrator is not higher accuracy per step. It is that leapfrog exactly conserves a nearby "shadow" Hamiltonian, so its energy error oscillates within a bounded band over arbitrarily long trajectories rather than drifting. HMC needs long trajectories to decorrelate; an integrator whose energy error grows with path length is useless no matter how small its local truncation error is.

Practical consequences

Tuning HMC is tuning \(\epsilon\) and the number of leapfrog steps \(L\). Too large an \(\epsilon\) and the shadow Hamiltonian argument fails, the integrator becomes unstable and acceptance collapses to near zero, sharply rather than gradually. Too small and you pay gradient evaluations for nothing. The mass matrix \(M\) acts as a preconditioner: setting \(M^{-1}\) to the target’s covariance equalises the timescales, exactly as in the Boltzmann picture where the momenta thermalise independently. Momentum resampling between trajectories is what makes the chain ergodic, pure Hamiltonian flow would stay on one energy level set forever.

Interview trap: "leapfrog is used because it is more accurate" is the wrong answer. It is second-order accurate, which is unremarkable; Runge–Kutta 4 is far more accurate per step and still ruins HMC, because it is not volume-preserving and its energy drifts. The properties that matter are exact volume preservation (Jacobian determinant 1, no correction in the acceptance ratio) and reversibility (symmetric proposal). Also note that HMC does not require exact energy conservation, the Metropolis step corrects whatever error remains.

Recap

  • Legendre transform: \(p = \partial L/\partial\dot q\), \(H = p\dot q - L\); for \(L = \tfrac12 m\dot q^2 - V\) this gives \(H = p^2/2m + V\), the energy.
  • Hamilton's equations are \(\dot q = \partial H/\partial p\) and \(\dot p = -\partial H/\partial q\), the minus is on the momentum equation.
  • The phase-space flow is divergence-free, so volume is conserved (Liouville), and \(dH/dt = 0\) for time-independent \(H\).
  • HMC's acceptance ratio reduces to \(e^{-\Delta H}\) precisely because volume preservation kills the Jacobian and reversibility kills the proposal ratio.
  • Leapfrog has determinant exactly 1 and bounded energy error; explicit Euler multiplies energy by \(1+\epsilon^2\) per step and rejects everything.

References

  1. Neal, R. M. MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo, Chapman & Hall/CRC, 2011.
  2. Duane, S., Kennedy, A. D., Pendleton, B. J., & Roweth, D. Hybrid Monte Carlo. Physics Letters B, 195(2):216–222, 1987.
  3. Betancourt, M. A Conceptual Introduction to Hamiltonian Monte Carlo. 2017.
  4. Hairer, E., Lubich, C., & Wanner, G. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, 2nd ed. Springer, 2006.

You May Also Enjoy

The Standard Library Modules Worth Knowing Before You pip install

8 minute read

Published:

Python ships with a counter, a queue, a memoiser, a CLI parser, a test framework’s smaller cousin, and a JSON-safe date library already installed, most of what people reach for a package to solve is one import away.

Files and Context Managers: Why with Is Not Optional

9 minute read

Published:

A file object that falls out of scope without with does get closed eventually, but ‘eventually’ means whenever the garbage collector gets around to it, which is not a promise any program handling more than a handful of files can live with.

Errors and Exceptions: EAFP, the Hierarchy, and Reading a Traceback

9 minute read

Published:

Exceptions in Python are not exceptional. They are a normal control-flow mechanism that the language leans on so heavily that the idiomatic style is to try the operation and handle the failure, rather than check first, which is faster, shorter, and free of race conditions.