  Computing the inverse covariance matrix (or precision matrix) of large data
vectors is crucial in weak lensing (and multi-probe) analyses of the large
scale structure of the universe. Analytically computed covariances are
noise-free and hence straightforward to invert, however the model
approximations might be insufficient for the statistical precision of future
cosmological data. Estimating covariances from numerical simulations improves
on these approximations, but the sample covariance estimator is inherently
noisy, which introduces uncertainties in the error bars on cosmological
parameters and also additional scatter in their best fit values. For future
surveys, reducing both effects to an acceptable level requires an unfeasibly
large number of simulations.
  In this paper we describe a way to expand the true precision matrix around a
covariance model and show how to estimate the leading order terms of this
expansion from simulations. This is especially powerful if the covariance
matrix is the sum of two contributions, $\smash{\mathbf{C} =
\mathbf{A}+\mathbf{B}}$, where $\smash{\mathbf{A}}$ is well understood
analytically and can be turned off in simulations (e.g. shape-noise for cosmic
shear) to yield a direct estimate of $\smash{\mathbf{B}}$. We test our method
in mock experiments resembling tomographic weak lensing data vectors from the
Dark Energy Survey (DES) and the Large Synoptic Survey Telecope (LSST). For DES
we find that $400$ N-body simulations are sufficient to achive negligible
statistical uncertainties on parameter constraints. For LSST this is achieved
with $2400$ simulations. The standard covariance estimator would require
>$10^5$ simulations to reach a similar precision. We extend our analysis to a
DES multi-probe case finding a similar performance.
