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 univariate 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 data. 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.