CARPool Covariance: Fast, unbiased covariance estimation for large-scale structure observables

The covariance matrix $\boldsymbol{\Sigma}$ of non-linear clustering statistics that are measured in current and upcoming surveys is of fundamental interest for comparing cosmological theory and data and a crucial ingredient for the likelihood approximations underlying widely used parameter inference and forecasting methods. The extreme number of simulations needed to estimate $\boldsymbol{\Sigma}$ to sufficient accuracy poses a severe challenge. Approximating $\boldsymbol{\Sigma}$ using inexpensive but biased surrogates introduces model error with respect to full simulations, especially in the non-linear regime of structure growth. To address this problem we develop a matrix generalization of Convergence Acceleration by Regression and Pooling (CARPool) to combine a small number of simulations with fast surrogates and obtain low-noise estimates of $\boldsymbol{\Sigma}$ that are unbiased by construction. Our numerical examples use CARPool to combine GADGET-III $N$-body simulations with fast surrogates computed using COmoving Lagrangian Acceleration (COLA). Even at the challenging redshift $z=0.5$, we find variance reductions of at least $\mathcal{O}(10^1)$ and up to $\mathcal{O}(10^4)$ for the elements of the matter power spectrum covariance matrix on scales $8.9\times 10^{-3}<k_\mathrm{max} <1.0$ $h {\rm Mpc^{-1}}$. We demonstrate comparable performance for the covariance of the matter bispectrum, the matter correlation function and probability density function of the matter density field. We compare eigenvalues, likelihoods, and Fisher matrices computed using the CARPool covariance estimate with the standard sample covariance estimators and generally find considerable improvement except in cases where $\Sigma$ is severely ill-conditioned.