Symplectic integration and dynamic Hamiltonian Monte Carlo

[ChrisRackauckas]

@Tamas_Papp

Stepsize adaptation usually happens with a heuristic. It can be problematic if the logdensity has highly varying curvature, as the sampler cannot adapt to the whole space well. The practical solution for this is jittering the stepsize and hoping for the best, and for truly nasty models think about reparametrizations (cf the folk theorem of statistical computing ).

What symplectic integrators are actually doing is something that’s global, as shown here:

What does "symplectic" mean in reference to numerical integrators, and does SciPy's odeint use them?

So this kind of stepsize jittering actually gets rid of the symplectic property (which is why I said adaptive timestepping for symplectic is very difficult). However, HMC is using a symplectic integrator because of its volume preservation property meaning less rejections, not because of its long-term integration properties, so it’s probably fine to do with HMC but not a true integration.

@Tamas_Papp

So the real question is: compared to leapfrog, what would be the accuracy of other algorithms traversing the same trajectory (we want to get far), and the same number of logdensity evaluations ? I don’t know enough about numerical ODE to answer this, but if someone suspects that they can do better than leapfrog on problems of varying curvature while not increasing the function evaluations drastically, that could be a big win as it would allow another layer of automation. Thinking about reparametrizations that “smooth” problems is probably the trickiest part of practical MCMC now for complex models.

That’s what methods like Highly stable symplectic methods · Issue #82 · SciML/OrdinaryDiffEq.jl · GitHub are supposed to do by increasing the stability region proportionally larger than the stepsize. I need to do some more research in this though, and I want to go to higher order.