Source-linked AI summary
Coupling and Convergence for Hamiltonian Monte Carlo
Nawaf Bou-Rabee, Andreas Eberle, Raphael Zimmer
TL;DR
HMC can suffer from the diffusive behavior of random-walk MCMC, motivating methods with faster state-space exploration. This paper introduces a velocity coupling and proves explicit contraction in a tailored L1-Wasserstein metric. The resulting bounds show that appropriately choosing the Hamiltonian duration can overcome diffusive behavior, while numerical and unadjusted variants retain scope limitations.
Problem
Random-walk MCMC can be slow because of diffusive behavior, while the mathematical properties of HMC are not fully understood.
Method
The paper couples velocities in two HMC copies and analyzes their transition kernel using an explicit Kantorovich L1-Wasserstein metric.
Results
The transition kernel is contractive with an explicit rate, and the rate can depend on K and R/T rather than explicitly on dimension.
Takeaways & Limitations
Choosing the duration T proportional to R can give a target approximation after total Hamiltonian-dynamics time of kinetic order O(R).
Takeaways & Limitations
For adjusted numerical HMC, the bounds require h^-1 of order d, and improving this dimension dependence remains open.
Abstract
from arXiv · showhide
Based on a new coupling approach, we prove that the transition step of the Hamiltonian Monte Carlo algorithm is contractive w.r.t. a carefully designed Kantorovich (L1 Wasserstein) distance. The lower bound for the contraction rate is explicit. Global convexity of the potential is not required, and thus multimodal target distributions are included. Explicit quantitative bounds for the number of steps required to approximate the stationary distribution up to a given error are a direct consequence of contractivity. These bounds show that HMC can overcome diffusive behaviour if the duration of the Hamiltonian dynamics is adjusted appropriately.
1. Introduction
HMC is introduced as a way to address the slow, diffusive behavior of random-walk MCMC, but its mathematical properties are not fully understood. The paper analyzes HMC through a new coupling that yields explicit contraction and convergence bounds without requiring global convexity.
- Random-walk MCMC can be slow because its meandering behavior produces diffusive sampling.
- HMC adds momentum to enable faster exploration, yet its mathematical properties remain incompletely understood.
- The paper couples velocities of two HMC copies so their positions approach each other with maximal probability in the free case.
- The resulting transition kernel is contractive in an explicit Kantorovich L1-Wasserstein metric equivalent to Euclidean distance.
- Under suitable regularity and outside-ball strong-convexity assumptions, the contraction rate depends on K and R/T rather than explicitly on dimension.
- Choosing T proportional to R can yield approximation after total Hamiltonian-dynamics time of kinetic order O(R), while numerical HMC extends the results under sufficiently small step sizes.
2.1. Hamiltonian Monte Carlo
HMC evolves positions using Hamiltonian dynamics after drawing Gaussian velocities, optionally followed by numerical integration and accept/reject correction. Exact and adjusted numerical HMC preserve the target measure under the stated assumptions, whereas unadjusted numerical HMC generally targets a different invariant measure.
- The exact flow solves dq_t/dt = p_t and dp_t/dt = −∇U(q_t) from initial state (x,v).
- Numerical HMC approximates the flow with a velocity Verlet integrator, with an accept/reject step defining adjusted versus unadjusted variants.
- The target measure is invariant for exact and adjusted numerical HMC, while unadjusted numerical HMC generally has a different invariant measure that approaches the target as h decreases.
2.2. Assumptions
The analysis assumes smoothness and bounded derivatives of the potential, together with either strong convexity outside a ball or a Lyapunov drift condition. These conditions support explicit bounds for exact and numerical HMC, with adjusted numerical HMC requiring additional finite-region and step-size restrictions.
- The potential U is assumed to be C4 with bounded second, third, and fourth derivatives.
- One main condition requires U to be strongly convex outside a Euclidean ball with radius R and constant K.
- The strong-convexity condition can be replaced by a Lyapunov drift condition involving a function Ψ and constants λ, α, and R2.
- For adjusted numerical HMC, global Lyapunov functions may fail, so contraction can require restricting the analysis to a finite ball and choosing h sufficiently small.
- Quadratic Lyapunov functions yield explicit drift bounds for exact and unadjusted HMC when the discretization conditions are satisfied.
- Exponential Lyapunov functions are available under a radial drift condition, with δ = κ/(4C + 8d + Q^2T^2) for exact or unadjusted HMC.
2.3. Coupling
The coupling uses shared randomness when positions are far apart and a velocity coupling with reflection otherwise, maximizing the chance of a contracting initial velocity difference.
- Far-apart states: For distant states, synchronous coupling reuses the same momentum and acceptance variables, while strong convexity supports contractivity.This construction applies when |x−y| ≥ 2R.
- Nearby states: For nearby states, the coupling uses a velocity difference ξ−η=−γz with maximal probability and reflection coupling otherwise.Here z=x−y and γ>0 is chosen appropriately.
- Transition construction: The coupled transitions use the same uniform random variable to decide acceptance in the adjusted case, while the unadjusted case accepts every proposal.The resulting second component is defined by either the proposed endpoint or the original position.
- Coupling rationale: A negatively proportional initial velocity difference makes the position difference contract over an initial time interval.This is the dynamical motivation for the coupling illustrated in Figure 1.
2.4. Numerical illustration of couplings
The numerical experiments visualize coupled trajectories on multimodal distributions and compare how quickly different coupling choices bring their components together.
- Experimental setup: The experiments use normal-mixture and Laplace-mixture target distributions, with the coupling applied globally and a small step size chosen so nearly all proposals are accepted.The normal mixture contains twenty two-dimensional Gaussian components; the Laplace mixture uses twenty two-dimensional regularized Laplace components.
- Figure 2: With T=1 and γ=1, Figure 2 depicts coupling components as colored dots over potential contours, with smaller markers representing later trajectory steps.Insets show the inter-component distance against the step index.
- Figure 3: Figure 3 estimates the average time to reach distance 10^-9 over 10^5 coupled samples and 100 duration values, comparing γ=T^-1 with synchronous coupling γ=0.The inverse-duration choice is motivated by Figure 1.
2.5. Contractivity
The paper establishes contractivity for exact and numerical HMC under explicit conditions, extending beyond globally convex potentials while identifying step-size and dimensionality boundaries.
- Exact HMC: The contraction condition requires LT^2 to remain below a specified threshold, and longer durations can fail because Hamiltonian flows may be periodic.A necessary duration bound can prevent such periodicities in the strongly convex setting.
- Exact HMC: Exact HMC admits global contraction under the stated assumptions, with a transition-rate bound for sufficiently small LT^2 and suitable coupling distance.The strongly convex result first gives contraction for sufficiently separated states; the general result extends this globally.
- Numerical HMC: Numerical HMC is contractive on finite balls for sufficiently small discretization steps, with larger balls requiring smaller steps.The contraction rate itself does not depend on the ball radius in the general numerical result.
- Limitations: The proofs require globally Lipschitz ∇U, while adjusted HMC bounds require h^-1 of order d and therefore introduce at least the same order of dimension dependence in computational complexity.The authors leave improvement of this adjusted-HMC dependence open.
- Beyond strong convexity: Under a global Lyapunov condition, exact and unadjusted numerical HMC retain global average contractivity with respect to a semimetric.Adjusted numerical HMC instead has contractivity on a given ball for sufficiently small step sizes.
- Rates and dimension: For fixed K and R/T, the contraction rate need not deteriorate as R increases when T is chosen proportional to R.The rate has no explicit dimension dependence, although model parameters can depend on dimension.
2.6. Quantitative bounds for distance to the invariant measure
The paper derives explicit Wasserstein-distance bounds for convergence to HMC’s invariant measure, using contractivity together with Lyapunov control. These results quantify the required number of steps and extend, under sufficiently small discretization, to numerical HMC.
- Exact HMC is globally contractive in a carefully designed Kantorovich (L1 Wasserstein) distance.
- For a target error ϵ, the number of exact-HMC steps can be chosen explicitly from the contraction bound.
- Choosing T proportional to R can make the required step count universal for a fixed initial error and tolerance, subject to LR^2 remaining bounded.
- The Wasserstein results also imply explicit bounds for biases, variances, and concentration of HMC ergodic averages.
- Numerical HMC combines local contractivity on a ball with Lyapunov bounds controlling exits from that ball.
- Adjusted numerical HMC achieves a similar step count to exact HMC when h is sufficiently small, whereas unadjusted HMC incurs additional invariant-measure bias.
3. A priori estimates
The a priori analysis establishes bounds for Hamiltonian and numerical flows, the velocity coupling, and acceptance-rejection discrepancies. These ingredients support average contractivity for exact and numerical HMC under explicit conditions.
- The analysis develops bounds for Hamiltonian flow, coupling, and acceptance-rejection probabilities used in the main contraction proof.
- For velocity Verlet, short-time flow differences are controlled by the corresponding constant-velocity motion, yielding contraction when u−v=−γ(x−y).
- Some numerical-flow bounds require restrictions on |v| or on the state and discretization scales, unlike the corresponding exact-flow result.
- Exact Hamiltonian dynamics conserve the Hamiltonian, while velocity Verlet incurs an energy error controlled by h^2 and the initial state magnitudes.
- The coupling makes ξ−η=−γz with maximal probability and uses reflection coupling otherwise, enabling control of the exceptional event.
- Acceptance-rejection discrepancies are bounded by terms proportional to T(1+T)h^2 and the position separation.
4. Proofs of a priori bounds
The proofs obtain the a priori estimates by comparing numerical trajectories with simpler motions, deriving differential inequalities, and controlling Hamiltonian and acceptance-rejection errors.
- Trajectory and momentum bounds are derived from the velocity Verlet recursion and Lipschitz control of ∇U.
- Under the stated short-time and separation conditions, the position difference remains controlled and the trajectories retain sufficient separation.
- For separated trajectories, the proof tracks a(t)=|z_t|^2 and b(t)=2z_t·w_t through a coupled differential system.
- The numerical Hamiltonian error is bounded by combining estimates for positions, velocities, and derivative flows.
- Acceptance-event estimates follow by expressing acceptance through Hamiltonian errors and applying the resulting trajectory bounds.
- The coupling estimates use reflection on the complementary event and moment bounds for the Gaussian velocity component.
5. Proofs of main results
The proofs establish contraction for exact and numerical HMC by combining coupled Hamiltonian proposals with bounds on rejection and discretization errors. Exact HMC inherits the same contraction rate, while numerical HMC requires sufficiently small step sizes.
- Exact and numerical HMC: The coupling event A(x) ∩ A(y) occurs with probability at least 4/5 under the selected step-size conditions.This probability bound is used in the contraction estimate for numerical HMC.
- Exact and numerical HMC: For adjusted numerical HMC, acceptance probabilities and rejection-event mismatches are bounded so perturbation terms remain smaller than the contraction term.The proof controls these terms by choosing h below several parameter-dependent thresholds.
- Global contraction: The proof combines separate distance regimes to obtain global contraction for adjusted numerical HMC when h is sufficiently small.The global result follows by combining the bounds for |x − y| ≥ 2R and |x − y| < 2R.
- Global contraction: For exact HMC, rejection events are absent, so the contraction bound holds globally with the same rate c.The exact-HMC proof is obtained from the numerical-HMC argument after removing rejection-related terms.
6. Proofs of results in Section 2.6
This section turns local contractivity into convergence bounds by using couplings, stopped supermartingales, and Lyapunov control of excursions outside the contractive region. The resulting estimates provide explicit choices of iteration count and Lyapunov threshold for a target error.
- Coupling and supermartingales: A locally contractive coupling yields a stopped process M_n = e^{c(n∧T)}ρ(X_{n∧T}, Y_{n∧T}) that is a non-negative supermartingale.The stopping time records the first exit of either chain from the contractive set.
- Exact HMC convergence: The exact-HMC error bound follows by applying the local coupling result to arbitrary initial laws and then taking the infimum over couplings.The metric ρ is equivalent to Euclidean distance through e^{-aR1}|x − y| ≤ ρ(x, y) ≤ |x − y|.
- Numerical HMC convergence: Choosing n so the contraction contribution is at most ϵ/2 and C sufficiently large makes the total distance satisfy ∆(n) ≤ ϵ.The argument explicitly separates the within-region contraction term from the exit-related term.
Appendix A: Explicit Lyapunov functions
The appendix verifies Lyapunov conditions for exact, unadjusted, and adjusted numerical HMC using explicit energy-growth estimates. These estimates establish stability globally in some cases and on bounded regions for adjusted numerical HMC.
- Exact and unadjusted HMC: For exact and unadjusted HMC, the Lyapunov drift condition holds globally when n is sufficiently large and L(T^2 + hT) is sufficiently small.The proof uses Ψ(x) = |x|^2 for these algorithms.
- Adjusted numerical HMC: For adjusted numerical HMC, the corresponding Lyapunov assertion is restricted to |x| ≤ R2 because discretization and rejection control may require smaller step sizes on larger balls.The restriction depends on a finite radius R2 and additional step-size conditions.
- Explicit Lyapunov functions: The appendix also constructs exponential-growth Lyapunov functions Ψ(x) = g(|x|^2) for exact HMC under the stated assumptions.The function is chosen to agree with exp(δ|x|) outside a bounded region.
- Explicit Lyapunov functions: For numerical HMC, the Lyapunov verification combines bounds on the numerical Hamiltonian energy with control of rejection probabilities.The resulting conditions support the Lyapunov assumptions used in the convergence proof.