Bayesian Control Variates for optimal covariance estimation with pairs of simulations and surrogates

Predictions of the mean and covariance matrix of summary statistics are critical for confronting cosmological theories with observations, not least for likelihood approximations and parameter inference. The price to pay for accurate estimates is the extreme cost of running $N$-body and hydrodynamics simulations. Approximate solvers, or surrogates, greatly reduce the computational cost but can introduce significant biases, for example in the non-linear regime of cosmic structure growth. We propose "CARPool Bayes", an approach to solve the inference problem for both the means and covariances using a combination of simulations and surrogates. Our framework allows incorporating prior information for the mean and covariance. We derive closed-form solutions for Maximum A Posteriori covariance estimates that are efficient Bayesian shrinkage estimators, guarantee positive semi-definiteness, and can optionally leverage analytical covariance approximations. We discuss choices of the prior and propose a simple procedure for obtaining optimal prior hyperparameter values with a small set of test simulations. We test our method by estimating the covariances of clustering statistics of GADGET-III $N$-body simulations at redshift $z=0.5$ using surrogates from a 100-1000$\times$ faster particle-mesh code. Taking the sample covariance from 15,000 simulations as the truth, and using an empirical Bayes prior with diagonal blocks, our estimator produces nearly identical Fisher matrix contours for $\Lambda$CDM parameters using only $15$ simulations of the non-linear dark matter power spectrum. In this case the number of simulations is so small that the sample covariance would be degenerate. We show cases where even with a na\"ive prior our method still improves the estimate. Our framework is applicable to a wide range of cosmological and astrophysical problems where fast surrogates are available.