A Conceptual Introduction to Hamiltonian Monte Carlo

1. COMPUTING EXPECTATIONS BY EXPLORING PROBABILITY DISTRIBUTIONS

• Expectation computation begins with a target distribution π on a D-dimensional sample space and its function expectations Eπ[f], represented as integrals after selecting a smooth parameterization.Different parameterizations produce different density representations but identical expectation values; the practical challenge is approximating these integrals under limited computational power.Source: §1 · pp. 3–4, §1.1 · p. 4

1.1 Computing Expectations in Practice

• Efficient expectation estimation requires numerical approximations that use every target-density and function evaluation effectively, because finite computation and high-dimensional geometry limit accuracy.Expectations become integrals after choosing a parameterization, but their values remain invariant across parameterizations; algorithms should therefore avoid depending on irrelevant density-shape details.Source: §1.1 · pp. 4–5

1.2 Parsimonious Expectation Computation

• Focusing only on the density mode is misleading because expectations accumulate over volume, and methods should generally support many sufficiently uniform target functions rather than one function-specific integrand.Exploiting a single function’s structure can help when only one expectation matters, but the general framework prioritizes regions determined by density and volume together.Source: §1.2 · p. 5

1.3 The Geometry of High-Dimensional Spaces

• In high-dimensional spaces, volume grows exponentially away from the mode and can compensate for lower density, so significant regions must be identified by balancing density against volume.Rectangular and spherical constructions both show that neighborhoods around the mode become negligible relative to surrounding volume as dimension increases.Source: §1.3 · pp. 6–7

1.4 The Geometry of High-Dimensional Probability Distributions

• The only regions making substantial contributions to high-dimensional expectations form the typical set, where sufficiently large density and volume balance between the mode and the tails.The mode has high density but little volume, while distant regions have large volume but vanishing density; the typical set is invariant under reparameterization.Source: §1.3 · p. 6, §1.4 · pp. 7–8
• As dimension increases, the typical set narrows through concentration of measure, so evaluating outside it has negligible effect and wastes computational resources.Averaging over the typical set can accurately estimate expectations, whereas naive quadrature requires N^D grid points and is unlikely to intersect the narrow contributing region.Source: §1.4 · pp. 7–8
• Markov chain Monte Carlo offers a generic strategy for quantifying the typical set of arbitrary target distributions after identifying where expectation contributions concentrate.Understanding the target region is only the first step; constructing an algorithm that can explore it remains a separate problem.Source: §1.4 · p. 8

2. MARKOV CHAIN MONTE CARLO

• MCMC uses a Markov chain to stochastically explore the typical set, but finite-time chains may fail to explore it adequately.A properly designed chain eventually explores any distribution given sufficient resources, while practical analyses require reliable behavior within finite computation.Source: §2 · pp. 8–9

2.1 Estimating Expectations with Markov Chains

• A target-preserving Markov transition eventually drives chains from any initial point into and across the typical set, enabling expectation estimates from their history.As chains run longer, their estimators converge, up to technicalities, to the true expectations.Source: §2.1 · pp. 10–11, Figure 6 · p. 10

2.2 Ideal Behavior

• Under ideal conditions, chains pass through convergence, initial typical-set exploration, and refinement phases, with estimator bias falling rapidly before precision improves more slowly.The effective sample size measures the number of exact samples with equivalent precision, or equivalently the chain’s effective sojourns across the typical set.Source: §2.2 · p. 11, Figure 7 · p. 12

2.3 Pathological Behavior

• Pathological high-curvature regions can be missed, causing chains to become stuck and finite-time estimators to oscillate around the true expectations with substantial bias.Geometric ergodicity formally guarantees ideal behavior, but is difficult to verify; split R-hat near 1 is an empirical diagnostic for suspicious chain inconsistency.Source: §2.3 · pp. 13, 15

2.4 The Metropolis-Hastings Algorithm

• Random Walk Metropolis combines stochastic proposals with acceptance correction, but its simplicity hides poor scaling: high-dimensional chains move slowly, yielding large autocorrelations and imprecise estimators.Large proposals are usually rejected outside the typical set, whereas small proposals improve acceptance but barely move the chain; Metropolis-Hastings is attributed to Metropolis et al. (1953) and Hastings (1970).Source: §2.4 · pp. 15–16

3. THE FOUNDATIONS OF HAMILTONIAN MONTE CARLO

• Random-walk proposals fail in high dimensions because large jumps leave the typical set while small accepted jumps explore only locally, motivating geometry-aware transitions.Effective transitions should follow high-probability contours coherently into unexplored regions rather than diffuse around the initial point.Source: §3 · pp. 16–17, Figure 10 · p. 17

3.1 Informing Effective Markov Transitions

• Hamiltonian Monte Carlo uses auxiliary momentum and differential geometry to twist gradient information into parameterization-invariant directions aligned with the typical set.The gradient alone points toward the mode; adding appropriately structured momentum enables conservative dynamics and stable orbital trajectories.Source: §3.1 · pp. 19, 20, 22

3.2 Phase Space and Hamilton’s Equations

• The canonical phase-space distribution preserves the target distribution under marginalization, while Hamilton’s equations generate rapidly moving trajectories constrained to its typical set.Momentum variables expand the D-dimensional parameter space to 2D-dimensional phase space, and projecting phase-space trajectories recovers efficient exploration of the target typical set.Source: §3.2 · pp. 23, 25

3.3 The Idealized Hamiltonian Markov Transition

• The idealized Hamiltonian Markov transition samples momentum, integrates Hamiltonian dynamics, and projects back to parameter space, rapidly exploring the target while remaining in its typical set.Starting from the target typical set, momentum sampling lifts the state into phase space; the resulting trajectories move away from previously explored neighborhoods, and the composed transition rapidly approaches the typical set from the tails.Source: §3.3 · p. 26

4. EFFICIENT HAMILTONIAN MONTE CARLO

• Hamiltonian Monte Carlo performance can vary dramatically because kinetic energy and integration time define infinitely many transitions, making their tuning essential.Well-chosen degrees of freedom support challenging high-dimensional problems, whereas poor choices can severely degrade performance.Source: §4 · p. 26

4.1 The Natural Geometry of Phase Space

• The phase space foliates into Hamiltonian energy level sets, so each transition alternates between deterministic exploration within a level set and stochastic movement across level sets.This decomposition enables separate efficiency criteria for trajectory exploration and energy diffusion.Source: §4.1 · pp. 27–28, Figure 22 · p. 29
• Efficient tuning requires regular energy level sets for rapid trajectory exploration and an energy transition distribution that matches the marginal energy distribution for rapid diffusion.The analysis constructs criteria for identifying effective choices of Hamiltonian degrees of freedom, including for complex targets.Source: §4.1 · pp. 28, 30

4.2.1 Euclidean-Gaussian Kinetic Energies

• Euclidean-Gaussian kinetic energies are optimized when the inverse Euclidean metric resembles the target covariance, which de-correlates the target and makes energy level sets more uniform.An empirical covariance estimate can be adapted during warm-up to obtain a near-optimal metric, but global rescaling cannot remove spatially varying curvature unless the target is exactly Gaussian.Source: §4.2.1 · p. 31, §4.2.2 · p. 31

4.2.2 Riemannian-Gaussian Kinetic Energies

• Riemannian-Gaussian kinetic energies use position-dependent metrics to correct local correlations, while the metric’s log-determinant can improve energy-transition matching for complex targets.A Hessian-like metric can produce highly uniform level sets, but SoftAbs regularization remains heuristic and its formalization is an active research topic (Holmes, Rubinstein-Salzedo and Seiler, 2014); non-Gaussian kinetic energies tend to perform poorly in high dimensions.Source: §4.2.2 · p. 32, §4.2.3 · p. 32

4.3 Optimizing the Choice of Integration Time

• Optimal integration time balances rapid early trajectory exploration against diminishing returns, and it must generally vary across trajectories because it scales with level-set energy.For heavy-tailed targets, β < 2, optimal integration times grow toward the tails; dynamic termination such as the No-U-Turn criterion (Betancourt, 2013b) is therefore needed, while exhaustive criteria can be more robust but require careful tuning (Betancourt, 2016a).Source: §4.3 · pp. 33–35

Numerical integration

• Numerical integration is the main implementation obstacle because Hamilton’s equations are rarely solvable exactly, and inaccuracies can quickly compromise a well-tuned Hamiltonian transition.The effectiveness of the implementation depends on how accurately the ordinary differential equations on phase space are solved.Source: §5 · p. 36
• Generic numerical integrators accumulate coherent error over long trajectories, driving them away from the true trajectory and typical set, especially as phase-space dimension increases.Their drift limits them to short trajectories that explore energy level sets inefficiently; Figure 27 illustrates the departure from exact solutions.Source: §5 · p. 36, Figure 27 · p. 36

Symplectic integrators

• Symplectic integrators use phase-space geometry to enable high-performance Hamiltonian Monte Carlo implementations despite numerical integration errors.The section presents their practical properties, error correction, and guidance for choosing an integrator for a given problem.Source: §5 · p. 36
• Symplectic integrators produce incompressible numerical trajectories whose errors oscillate around the exact energy level set rather than accumulating coherently over longer integrations.This geometric property makes them robust to drift and supports their use for high-performance implementations.Source: Figure 28 · p. 37

5.1 Symplectic Integrators

• Symplectic integrators preserve phase-space volume and keep numerical trajectories oscillating near exact energy levels, but unresolved high curvature can cause detectable divergence.These properties explain both long-time stability and a practical failure mode that rapidly drives trajectories toward infinite energies.Source: §5.1 · pp. 37–38, Figure 29 · p. 38

5.2 Correcting for Symplectic Integrator Error

• A momentum flip is essential because unreversed deterministic trajectories are nonreversible, making ordinary Metropolis-Hastings proposals always reject.Proposing all trajectory points remains reversible but can be suboptimal when a selected state has large numerical error.Source: Figure 30 · p. 39, §5.2 · pp. 39–40
• Metropolis-Hastings correction makes symplectic-integrator transitions exact, with acceptance degrading only negligibly as target dimension increases.A momentum flip restores reversibility, enabling a valid correction for numerical error.Source: §5.2 · pp. 39–40, Figure 31 · p. 40

5.3 Optimal Choice of Symplectic Integrator

• Integrator choice balances accuracy against cost: smaller step sizes improve exploration but increase L = T/ϵ, while higher order increases both accuracy and expense.Comparing exact and modified Hamiltonian level sets yields cost–acceptance bounds that support adaptation toward average acceptance between 0.6 and 0.8.Source: §5.3 · p. 41, Figure 33 · p. 43
• The formal cost–acceptance analysis applies directly to simple fixed-time implementations, while corresponding guarantees for more sophisticated HMC methods need further work.This limits how broadly the tuning results can currently be generalized.Source: §5.3 · pp. 41, 43
• Higher-order integrators may outperform leapfrog in higher dimensions, but poorly chosen extra parameters can eliminate their gains and require new adaptation methods.The cited work motivates higher-order methods, while the appendix identifies parameter interactions and automated adaptation as unresolved next steps.Source: §5.3 · p. 43

6. THE ROBUSTNESS OF HAMILTONIAN MONTE CARLO

• HMC implementations are geometrically ergodic for a target class significantly larger than that of Random Walk Metropolis, while gradients support robust high-dimensional exploration.The analysis also identifies target-distribution pathologies that obstruct geometric ergodicity and motivate HMC-specific diagnostics.Source: §6 · p. 44

6.1 Diagnosing Poorly-Chosen Kinetic Energies

• Poorly chosen kinetic energies can make energy transitions too slow to explore heavy-tailed targets completely; an energy Bayesian fraction of missing information below 0.3 has proven problematic empirically.Visualizing the marginal energy and energy-transition densities, including across diffuse-start chains, reveals inefficient exploration of energy level sets.Source: §6.1 · p. 44, Figure 34 · p. 45

6.2 Diagnosing Regions of High Curvature

• High-curvature neighborhoods can cause incomplete exploration and biased estimators, while divergent numerical trajectories provide an immediate diagnostic of these pathological regions.Mild divergences may be removed by reducing step size; persistent divergences require regularization such as stronger priors or alternative parameterizations.Source: §6.2 · p. 45

6.3 Limitations of Diagnostics

• Energy and divergence diagnostics identify many pathologies but are necessary, not sufficient, for estimator validity because subtler failures may remain in increasingly complex targets.Divergences may reflect high-curvature neighborhoods that frustrate both geometric ergodicity and symplectic integration.Source: §6.3 · pp. 45–46, Figure 35 · p. 46
• Combining HMC-specific diagnostics with split R̂ from multiple parallel chains makes HMC a uniquely robust method in practice.The claim concerns practical robustness despite the diagnostics’ inability to guarantee estimator validity.Source: §6.3 · pp. 45–46

• Hamiltonian Monte Carlo uses typical-set geometry to generate coherent exploration, improving computational efficiency and estimator-validity guarantees over other Markov chain Monte Carlo algorithms.Geometric analysis also supports principled automatic construction of optimal implementations, exemplified by the success of Stan.Source: §7 · p. 46
• The typical set concentrates near a complex surface in parameter space, making differential geometry a natural framework for quantifying and exploiting it.Source: §7 · p. 46
• Future work includes improving geometric ergodicity analysis, refining implementation techniques, and using non-Gaussian or non-Euclidean kinetic energies.Source: §7 · pp. 46–47
• Differential-geometric techniques may generalize Hamiltonian Monte Carlo principles to harder computational problems and parameter spaces, including discrete, tree, and Hilbert spaces.Geometric methods have already generalized coherent exploration to thermodynamic algorithms in Adiabatic Monte Carlo (Betancourt, 2014).Source: §7 · p. 47

• HMC alternates exploration within an energy level set with momentum resampling to generate efficient exploration of the canonical distribution.The method leverages the decomposition into a microcanonical distribution over energy level sets and a marginal energy distribution.Source: Appendix A · p. 47
• Hamiltonian trajectories preserve energy and are typically ergodic, so longer integration explores the entire energy level set and approaches the microcanonical distribution.Source: Appendix A · p. 47
• Because analytic Hamiltonian trajectories are generally unavailable, practical implementations use symplectic numerical integrators whose errors can be corrected probabilistically.Only a few exceptional systems permit analytic trajectory generation; symplectic integrators provide small numerical errors with straightforward corrections.Source: Appendix A · pp. 47–48
• The appendix develops implementations for fixed-length static trajectories and dynamically growing trajectories, then discusses their relationship to NUTS and Stan.Source: Appendix A · p. 48

A.1 Notation

• A phase-space point is denoted z = (q, p), and a symplectic integrator generates discrete trajectories by integrating forward or backward by ϵ.The integrator reduces phase space to points reachable through these incremental numerical integrations.Source: §A.1 · p. 48
• From an initial point z0, zn denotes the point reached after n forward integrations, while z−n denotes the point reached after n backward integrations.For L = 3, example trajectories include {z0, z1, z2}, {z−1, z0, z1}, and {z−2, z−1, z0}.Source: §A.1 · p. 48
• The notation T^L_z denotes all length-L symplectic-integrator trajectories containing z, while T^L_{z,z′} denotes those containing both z and z′.These sets organize numerical trajectories according to the points they contain.Source: §A.1 · p. 48

A.2 Static Implementations

• Canonical-density weighting exactly corrects the sampling error caused by non-energy-preserving symplectic integrators.Uniform sampling along a numerical trajectory does not produce the canonical distribution.Source: §A.2 · p. 48
• The transition is constructed sequentially by sampling a trajectory containing the initial point and then sampling a state from that trajectory.This yields the joint transition T(z′ | t) T(t | z).Source: §A.2 · p. 48
• For static implementations, only length-L trajectories containing the initial state z are considered, so T(t | z) is supported on T_L.Source: §A.2.1 · p. 48

• Trajectory transitions must be reversible: the probability of selecting a trajectory must be the same from every state on that trajectory.Source: §A.2.1 · p. 49
• Uniformly sampling all trajectories containing the initial point is the simplest practical way to satisfy the required condition.Source: §A.2.1 · p. 49
• After sampling a numerical trajectory, any process with transition probabilities proportional to canonical densities can generate an exact canonical-distribution sample.Source: §A.2.2 · p. 49

• Exact trajectories would yield uniform transitions over states on the initial energy level, whereas numerical error disfavors moves toward higher energies.The non-uniformity arises from deviations in the numerical trajectory’s energy.Source: §A.2.2 · p. 49
• A slice sampler can recover the desired transitions by introducing u ∼ U(0, π(z)) and sampling uniformly among trajectory states with π(z′) > u.Source: §A.2.2 · p. 49
• Direct multinomial sampling over trajectory states with the desired probabilities avoids auxiliary variables and is straightforward to implement.The Rao–Blackwell Theorem establishes that omitting auxiliary random variables improves performance.Source: §A.2.2 · p. 49
• Symplectic integrators support high practical performance by exactly preserving phase-space volume and enabling exact canonical samples for canonical-density-preserving transitions.Source: §A.2.3 · p. 50

A.2.3 Invariance of the Canonical Distribution

• Explicit calculation shows that the stated trajectory properties suffice to preserve invariance of the canonical distribution.Source: §A.2.3 · p. 50

A.3.1 Uniform Progressive Sampling

• Progressive sampling avoids storing entire trajectories by mixing old and newly integrated components, while recursive doubling builds longer trajectories and samples concurrently.The mixture transition combines component distributions, and binary-tree structure supports recursive multiplicative expansion.Source: §A.3 · p. 50, §A.3.1 · pp. 50–51

A.3.2 Biased Progressive Sampling

• Biased progressive sampling favors newly added trajectory components, yet averaging over possible initial states cancels the induced repulsion and preserves the restricted canonical distribution.Source: §A.3.2 · p. 54

A.4 Dynamic Implementations

• Dynamic implementations adapt integration time across energy level sets because static schemes cannot efficiently explore levels requiring substantially different trajectory lengths.They expand a trajectory until a selected termination criterion indicates sufficient exploration, while preserving reversibility remains essential.Source: §A.4 · p. 55

A.4.1 Maintaining Reversibility

• Reversible dynamic expansion rejects trajectories containing prematurely terminating sub-trajectories, with multiplicative expansion reducing checking costs relative to additive expansion.Additive checking is quadratic in final trajectory length, whereas multiplicative checking is quadratic only in its logarithm and can proceed concurrently.Source: §A.4.1 · p. 55, Figure 39 · p. 56, Figure 40 · p. 56

A.4.2 Dynamic Termination Criteria

• The original No-U-Turn criterion stops when endpoint momenta become anti-aligned with the boundary-connecting direction, typically after the trajectory spans the energy-level-set width.Source: §A.4.2 · pp. 56–57, Figure 42 · p. 57
• The generalized No-U-Turn criterion applies on Riemannian manifolds by replacing Euclidean point subtraction with an integral approximated from numerical-trajectory momenta.The approximation produces almost identical results to the original criterion because symplectic integrators are highly accurate (Betancourt, 2013b).Source: §A.4.2 · pp. 57–58

A.5 The No-U-Turn Sampler and the Current State of Stan

• Stan’s Hamiltonian Monte Carlo uses generalized No-U-Turn termination and multinomial trajectory sampling, significantly improving overall performance over slice sampling.The original No-U-Turn sampler used multiplicative expansion, slice sampling, and uniform and biased progressive sampling.Source: §A.5 · pp. 58–59