Higher Order Hamiltonian Monte Carlo Sampling for Cosmological Large-Scale Structure Analysis

We investigate higher order symplectic integration strategies within Bayesian cosmic density field reconstruction methods. In particular, we study the fourth order discretisation of Hamiltonian equations of motion. This is achieved through an operator formalism, in which the original leap-frog algorithm is recursively applied in a combination of two forward time integration steps with an intermediate backward step with appropriate step-sizes. We restrict this study to the lognormal-Poisson model applied to a full volume halo catalogue in real space on a cubical mesh of 1250 h^{-1} Mpc side and 256^3 cells. Hence, we neglect selection effects, redshift space distortions, and displacements. We note that those observational and cosmic evolution effects can be accounted for in subsequent Gibbs-sampling steps within the COSMIC BIRTH algorithm. We find that going from the usual second to fourth order in the leap-frog scheme improves the convergence by a factor of ~20 in computing time, increasing the acceptance rate from 52 to 79%. Moreover, we obtain a correlation length of about 10 iterations, as opposed to ~300. This gain in computational efficiency is crucial to go towards a full Bayesian analysis of the cosmological large-scale structure for upcoming galaxy surveys.