Hamiltonian Dynamics: Phase Space, Liouville, and Why HMC Works
Published:
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:
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
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
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:
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.
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.
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
- Neal, R. M. MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo, Chapman & Hall/CRC, 2011.
- Duane, S., Kennedy, A. D., Pendleton, B. J., & Roweth, D. Hybrid Monte Carlo. Physics Letters B, 195(2):216–222, 1987.
- Betancourt, M. A Conceptual Introduction to Hamiltonian Monte Carlo. 2017.
- Hairer, E., Lubich, C., & Wanner, G. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, 2nd ed. Springer, 2006.