A Conceptual Introduction to Hamiltonian Monte Carlo
arxiv.org
arxiv.org
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!
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.
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.
Not quite. You explore the level sets of the augmented system of kinetic-plus-potential energy. The use of a Gaussian distribution as auxiliary variable is arbitrary, so it's not surprising that it's suboptimal.
That the gradient information is needed comes from physical intuition. I roll a marble with along an inclined plane with velocity V. What determines the trajectory of the marble? The gradient of the plane. Now since differentiable functions look locally like planes, that's more or less all you need, at least instantaneously.
I'm sure you know more about this topic than I do, but you'll have to work quite a bit harder to convince me that the hockey puck intuition is "incredibly shallow".
Even in your marble system, for example, the gradient of the plane does not affect the position of the marble directly. The gradient affects the velocity which then moves the marble. If the interactions aren't just right then you no longer have a physical system with energy conservation and all the other nice properties one might expect.
In both the physical system of the marble and the pseudo-physical system created by HMC there is a deeper mathematical structure that relates gradients to motion to probability in exactly the right way. If that structure is compromised, say in a naive implementation of HMC, then all of the magic is lost and you don't get the performance that you might expect.
It may sound like I'm being pedantic here, but it's only be understanding and respecting the underlying math far beneath the "hockey puck" description that we've been able to take HMC to the next level, especially in tools like Stan, with faster and more robust performance as well as sensitive diagnostics of failures.
Take a sample from an N-dimensional standard Gaussian. The square of its Its distance from the origin is x_1^2 + x_2^2 + ... + x_N^2. Each component of this sum has mean 1, so the law of large numbers says for large N, the sum will be about equal to N.
On other words, I take a sample from a high-dimensional Gaussian and with very high probability, it will live a distance sqrt(N) from the origin. Which is counterintuitive as you point out, since the most likely value is the 0 vector.
Edit: I'm not sure if your comment contradicts or complements what I wrote. As I said, taking a hypersphere of radius 1 the "inner" hypersphere of radius 0.5 accounts for just 1/2^n of the hypervolume. Would you say that the "outer" shell is a neighbourhood? How is that defined? What is the advantage of doing so instead of using the whole (only marginally bigger) hypersphere?
Defined in this way, you could of course leave the origin in the typical set. But then there would be another, smaller, typical set that would do the job just as well.
> A grid of length N distributed uniformly in a D-dimensional parameter space requires N^D points and hence N^D evaluations of the integrand. Unless N is incredibly large, however, it is unlikely that any of these points will intersect the narrow typical set, and the exponentially-growing cost of averaging over the grid yields worse and worse approximations to expectations.
In fact any point within the hypersphere is as good (or better) as any point in that narrow typical set not including the mode. The problem is not missing the "narrow typical set", is missing the whole hypersphere (because the volume of the hypersphere is small). I fail to see where in that introduction the fact that the typical set excludes the mode of the distribution makes any difference.
The typical set is a vaguely defined notation as you cannot define an explicit neighborhood without defining some explicit probability threshold which is often done in information theory of discrete systems. The point here, however, is not to define an explicit neighborhood but rather to motivate that the whichever neighborhood of high-probability mass you would define would not near the mode and instead extend out into a compact neighborhood around the mode.
You are welcome to define a neighborhood as the convex hull of the typical set, which would then include the mode (subject to some technical conditions). That wouldn't add any appreciable probability mass, however, but would introduce a whole bunch of space that your statistical algorithm would have to explore at additional computational cost.
No, one evaluation near the mode is more efficient than one evaluation at the boundary (because the density is higher).
I agree that the region immediately around the mode is naturally "avoided" (because the volume is small). But your paper makes it look like we increase efficiency by explicitely avoiding it and concentrating somewhere else. That's what I found confusing.
Incorrect in general. Firstly, one evaluation anywhere does not yield any reasonably accurate estimate of expectations. We'll always need an ensemble of evaluations, in which case the relevant question really is "how should I distribute these evaluations". And for any scalable algorithm none of them will be near the mode.
This is often hard to grok because people implicit fall back to the example of a gaussian density where the mode and the Hessian at the mode fully characterize the density which can then be used to compute analytic integrals. But for a general target distribution we do not have any of that structure and instead have to consider general computational strategies.
"I agree that the region immediately around the mode is naturally "avoided" (because the volume is small). But your paper makes it look like we increase efficiency by explicitely avoiding it and concentrating somewhere else. That's what I found confusing."
In high-dimensions the probability mass of any well-behaved probability distribution concentrates in (or _around_ if you want to acknowledge the fuzziness) the typical set. Hence accurate estimation of expectations requires quantifying the typical set. Any evaluation outside of the typical set is wasted because it offers increasingly negligible contributions to the integrals -- a few additional evaluations may not drive the cost up appreciably but they will still be wasted.
So it's not that the only thing that matters is avoiding the mode. Rather what matters is that, contrary to many's intuitions, the neighborhood around the mode does inform expectation values and so exploring that neighborhood is insufficient for estimating expectations. That then motivates the question of what neighborhoods do matter, which is answered by concentration of measure and the existence of the typical set.
And in practice, we don't actually do any of this explicitly. Instead we construct algorithms that somehow quantify probability mass (MCMC, VB, etc) and they will implicitly avoiding the mode and end up working with the typical set.
Why not distributing them according to the probability distribution? If I sample one million points and one happens to be near the mode, this doesn't mean any additional cost. This is no worse than sampling one million points none of which happens to be near the mode.
> Any evaluation outside of the typical set is wasted because it offers increasingly negligible contributions to the integrals
This is not a reason to exclude the mode from the typical set, because one evaluation at the mode offers a larger contribution to the integrals than one evaluation elsewhere.
> Instead we construct algorithms that somehow quantify probability mass (MCMC, VB, etc) and they will implicitly avoiding the mode and end up working with the typical set.
The algorithms don't have to avoid the mode, they just have to sample from it fairly (not much, but more than from any other region of the same size in the rest of the typical set).
By construct samples from a distribution will concentrate in neighborhoods of high probability mass, and hence the typical set.
Excluding the region with highest probability density from the typical set is a bit like saying that the population of New York is concentrated outside NYC because most people lives elsewhere.
https://youtu.be/DJ0c7Bm5Djk?t=16810
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/, which works well on many problems.
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.
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.
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.
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.
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.
Suppose I'll have to try it now.