  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.
