Adaptive Bayesian Sampling with Monte Carlo EM
We present a novel technique for learning the mass matrices in samplers obtained from discretized dynamics that preserve some energy function. Existing adaptive samplers use Riemannian preconditioning techniques, where the mass matrices are functions of the parameters being sampled. This leads to significant complexities in the energy reformulations and resultant dynamics, often leading to implicit systems of equations and requiring inversion of high-dimensional matrices in the leapfrog steps. Our approach provides a simpler alternative, by using existing dynamics in the sampling step of a Monte Carlo EM framework, and learning the mass matrices in the M step with a novel online technique. We also propose a way to adaptively set the number of samples gathered in the E step, using sampling error estimates from the leapfrog dynamics. Along with a novel stochastic sampler based on Nosé-Poincaré dynamics, we use this framework with standard Hamiltonian Monte Carlo (HMC) as well as newer stochastic algorithms such as SGHMC and SGNHT, and show strong performance on synthetic and real high-dimensional sampling scenarios; we achieve sampling accuracies comparable to Riemannian samplers while being significantly faster.
This paper considers the problem of learning the mass matrix within Hamiltonian type MCMC samplers.  The authors do so by proposing the use of an online EM style algorithm with a Monte Carlo sampling step that adjusts according to the sampling error estimates.
The issue of how to sample efficiently from complex target densities is a very important one and the authors make a very nice contribution in this paper.  It is rather densely written but the topics are well explained with thorough referencing to the existing literature. The idea of making use of EM to estimate the mass matrix in these types of algorithms is very interesting, the proposed algorithms are explained well and appear to perform well in practice.  I therefore only have the following relatively minor comments:
Line 128:  “we perform generalised leapfrog steps” - it seems to me that standard leapfrog integrator is used here, rather than the generalised one?
Algorithm 1:  I find the use of \epsilon in the algorithm description slightly confusing.  Perhaps you could use an alternative notation here to make clear that these integration steps are only the proposed points for the subsequent Metropolis-Hastings step?
Line 328:  typo “we put model the columns”
Section 4.1:  Why was the observed FI used here instead of the expected FI?  I think the expected FI might give better performance.  More generally, this example seems to be a bit of a straw man, given that there isn’t any strong correlation here in the posterior density - a simplified Riemannian MALA sampler would likely outperform any of these HMC methods, which makes the comparison rather moot.  More interesting would be a banana shaped distribution, which is where Riemannian MCMC methods are really a bit more useful.
More generally, it is not clear to me how well the EM versions would work with strong correlations that change over the posterior.  Since it will be converging to a single mass matrix, then I suspect they won’t work well in strongly banana-like cases.  This isn’t explored in the paper at all, which is a shame - I think an example like this would be much more illustrative and informative than the current (trivial) Gaussian example.