Hacker News Reader

Comment by xioxox | Hacker News Reader

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/, which works well on many problems.

View on HN

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.

"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.

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).

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.

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.

It will result in changes in total 'energy', but that's fine if you account for this in the acceptance probability.

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.

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.

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.

Reply on news.ycombinator.com