Hamiltonian sampling

One of the problems I am interested is the power spectrum estimation of cosmological data in the sky. Hamiltonian can be very effective here in sampling the large dimensional (~10^6) posterior probability distributions. If the posterior distribution is unimodal this technique can be very effective. I have written several codes (in C,fortran and CUDA) which can be found here. CUDA version actually outperforms CPU version as long as as one can keep all the data inside the GPU memory.

Here are some papers I found very useful for understanding the fundamental concepts of Hamiltonian Monte Carlo (HMC)

1. Neal, 2011, MCMC using Hamiltonian dynamics
2. Beskos et al. 2010, The Acceptance Probability of the Hybrid Monte Carlo Method in High-Dimensional Problems
3. Girolami & Calderhead, 2011, Riemann manifold Langevin and Hamiltonian Monte Carlo methods
4. Hanson, 2001, Markov chain Monte Carlo posterior sampling with the Hamiltonian method
5. Hoffman & Gelman, 2014, The No-U-Turn Sampler: Adaptively Setting Path Lengths in Hamiltonian Monte Carlo

I started with a fortran code which can be found here, then moved to a c version called Ellipsis and finally created a CUDA version. For the CUDA version I get really big speed up. You can see the plots below.