Hamiltonian Monte Carlo

Short description: Sampling algorithm

Hamiltonian Monte Carlo sampling a two-dimensional probability distribution

The Hamiltonian Monte Carlo algorithm (originally known as hybrid Monte Carlo) is a Markov chain Monte Carlo method for obtaining a sequence of random samples whose distribution converges to a target probability distribution that is difficult to sample directly. This sequence can be used to estimate integrals of the target distribution, such as expected values and moments.

Hamiltonian Monte Carlo corresponds to an instance of the Metropolis–Hastings algorithm, with a Hamiltonian dynamics evolution simulated using a time-reversible and volume-preserving numerical integrator (typically the leapfrog integrator) to propose a move to a new point in the state space. Compared to using a Gaussian random walk proposal distribution in the Metropolis–Hastings algorithm, Hamiltonian Monte Carlo reduces the correlation between successive sampled states by proposing moves to distant states which maintain a high probability of acceptance due to the approximate energy conserving properties of the simulated Hamiltonian dynamic when using a symplectic integrator.[citation needed] The reduced correlation means fewer Markov chain samples are needed to approximate integrals with respect to the target probability distribution for a given Monte Carlo error.

The algorithm was originally proposed by Simon Duane, Anthony Kennedy, Brian Pendleton and Duncan Roweth in 1987 for calculations in lattice quantum chromodynamics.[1] In 1996, Radford M. Neal showed how the method could be used for a broader class of statistical problems, in particular artificial neural networks.[2] But the burden of having to provide the algorithm with gradients of the model graph delayed its wider adoption in statistics and other quantitative disciplines, until in the mid-2010s the developers of Stan implemented HMC in combination with automatic differentiation.[3]

## Contents

* 1 Algorithm
* 2 No U-Turn Sampler
* 3 See also
* 4 References
* 5 Further reading
* 6 External links

## Algorithm

Suppose the target distribution to sample is f(𝐱) for 𝐱∈ℝd (d≥1) and a chain of samples 𝐗0,𝐗1,𝐗2,… is required.

The Hamilton's equations are

dxidt=∂H∂pianddpidt=−∂H∂xi

where xi and pi are the ith component of the position and momentum vector respectively and H is the Hamiltonian. Let M be a mass matrix which is symmetric and positive definite, then the Hamiltonian is

H(𝐱,𝐩)=U(𝐱)+12𝐩TM−1𝐩

where U(𝐱) is the potential energy. The potential energy for a target is given as

U(𝐱)=−ln⁡f(𝐱)

which comes from the Boltzmann's factor. Note that the Hamiltonian H is dimensionless in this formulation because the exponential probability weight exp⁡(−H) has to be well defined. For example, in simulations at finite temperature T the factor kBT (with the Boltzmann constant kB) is directly absorbed into U and M.

The algorithm requires a positive integer for number of leapfrog steps L and a positive number for the step size Δt. Suppose the chain is at 𝐗n=𝐱n. Let 𝐱n(0)=𝐱n. First, a random Gaussian momentum 𝐩n(0) is drawn from N(𝟎,M). Next, the particle will run under Hamiltonian dynamics for time LΔt, this is done by solving the Hamilton's equations numerically using the leapfrog algorithm. The position and momentum vectors after time Δt using the leapfrog algorithm are:[4]

𝐩n(t+Δt2)=𝐩n(t)−Δt2∇U(𝐱)|𝐱=𝐱n(t)
𝐱n(t+Δt)=𝐱n(t)+ΔtM−1𝐩n(t+Δt2)
𝐩n(t+Δt)=𝐩n(t+Δt2)−Δt2∇U(𝐱)|𝐱=𝐱n(t+Δt)

These equations are to be applied to 𝐱n(0) and 𝐩n(0) L times to obtain 𝐱n(LΔt) and 𝐩n(LΔt).

The leapfrog algorithm is an approximate solution to the motion of non-interacting classical particles. If exact, the solution will never change the initial randomly-generated energy distribution, as energy is conserved for each particle in the presence of a classical potential energy field. In order to reach a thermodynamic equilibrium distribution, particles must have some sort of interaction with, for example, a surrounding heat bath, so that the entire system can take on different energies with probabilities according to the Boltzmann distribution.

One way to move the system towards a thermodynamic equilibrium distribution is to change the state of the particles using the Metropolis–Hastings algorithm. So first, one applies the leapfrog step, then a Metropolis-Hastings step.

The transition from 𝐗n=𝐱n to 𝐗n+1 is

𝐗n+1|𝐗n=𝐱n={𝐱n(LΔt)with probability α(𝐱n(0),𝐱n(LΔt))𝐱n(0)otherwise

where

α(𝐱n(0),𝐱n(LΔt))=min(1,exp⁡[−H(𝐱n(LΔt),𝐩n(LΔt))]exp⁡[−H(𝐱n(0),𝐩n(0))]).

A full update consists of first randomly sampling the momenta 𝐩 (independently of the previous iterations), then integrating the equations of motion (e.g. with leapfrog), and finally obtaining the new configuration from the Metropolis-Hastings accept/reject step. This updating mechanism is repeated to obtain 𝐗n+1,𝐗n+2,𝐗n+3,….

## No U-Turn Sampler

The No U-Turn Sampler (NUTS)[5] is an extension by controlling the number of steps L automatically. Tuning L is critical. For example, in the one dimensional N(0,1/k) case, the potential is U(x)=kx2/2 which corresponds to the potential of a simple harmonic oscillator. For L too large, the particle will oscillate and thus waste computational time. For L too small, the particle will behave like a random walk.

Loosely, NUTS runs the Hamiltonian dynamics both forwards and backwards in time randomly until a U-Turn condition is satisfied. When that happens, a random point from the path is chosen for the MCMC sample and the process is repeated from that new point.

In detail, a binary tree is constructed to trace the path of the leap frog steps. To produce a MCMC sample, an iterative procedure is conducted. A slice variable Un∼Uniform(0,exp⁡(−H[𝐱n(0),𝐩n(0)])) is sampled. Let 𝐱n+ and 𝐩n+ be the position and momentum of the forward particle respectively. Similarly, 𝐱n− and 𝐩n− for the backward particle. In each iteration, the binary tree selects at random uniformly to move the forward particle forwards in time or the backward particle backwards in time. Also for each iteration, the number of leap frog steps increase by a factor of 2. For example, in the first iteration, the forward particle moves forwards in time using 1 leap frog step. In the next iteration, the backward particle moves backwards in time using 2 leap frog steps.

The iterative procedure continues until the U-Turn condition is met, that is

(𝐱n+−𝐱n−)⋅𝐩n−<0or.(𝐱n+−𝐱n−)⋅𝐩n+<0

or when the Hamiltonian becomes inaccurate

exp⁡[−H(𝐱n+,𝐩n+)+δ]<Un

exp⁡[−H(𝐱n−,𝐩n−)+δ]<Un

where, for example, δ=1000.

Once the U-Turn condition is met, the next MCMC sample, 𝐱n+1, is obtained by sampling uniformly the leap frog path traced out by the binary tree {𝐱n−,…,𝐱n(−Δt),𝐱n(0),𝐱n(Δt),…,𝐱n+} which satisfies

Un<exp⁡[−H(x𝐧+𝟏,p𝐧+𝟏)]

This is usually satisfied if the remaining HMC parameters are sensible.[citation needed]

## See also

* Dynamic Monte Carlo method
* Software for Monte Carlo molecular modeling
* Stan, a probabilistic programing language implementing HMC.
* PyMC, a probabilistic programming language implementing HMC.
* Metropolis-adjusted Langevin algorithm
* Turing.jl, a probabilistic programming language in Julia implementing HMC.

## References

1. ↑ Duane, Simon; Kennedy, Anthony D.; Pendleton, Brian J.; Roweth, Duncan (1987). "Hybrid Monte Carlo". Physics Letters B 195 (2): 216–222. doi:10.1016/0370-2693(87)91197-X. Bibcode: 1987PhLB..195..216D.
2. ↑ Neal, Radford M. (1996). "Monte Carlo Implementation". Bayesian Learning for Neural Networks. Lecture Notes in Statistics. 118. Springer. pp. 55–98. doi:10.1007/978-1-4612-0745-0_3. ISBN 0-387-94724-8.
3. ↑ Gelman, Andrew; Lee, Daniel; Guo, Jiqiang (2015). "Stan: A Probabilistic Programming Language for Bayesian Inference and Optimization". Journal of Educational and Behavioral Statistics 40 (5): 530–543. doi:10.3102/1076998615606113.
4. ↑ Betancourt, Michael (2018-07-15). "A Conceptual Introduction to Hamiltonian Monte Carlo". arXiv:1701.02434 [stat.ME].
5. ↑ Hoffman, Matthew D; Gelman, Andrew (2014). "The No-U-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo". Journal of Machine Learning Research 15 (1): 1593–1623. http://jmlr.org/papers/v15/hoffman14a.html. Retrieved 2024-03-28.

## Further reading

* Betancourt, Michael; Girolami, Mark (2015). "Hamiltonian Monte Carlo for Hierarchical Models". in Upadhyay, Satyanshu Kumar; Singh, Umesh; Dey, Dipak K. et al.. Current Trends in Bayesian Methodology with Applications. CRC Press. pp. 79–101. ISBN 978-1-4822-3511-1.
* Betancourt, Michael (2018). "A Conceptual Introduction to Hamiltonian Monte Carlo". arXiv:1701.02434 [stat.ME].
* Barbu, Adrian; Zhu, Song-Chun (2020). "Hamiltonian and Langevin Monte Carlo". Monte Carlo Methods. Singapore: Springer. pp. 281–326. ISBN 978-981-13-2970-8.
* Neal, Radford M (2011). "MCMC Using Hamiltonian Dynamics". in Steve Brooks. Handbook of Markov Chain Monte Carlo. Chapman and Hall/CRC. ISBN 9781420079418. http://www.mcmchandbook.net/HandbookChapter5.pdf.

## External links

* Betancourt, Michael (17 November 2014). "Efficient Bayesian inference with Hamiltonian Monte Carlo". MLSS Iceland 2014. https://www.youtube.com/watch?v=pHsuIaPbNbY&list=PLqdbxUnkqOw2nKn7VxYqIrKWcqRkQYOsF&index=11.
* McElreath, Richard (25 January 2022). "Markov chain Monte Carlo". Statistical Rethinking 2022. https://www.youtube.com/watch?v=Qqz5AJjyugM&list=PLDcUM9US4XdMROZ57-OIRtIK0aOynbgZN&index=9.
* Hamiltonian Monte Carlo from scratch
* Optimization and Monte Carlo Methods

0.00

(0 votes)

Original source: https://en.wikipedia.org/wiki/Hamiltonian Monte Carlo. Read more

|

Retrieved from "https://handwiki.org/wiki/index.php?title=Hamiltonian_Monte_Carlo&oldid=3980594"