PaperlayerView exact revision on arXiv

Source-linked AI summary

A Conceptual Introduction to Hamiltonian Monte Carlo

Michael Betancourt

arXiv:1701.02434v2stat.ME

TL;DR

This review addresses why Hamiltonian Monte Carlo succeeds on difficult, high-dimensional problems and how it should be applied in practice. It develops an intuitive account of the geometry underlying HMC, its construction, implementation, and diagnostics, concluding that coherent trajectories efficiently explore the target while practical performance remains constrained by geometric and numerical issues.

  • Problem

    Efficient statistical computation is difficult because finite computational resources and high-dimensional geometry limit how effectively target-density evaluations can approximate expectations.

  • Method

    The paper develops an intuitive, geometry-based account of Hamiltonian Monte Carlo, including its transition construction, geometric adaptations, and practical diagnostics.

  • Results

    Hamiltonian trajectories rapidly explore the typical set, while momentum sampling and projection move chains from the tails toward it for sufficiently well-behaved targets.

  • Takeaways & Limitations

    Understanding the target geometry motivates HMC and its efficient practical implementations, while also clarifying when the method can fail or require additional corrections.

  • Takeaways & Limitations

    Generic numerical integrators accumulate coherent drift on long trajectories, with drift worsening in higher-dimensional phase space and limiting them to short, inefficient trajectories.

Abstract

from arXiv · show

Hamiltonian Monte Carlo has proven a remarkable empirical success, but only recently have we begun to develop a rigorous understanding of why it performs so well on difficult problems and how it is best applied in practice. Unfortunately, that understanding is confined within the mathematics of differential geometry which has limited its dissemination, especially to the applied communities for which it is particularly important. In this review I provide a comprehensive conceptual account of these theoretical foundations, focusing on developing a principled intuition behind the method and its optimal implementations rather of any exhaustive rigor. Whether a practitioner or a statistician, the dedicated reader will acquire a solid grasp of how Hamiltonian Monte Carlo works, when it succeeds, and, perhaps most importantly, when it fails.

1. COMPUTING EXPECTATIONS BY EXPLORING PROBABILITY DISTRIBUTIONS

Computing expectations requires numerical methods because nontrivial integrals cannot generally be evaluated analytically, while high-dimensional geometry makes efficient evaluation difficult. The density–volume balance concentrates meaningful contributions in the typical set, so computation should focus there rather than near the mode or throughout the full space.

  • 1.1 Computing Expectations in Practice: Expectations over target distributions are represented as integrals, but nontrivial integrals require finite-resource numerical approximations.The target distribution is defined on a smooth parameter space, with expectations expressed through a probability density and integrals.
  • 1.1 Computing Expectations in Practice: Efficient methods must use each evaluation of the target density and relevant functions because high-dimensional geometry frustrates computation.Wasting evaluations in regions with negligible contributions limits the utility of computational methods.
  • 1.2 Parsimonious Expectation Computation: Focusing only on regions where density or target functions are largest is inadequate because expectation contributions depend on both density and volume.The mode has high density but little volume, while distant regions have large volume but vanishing density.
  • 1.3 The Geometry of High-Dimensional Spaces: In high dimensions, volume away from the mode grows rapidly and can compensate for smaller density, shifting substantial expectation contributions away from the mode.The volume disparity grows exponentially with dimension, making neighborhoods around the mode increasingly negligible.
  • 1.4 The Geometry of High-Dimensional Probability Distributions: The typical set is the intermediate region where density and volume are jointly large enough to contribute significantly, and it is invariant to reparameterization.Regions immediately around the mode and far in the tails contribute little, while the typical set narrows as dimension increases.
  • 1.4 The Geometry of High-Dimensional Probability Distributions: Naive quadrature scales poorly because a uniform D-dimensional grid requires N^D integrand evaluations and is unlikely to intersect the narrow typical set.Framing algorithms by how they quantify the typical set helps assess their practical performance in high dimensions.

2. MARKOV CHAIN MONTE CARLO

Markov chain Monte Carlo explores the typical set through successive Markov transitions, enabling expectation estimates when the chain adequately explores the target distribution. In finite computation, slow or incomplete exploration can bias estimators, reduce effective precision, and undermine ideal convergence guarantees.

  • Estimating expectations: Markov chain Monte Carlo stochastically explores the typical set to generate samples for estimating expectations.A properly designed chain can eventually explore the typical set, but finite-time exploration is the practical challenge.
  • Markov chains: A Markov chain is generated by repeatedly applying a transition density T(q′ | q) that determines the next state from the current state.When the transition preserves the target distribution, ensembles remain distributed according to that target.
  • Estimating expectations: Target-preserving transitions drive chains toward and across the typical set regardless of their initial position, making their history a quantification of that set.Averaging the target function over the chain history estimates expectations across the typical set and, hence, the full parameter space.
  • Ideal behavior: With sufficient exploration, Markov chain Monte Carlo estimators converge toward true expectations, but this asymptotic result is difficult to guarantee with finite computational resources.Ideal behavior includes initial bias during convergence, rapid improvement after entering the typical set, and eventual Central Limit Theorem behavior.
  • Ideal behavior: Effective sample size measures the number of independent target samples needed to match a chain’s estimator precision and depends on autocorrelation.It can also be interpreted as the number of sojourns the chain makes across the typical set.
  • Pathological behavior: High curvature and pathological regions can cause incomplete exploration, estimator bias, oscillations around true expectations, and failure of Central Limit Theorems.Random Walk Metropolis additionally explores the typical set extremely slowly except in the lowest-dimensional spaces, producing high autocorrelation and imprecise or biased estimates.

3. THE FOUNDATIONS OF HAMILTONIAN MONTE CARLO

Hamiltonian Monte Carlo addresses the failure of uninformed or diffusive transitions in high-dimensional typical sets by generating coherent exploration from geometric information. Its construction uses gradients, auxiliary momentum, and Hamiltonian dynamics to move rapidly through the target distribution while preserving the relevant typical-set structure.

  • Random Walk Metropolis fails in high dimensions because most proposal directions leave the typical set, while smaller steps remain accepted but diffuse locally.Large proposal variances are rejected; small variances concentrate transitions near the initial point.
  • Effective transitions must follow contours of high probability mass to move coherently into new, unexplored neighborhoods of the typical set.Aligned vector fields provide directions that act as guide posts through the typical set.
  • The target-density gradient contains geometric information but points toward parameterization-sensitive regions such as the mode rather than along the typical set.Additional geometric constraints are needed to remove parameterization dependence and twist directions into alignment with the typical set.
  • Hamiltonian Monte Carlo introduces auxiliary momentum with a probabilistic structure that generates conservative dynamics and coherent exploration.The paper presents HMC as the unique procedure for automatically generating this exploration for sufficiently well-behaved target distributions.
  • Hamilton’s equations channel the gradient through momentum, producing trajectories that rapidly traverse phase space while remaining constrained to the canonical typical set.Projecting these trajectories back to parameter space yields exploration of the target distribution.
  • The resulting Hamiltonian Markov transition combines random trajectory generation, momentum sampling, and projection to explore the target efficiently from tails and across the typical set.The trajectories push transitions away from previously explored neighborhoods while staying within the joint typical set.

4. EFFICIENT HAMILTONIAN MONTE CARLO

Efficient Hamiltonian Monte Carlo separates deterministic exploration within energy level sets from stochastic movement between them, making tuning central to performance. The review relates kinetic energy, integration time, and geometry to practical efficiency, while emphasizing that poorly chosen configurations can fail dramatically.

  • Microcanonical geometry: Hamiltonian Monte Carlo transitions decompose into deterministic trajectories exploring individual energy level sets and stochastic jumps between those level sets.The random lift and projection correspond to momentum resampling, while Hamiltonian trajectories remain within energy level sets.
  • Tuning the transition: Kinetic energy determines the geometry of energy level sets, while integration time controls how completely trajectories explore them and the computational cost of each transition.Short trajectories underuse coherent exploration; longer trajectories cost proportionally more integration time.
  • Tuning the transition: The efficiency of stochastic exploration depends on how quickly momentum resampling diffuses across energies relative to the marginal energy distribution.A narrow energy transition distribution relative to the marginal energy distribution slows diffusion across typical energies.
  • Riemannian geometry: A well-chosen metric can correct spatially varying correlations, producing more uniform level sets and more efficient exploration than global rotation and rescaling.Riemannian metrics vary across parameter space and can resemble the target Hessian; technical care is needed to keep the inverse metric well-behaved.
  • Integration-time selection: The No-U-Turn criterion has shown empirical success across diverse target distributions, while exhaustive termination criteria may be more robust but require careful tuning.The optimal integration time can vary with the energy level, especially for heavy-tailed targets where it grows toward the tails.

5. IMPLEMENTING HAMILTONIAN MONTE CARLO IN PRACTICE

Practical Hamiltonian Monte Carlo requires numerical integration because Hamilton’s equations are rarely solvable exactly. Generic integrators drift increasingly with trajectory length and dimension, motivating symplectic methods that preserve phase-space volume and remain near the exact energy level set.

  • Numerical integration: Most Hamiltonian trajectories must be generated numerically because Hamilton’s equations cannot be solved exactly except in trivial cases.Numerical accuracy directly affects the effectiveness of the resulting Hamiltonian Monte Carlo implementation.
  • Numerical integration: Generic numerical integrators accumulate error coherently, causing trajectories to drift away from the true trajectory and typical set, especially in high-dimensional phase space.This drift limits generic solvers to short trajectories that explore energy level sets inefficiently.
  • Symplectic integration: Symplectic integrators preserve phase-space volume exactly, preventing coherent drift and keeping numerical trajectories oscillating near the exact energy level set over long integrations.Their geometric preservation makes them robust to the drift that limits ordinary differential-equation solvers.

Appendix A.

Symplectic integrators make Hamiltonian Monte Carlo practical by preserving phase-space volume and keeping numerical trajectories near energy level sets, but implementation requires exact bias correction and reversible proposals. Their performance can fail in high-curvature regions, and some theoretical results do not yet extend formally to more sophisticated implementations.

  • Symplectic integration: Symplectic integrators preserve phase-space volume, limiting energy error and keeping numerical trajectories near the exact energy level set over long integrations.This near-conservation explains their practical accuracy compared with generic numerical integration.
  • Symplectic integration: With position-independent momentum distributions, the leapfrog integrator provides a simple implementation using alternating discrete momentum and position updates.The step size ϵ determines the time discretization, while the update interleaving preserves volume exactly.
  • Failure modes: Finite-step integration can fail in high-curvature neighborhoods, where numerical trajectories diverge toward infinite energy and rapidly approach phase-space boundaries.Unlike gradual drift from generic integrators, this failure is distinctive and straightforward to identify in practice.
  • Metropolis-Hastings correction: Because integrator error biases Hamiltonian transitions, Metropolis-Hastings correction is required; evaluating the Hamiltonian makes the acceptance probability analytically available.The resulting correction yields exact samples from the canonical distribution, which projects to the target distribution.
  • Reversible proposals: Forward-only deterministic trajectories produce zero Metropolis-Hastings acceptance because their reverse proposal probabilities vanish, whereas flipping momentum makes the proposal reversible and correctable.The momentum negation step gives non-zero forward and backward proposal-density ratios.
  • Trajectory proposals: Uniformly proposing points along a numerical trajectory can be suboptimal because high-error states may be selected and rejected despite lower-error states elsewhere on the trajectory.Averaging proposals over trajectories requires integrating both forwards and backwards, with schemes also extending to dynamic trajectory lengths.
  • Integrator configuration: Integrator configuration depends on step size ϵ and order K, but formal analysis for more sophisticated implementations remains incomplete and higher-order methods introduce serious complications.Higher-order integrators may be more effective in higher-dimensional problems, although the supplied passage does not specify the complication.

6. THE ROBUSTNESS OF HAMILTONIAN MONTE CARLO

Hamiltonian Monte Carlo’s robustness can be analyzed through geometric ergodicity and the pathological target-distribution behaviors that obstruct it. Energy-distribution mismatch and divergent trajectories provide practical diagnostics, but these checks are necessary rather than sufficient.

  • Hamiltonian Monte Carlo is geometrically ergodic for a larger class of target distributions than non-gradient methods such as Random Walk Metropolis.
  • Pathologies: Heavy tails can require prolonged excursions through the tails, while poorly chosen kinetic energies can make energy-level exploration too slow.
  • Diagnostics: Visualizing the marginal energy and energy-transition densities identifies inefficient momentum resampling between energy levels.
  • Diagnostics: Values of the energy Bayesian fraction of missing information below 0.3 have proven problematic, although the exact threshold remains theoretically unformalized.
  • Diagnostics: High-curvature neighborhoods can destabilize symplectic integrators, making divergent numerical trajectories immediate indicators of pathological regions.
  • Limitations: These diagnostics are necessary but not sufficient for estimator validity, so subtle pathologies may remain; combining them with split R-hat and multiple chains improves robustness.

7. CONCLUSION

The paper concludes that Hamiltonian Monte Carlo exploits typical-set geometry to achieve coherent exploration, efficient computation, and stronger estimator-validity guarantees. It also identifies geometric perspectives as a foundation for improving implementations and extending these principles to more complex spaces.

  • Hamiltonian Monte Carlo exploits typical-set geometry to generate coherent exploration of smooth target distributions.
  • This geometric exploration provides better computational efficiency and stronger guarantees on the validity of resulting estimators.
  • Geometric analysis motivates principled strategies for automatically constructing efficient implementations, allowing users to focus on building better models.
  • Future work includes improving geometric-ergodicity analysis, refining implementation techniques, and studying non-Gaussian and non-Euclidean kinetic energies.
  • Geometric methods may extend coherent exploration to more challenging computational problems and intricately structured spaces.

APPENDIX A: TECHNICAL DETAILS OF PRACTICAL IMPLEMENTATIONS

Hamiltonian Monte Carlo decomposes exploration into movement along energy level sets and momentum resampling between energy levels. Practical implementations approximate ideal Hamiltonian trajectories with symplectic integrators and probabilistic corrections.

  • Hamiltonian systems decompose into a microcanonical distribution over energy level sets and a marginal energy distribution between levels.
  • Hamiltonian Monte Carlo alternates trajectory-based exploration within an energy level set with momentum resampling across the marginal energy distribution.
  • Hamiltonian trajectories preserve energy and, when integrated longer, typically explore their entire level set.
  • Because analytic trajectories are generally unavailable, practical implementations use symplectic integrators whose numerical errors can be corrected probabilistically.
  • The appendix reviews fixed-length and dynamically terminated numerical trajectories and relates these implementations to the No-U-Turn sampler and Stan.

A.1 Notation

The appendix represents Hamiltonian Monte Carlo states in phase space and defines numerical trajectories through repeated forward or backward symplectic-integration steps. It then introduces notation for trajectories containing specified points.

  • A phase-space point is represented as z = (q, p), and a symplectic integrator generates a numerical trajectory from it using time increment ϵ.
  • Forward and backward integration from an initial point z0 produces points zn and z−n, respectively.
  • For length L = 3, example trajectories are contiguous sets of three points that may be centered at different positions along the integrated path.
  • TLz denotes the set of all symplectic-integrator trajectories of length L containing point z.
  • The notation is extended to trajectories of length L containing both points z and z′.

A.2 Static Implementations

Static implementations sample a numerical trajectory containing the current state, then sample a state from that trajectory. Correcting nonuniform trajectory errors requires weighted state sampling.

  • Symplectic integration produces a numerical trajectory containing the initial point, from which a subsequent state is sampled.
  • Uniform sampling along a numerical trajectory is generally incorrect because symplectic integration does not exactly preserve energy.
  • The transition is represented sequentially as sampling a trajectory T(t | z) followed by sampling a state T(z′ | t).
  • For static trajectories, T(t | z) is restricted to trajectories of fixed length L that contain the initial state z.

A.2.1 Sampling a Trajectory

Valid trajectory sampling requires reversibility: every state in a trajectory must induce the same probability of selecting that trajectory. Uniform selection over containing trajectories provides a practical construction.

  • Reversibility requires the probability of selecting a trajectory to be identical from every state on that trajectory.
  • This reversibility condition is necessary for a valid correction of the trajectory-based transition.
  • In practice, the condition is most easily satisfied by sampling uniformly among trajectories that contain the initial point.
  • After selecting a numerical trajectory, a transition proportional to canonical densities can generate an exact canonical-distribution sample.

A.2.2 Sampling a State

State sampling must compensate for integration-induced energy errors, using canonical-density-weighted transitions rather than uniform selection. Direct multinomial sampling avoids auxiliary variables and improves performance.

  • Integration errors create nonuniform state probabilities that disfavor numerical deviations toward higher energies.
  • A slice sampler can obtain the desired transition by sampling a threshold u and selecting trajectory states whose canonical density exceeds u.
  • Direct multinomial sampling assigns trajectory states probabilities proportional to their canonical densities.
  • Rao-Blackwellization implies that omitting auxiliary random variables yields better performance for this transition.
  • Volume preservation by symplectic integrators supports exact canonical sampling when the transition preserves canonical density.

A.2.3 Invariance of the Canonical Distribution

The paper develops progressive and dynamic trajectory construction while preserving canonical invariance. Dynamic termination adapts integration length, but reversibility, geometric scope, and tuning requirements constrain the available criteria and implementations.

  • A.2.3 Invariance of the Canonical Distribution: An explicit calculation shows that the stated transition properties are sufficient to guarantee invariance of the canonical distribution.
  • A.2.3 Invariance of the Canonical Distribution: Progressive sampling interleaves trajectory construction with state sampling, reducing memory from storing the entire trajectory to retaining only a few states.
  • A.2.3 Invariance of the Canonical Distribution: Trajectories can be expanded additively or by repeatedly doubling their length, with binary-tree structure supporting recursive construction and concurrent state updates.
  • A.2.3 Invariance of the Canonical Distribution: Bernoulli mixing combines samples from old and new trajectory components while preserving the required joint transition distribution.
  • A.2.3 Invariance of the Canonical Distribution: Biasing progressive updates away from the initial state can introduce anticorrelation while preserving the canonical distribution after averaging over initial states.
  • A.4 Dynamic Implementations: Static integration times are fundamentally limited because different energy level sets require different times for effective exploration.
  • A.4 Dynamic Implementations: Dynamic schemes expand trajectories until a termination criterion is met, but naive implementations can violate reversibility and canonical invariance.
  • A.4.2 Dynamic Termination Criteria: The original No-U-Turn criterion is limited to Euclidean manifolds, whereas its generalized form is well-defined on Riemannian manifolds and can be approximated discretely.
Loading 1701.02434v2…