Matrix-free Large Scale Bayesian inference in cosmology

In this work we propose a new matrix-free implementation of the Wiener sampler which is traditionally applied to high dimensional analysis when signal covariances are unknown. Specifically, the proposed method addresses the problem of jointly inferring a high dimensional signal and its corresponding covariance matrix from a set of observations. Our method implements a Gibbs sampling adaptation of the previously presented messenger approach, permitting to cast the complex multivariate inference problem into a sequence of uni-variate random processes. In this fashion, the traditional requirement of inverting high dimensional matrices is completely eliminated from the inference process, resulting in an efficient algorithm that is trivial to implement. Using cosmic large scale structure data as a showcase, we demonstrate the capabilities of our Gibbs sampling approach by performing a joint analysis of three dimensional density fields and corresponding power-spectra from Gaussian mock catalogues. These tests clearly demonstrate the ability of the algorithm to accurately provide measurements of the three dimensional density field and its power-spectrum and corresponding uncertainty quantification. Moreover, these tests reveal excellent numerical and statistical efficiency which will generally render the proposed algorithm a valuable addition to the toolbox of large scale Bayesian inference in cosmology and astrophysics.