Fast optimal CMB power spectrum estimation with Hamiltonian sampling

We present a method for fast optimal estimation of the temperature angular power spectrum from observations of the cosmic microwave background. We employ a Hamiltonian Monte Carlo (HMC) sampler to obtain samples from the posterior probability distribution of all the power spectrum coefficients given a set of observations. We compare the properties of the HMC and the related Gibbs sampling approach on low-resolution simulations. The correlation lengths of the samples produced by Gibbs and HMC are generally comparable. At very high signal-to-noise (greater than or similar to 10) the Gibbs outperforms the HMC while for intermediate signal-to-noise (similar to 0.1-10) the HMC results in shorter correlation lengths than the Gibbs. At low signal-to-noise (less than or similar to 0.1) both samplers appear to perform similarly. We also demonstrate the method on high-resolution data by applying it to simulated Wilkinson Microwave Anisotropy Probe (WMAP) data. Analysis of a single map from a WMAP-resolution data set is possible in around 80 h on a high-end desktop computer. HMC uses around a factor of 3-4 fewer spherical harmonic transforms than is reported for a Gibbs sampler (with pre-conditioning) and therefore produces noticeable performance gains. HMC imposes few conditions on the distribution to be sampled and obviates the requirement of having to be able to sample efficiently from the conditional densities, thus providing us with an extremely flexible approach upon which to build.