7 ms·

A Conceptual Introduction to Hamiltonian Monte Carlo

155 points by gwern 10y ago

tmalsburg2 10y ago
Also worth checking out is the author's talk at this year's StanCon: https://youtu.be/DJ0c7Bm5Djk?t=16810 Title: Everything you should have learned about Markov Chain Monte Carlo.
- leblancfg 10y ago
  Great talk, thanks for the link!
- betanalpha 10y ago
  There are also talks focusing on HMC available online, the most recent of which are https://www.youtube.com/watch?v=VnNdhsm0rJQ https://www.youtube.com/watch?v=VnNdhsm0rJQ (more introductory) and https://icerm.brown.edu/video_archive/#/play/1107 https://icerm.brown.edu/video_archive/#/play/1107 (more mathematically formal).
xioxox 10y ago
As far as I know, you need to be able to compute the gradient of the function being explored to use this method. Is that correct? Most of my likelihood functions I use are nasty and complex. I've mostly been using an affine-invariant sampler, emcee, http://dan.iel.fm/emcee/current/ http://dan.iel.fm/emcee/current/, which works well on many problems.
- plafl 10y ago
  I have used the NUTS sampler, via PyMC. If your function is smooth you are OK, even if it's complex, because both PyMC and I think Stan too will auto differentiate your function. PyMC in particular uses theano under the hood.
- contravariant 10y ago
  You need the gradient of the log-likelihood in all Hamiltonian MCMC implementations I've seen so far (you don't need it in classical MCMC though). However I'm reasonably sure that the method will still work even if you use the gradient of some approximate log-likelihood (the proof only requires the dynamics to be reversible and preserve volume, although the closer you can get to an exact solution the better so use this information at your own risk).
  - yummyfajitas 10y ago
    I don't think this is right - small deviations from the true gradient will result in large deviations along the Hamiltonian path. The magic trick of HMC is that you can safely make large jumps in state space because you are following the right path - doesn't feel safe to break this assumption. But I'm not an expert on HMC. My intuition comes mainly from my past career as a physicist. So if someone knows why it would be safe I'd love to hear.
    - imurray 10y ago
      Implementations of HMC are approximating ideal Hamiltonian dynamics by discretizing time. The "symplectic" integrators that are used have the property that they correspond to simulating continuous-time Hamiltonian dynamics exactly on some other energy surface, which is not too far from the one you had in mind. So, in some sense, an approximation like the grandparent had in mind is the usual case. You are right that a bad approximation of the surface could lead to proposals that get rejected a lot. But at least the dynamics are likely to be stable. It could be worth trying.
    - contravariant 10y ago
      It will result in changes in total 'energy', but that's fine if you account for this in the acceptance probability.
    - betanalpha 10y ago
      This intuition is correct -- symplectic integrators admit highly accurately energy conservation only when you have access to the exact gradients. When the gradients are inexact then you get bias in the integrator that pulls you away from the target energy level set and decreases the acceptance probabilities (or rather the transition probabilities away from the initial point if you using a more sophisticated scheme as advocated in the paper). Note that symplectic integrator trajectories stay near the target energy level set but they need not stay near the true trajectory itself. In general the two will diverge pretty quickly (a manifestation of Hamiltonian chaos) within the level sets, but that's not a problem as we don't really care how they explore the level sets, just that they're staying on them.
- imurray 10y ago
  Yes, you need the gradient. And despite what the autodiff people say, that might be hard and/or not very useful in practice, if your likelihood is based on running a lengthy astrophysical simulation. You can run HMC on any differentiable surrogate surface surface. It's only the accept/reject step that has to use your real model and data. You might consider that if you have a less hairy approximation to your model. Emcee is popular in astronomy, partly because of network effects, partly because it's a really nice package, and partly because the method's a great fit for some of the posteriors you have: often unimodal-ish even if the likelihood is nasty to compute, and often in not that many dimensions. It would fall flat on its face for something like a Bayesian neural network.
  - betanalpha 10y ago
    "And despite what the autodiff people say, that might be hard and/or not very useful in practice, if your likelihood is based on running a lengthy astrophysical simulation." Only if you treat the simulation as an un-interrogatable black box. All of those astro simulations are ODE or PDE systems of one kind or another which can be very cleanly autodiffed to calculate exact gradients. N-bodies are particularly well-suited to this approach. "You can run HMC on any differentiable surrogate surface surface. It's only the accept/reject step that has to use your real model and data. You might consider that if you have a less hairy approximation to your model." You _can_ but it won't work above O(10) dimensions. The problem is that the more dimension you have the more any surrogate surface (and its gradients) will drift away from the true surface and worse your acceptance probabilities will be. With exact gradients the optimal performance of HMC will drop negligibly with dimension (you need to get up O(2^100) or so before you really start to see problems), assuming the target distribution itself doesn't get more complex its dimension increases.
- contravariant 10y ago
  I've already responded regarding the need to calculate the gradient, but after reading up on the affine-invariant sampler I'm wondering if it wouldn't be an idea to use a hybrid setup where you have several chains running alongside like the affine-invariant sampler and use those to approximately solve the Hamiltonian dynamics... not sure what will happen, but it removes some of the 'arbitrarity' in the proposal of steps for the affine-invariant method, and actually uses the already computed values of the distribution, which should (in theory) be better. Suppose I'll have to try it now.
soVeryTired 10y ago
Radford Neal's expository paper [1] on HMC is a beautiful overview of the area, and honestly one of the most enlightening papers I've ever read. Conceptually, HMC really isn't that difficult. The broad idea is as follows: Borrowing from statistical physics, it's possible to link energy and probability as E(x) = exp(-P(x)). So you can interpret your MCMC target as an "energy landscape". Highpoints in the landscape correspond to low probabilities in the target distribution, and valleys to high-probability regions. Now add a source of kinetic energy (an auxiliary probability distribution), and think of your MCMC sampler as a hockey puck sliding frictionlessly through the energy landscape. Since kinetic + potential energy is conserved, it follows that you move through contours of your joint (target, auxiliary) distribution. That means you never have to reject a sample, regardless of how far the hockey puck travels. This interpretation also shows why you need to use gradient information. To know where the hockey puck is going to move next, you need to know how steep the walls of your energy landscape are: that is, you need the gradient! [1] https://arxiv.org/pdf/1206.1901.pdf https://arxiv.org/pdf/1206.1901.pdf
- sndean 10y ago
  > So you can interpret your MCMC target as an "energy landscape". Highpoints in the landscape correspond to low probabilities in the target distribution, and valleys to high-probability regions. Thinking of things in terms of energy landscapes was the only I could understand this stuff during grad school. Plus they had some of the coolest visualizations (e.g., Fig 2 from [1]) that kept me interested. [1] http://www.pnas.org/content/101/41/14766.full.pdf http://www.pnas.org/content/101/41/14766.full.pdf
- betanalpha 10y ago
  So V(x) = exp(-pi(x)) and you use HMC to explore those potential energy level sets? But what happens when you change coordinates and get a different density, pi(x), a different potential energy, V(x), and hence different level sets? And why does the gradient help as it points transverse to the level sets and not along them which you would need to move in the right direction? Maybe the momenta fixes that somehow? So you just add any IID gaussian variables to the target? You need special variables? Even so, why IID gaussians? The problem with the "hockey puck" intuition is that it is incredibly shallow and obfuscates the underlying principles that really make the method work. Unfortunately the false intuition also motivates bad tuning strategies and method extensions which have real practical consequences. Radford's review was really, really important but it was written when we really didn't know what was going on at a foundational level. The linked paper has the benefit of the last five years of research where we've figured those foundations out and consequently can identify the principled important to robust use of the method in practice.

155 points by gwern 10y ago

Also worth checking out is the author's talk at this year's StanCon: https://youtu.be/DJ0c7Bm5Djk?t=16810 Title: Everything you should have learned about Markov Chain Monte Carlo.

leblancfg 10y ago
Great talk, thanks for the link!
betanalpha 10y ago
There are also talks focusing on HMC available online, the most recent of which are https://www.youtube.com/watch?v=VnNdhsm0rJQ https://www.youtube.com/watch?v=VnNdhsm0rJQ (more introductory) and https://icerm.brown.edu/video_archive/#/play/1107 https://icerm.brown.edu/video_archive/#/play/1107 (more mathematically formal).

xioxox 10y ago

As far as I know, you need to be able to compute the gradient of the function being explored to use this method. Is that correct? Most of my likelihood functions I use are nasty and complex. I've mostly been using an affine-invariant sampler, emcee, http://dan.iel.fm/emcee/current/ http://dan.iel.fm/emcee/current/, which works well on many problems.

plafl 10y ago
I have used the NUTS sampler, via PyMC. If your function is smooth you are OK, even if it's complex, because both PyMC and I think Stan too will auto differentiate your function. PyMC in particular uses theano under the hood.
contravariant 10y ago
You need the gradient of the log-likelihood in all Hamiltonian MCMC implementations I've seen so far (you don't need it in classical MCMC though). However I'm reasonably sure that the method will still work even if you use the gradient of some approximate log-likelihood (the proof only requires the dynamics to be reversible and preserve volume, although the closer you can get to an exact solution the better so use this information at your own risk).
- yummyfajitas 10y ago
  I don't think this is right - small deviations from the true gradient will result in large deviations along the Hamiltonian path. The magic trick of HMC is that you can safely make large jumps in state space because you are following the right path - doesn't feel safe to break this assumption. But I'm not an expert on HMC. My intuition comes mainly from my past career as a physicist. So if someone knows why it would be safe I'd love to hear.
  - imurray 10y ago
    Implementations of HMC are approximating ideal Hamiltonian dynamics by discretizing time. The "symplectic" integrators that are used have the property that they correspond to simulating continuous-time Hamiltonian dynamics exactly on some other energy surface, which is not too far from the one you had in mind. So, in some sense, an approximation like the grandparent had in mind is the usual case. You are right that a bad approximation of the surface could lead to proposals that get rejected a lot. But at least the dynamics are likely to be stable. It could be worth trying.
  - contravariant 10y ago
    It will result in changes in total 'energy', but that's fine if you account for this in the acceptance probability.
  - betanalpha 10y ago
    This intuition is correct -- symplectic integrators admit highly accurately energy conservation only when you have access to the exact gradients. When the gradients are inexact then you get bias in the integrator that pulls you away from the target energy level set and decreases the acceptance probabilities (or rather the transition probabilities away from the initial point if you using a more sophisticated scheme as advocated in the paper). Note that symplectic integrator trajectories stay near the target energy level set but they need not stay near the true trajectory itself. In general the two will diverge pretty quickly (a manifestation of Hamiltonian chaos) within the level sets, but that's not a problem as we don't really care how they explore the level sets, just that they're staying on them.
imurray 10y ago
Yes, you need the gradient. And despite what the autodiff people say, that might be hard and/or not very useful in practice, if your likelihood is based on running a lengthy astrophysical simulation. You can run HMC on any differentiable surrogate surface surface. It's only the accept/reject step that has to use your real model and data. You might consider that if you have a less hairy approximation to your model. Emcee is popular in astronomy, partly because of network effects, partly because it's a really nice package, and partly because the method's a great fit for some of the posteriors you have: often unimodal-ish even if the likelihood is nasty to compute, and often in not that many dimensions. It would fall flat on its face for something like a Bayesian neural network.
- betanalpha 10y ago
  "And despite what the autodiff people say, that might be hard and/or not very useful in practice, if your likelihood is based on running a lengthy astrophysical simulation." Only if you treat the simulation as an un-interrogatable black box. All of those astro simulations are ODE or PDE systems of one kind or another which can be very cleanly autodiffed to calculate exact gradients. N-bodies are particularly well-suited to this approach. "You can run HMC on any differentiable surrogate surface surface. It's only the accept/reject step that has to use your real model and data. You might consider that if you have a less hairy approximation to your model." You _can_ but it won't work above O(10) dimensions. The problem is that the more dimension you have the more any surrogate surface (and its gradients) will drift away from the true surface and worse your acceptance probabilities will be. With exact gradients the optimal performance of HMC will drop negligibly with dimension (you need to get up O(2^100) or so before you really start to see problems), assuming the target distribution itself doesn't get more complex its dimension increases.
contravariant 10y ago
I've already responded regarding the need to calculate the gradient, but after reading up on the affine-invariant sampler I'm wondering if it wouldn't be an idea to use a hybrid setup where you have several chains running alongside like the affine-invariant sampler and use those to approximately solve the Hamiltonian dynamics... not sure what will happen, but it removes some of the 'arbitrarity' in the proposal of steps for the affine-invariant method, and actually uses the already computed values of the distribution, which should (in theory) be better. Suppose I'll have to try it now.

soVeryTired 10y ago

Radford Neal's expository paper [1] on HMC is a beautiful overview of the area, and honestly one of the most enlightening papers I've ever read. Conceptually, HMC really isn't that difficult. The broad idea is as follows: Borrowing from statistical physics, it's possible to link energy and probability as E(x) = exp(-P(x)). So you can interpret your MCMC target as an "energy landscape". Highpoints in the landscape correspond to low probabilities in the target distribution, and valleys to high-probability regions. Now add a source of kinetic energy (an auxiliary probability distribution), and think of your MCMC sampler as a hockey puck sliding frictionlessly through the energy landscape. Since kinetic + potential energy is conserved, it follows that you move through contours of your joint (target, auxiliary) distribution. That means you never have to reject a sample, regardless of how far the hockey puck travels. This interpretation also shows why you need to use gradient information. To know where the hockey puck is going to move next, you need to know how steep the walls of your energy landscape are: that is, you need the gradient! [1] https://arxiv.org/pdf/1206.1901.pdf https://arxiv.org/pdf/1206.1901.pdf

sndean 10y ago
> So you can interpret your MCMC target as an "energy landscape". Highpoints in the landscape correspond to low probabilities in the target distribution, and valleys to high-probability regions. Thinking of things in terms of energy landscapes was the only I could understand this stuff during grad school. Plus they had some of the coolest visualizations (e.g., Fig 2 from [1]) that kept me interested. [1] http://www.pnas.org/content/101/41/14766.full.pdf http://www.pnas.org/content/101/41/14766.full.pdf
betanalpha 10y ago
So V(x) = exp(-pi(x)) and you use HMC to explore those potential energy level sets? But what happens when you change coordinates and get a different density, pi(x), a different potential energy, V(x), and hence different level sets? And why does the gradient help as it points transverse to the level sets and not along them which you would need to move in the right direction? Maybe the momenta fixes that somehow? So you just add any IID gaussian variables to the target? You need special variables? Even so, why IID gaussians? The problem with the "hockey puck" intuition is that it is incredibly shallow and obfuscates the underlying principles that really make the method work. Unfortunately the false intuition also motivates bad tuning strategies and method extensions which have real practical consequences. Radford's review was really, really important but it was written when we really didn't know what was going on at a foundational level. The linked paper has the benefit of the last five years of research where we've figured those foundations out and consequently can identify the principled important to robust use of the method in practice.