Precision matrix expansion - efficient use of numerical simulations in estimating errors on cosmological parameters

TL;DR

This paper introduces a precision matrix expansion method that efficiently estimates errors on cosmological parameters using fewer simulations by leveraging analytical models and simulation data, improving accuracy for large-scale structure analyses.

Contribution

The paper presents a novel expansion technique for the precision matrix that reduces the number of simulations needed to accurately estimate cosmological parameter errors.

Findings

400 simulations suffice for DES to achieve negligible uncertainties

2400 simulations are enough for LSST to reach similar precision

Standard covariance estimation requires over 100,000 simulations for comparable accuracy

Abstract

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.

Tables

Table 1. Table 1: Number of tomographic bins, total number of data points and type of data vector for the different setups used to test the performance of precision matrix expansion.

setup	survey	lens bins	source bins	$N_{\mathrm{data}}$	data
Ia	DES	0	5	450	real space
Ib	DES	3	5	630	real space
II	LSST	0	10	2200	Fourier

Equations

$$p(\mathbf{\pi }|\hat{\mathbf{\xi }})\sim \exp \left(-\frac{1}{2}{\chi }^2\left[\mathbf{\pi }|\hat{\mathbf{\xi }},\mathbf{C}\right]\right)p(\mathbf{\pi })$$

$${\chi }^2\left[\mathbf{\pi }|\hat{\mathbf{\xi }},\mathbf{C}\right]={\left(\hat{\mathbf{\xi }}-\mathbf{\xi }[\mathbf{\pi }]\right)}^T{\mathbf{C}}^{-1}\left(\hat{\mathbf{\xi }}-\mathbf{\xi }[\mathbf{\pi }]\right)$$

$$\hat{\mathbf{C}}:=\frac{1}{\nu }\sum _{i=1}^{N_s}\left({\hat{\mathbf{\xi }}}_i-\bar{\mathbf{\xi }}\right){\left({\hat{\mathbf{\xi }}}_i-\bar{\mathbf{\xi }}\right)}^T,$$

$$\hat{\mathbf{\Psi }}=\frac{\nu -N_d-1}{\nu }{\hat{\mathbf{C}}}^{-1}$$

$${\hat{\mathbf{\pi }}}_{\mathrm{ML}}=\underset{\mathbf{\pi }}{\min }\left\{{\left(\hat{\mathbf{\xi }}-\mathbf{\xi }[\mathbf{\pi }]\right)}^T\hat{\mathbf{\Psi }}\left(\hat{\mathbf{\xi }}-\mathbf{\xi }[\mathbf{\pi }]\right)\right\}.$$

$$\mathbf{C}=\mathbf{A}+\mathbf{B},$$

$$\mathbf{C}=\mathbf{M}+(\mathbf{B}-{\mathbf{B}}_m),$$

$$\mathbf{C}=\left(1+\mathbf{X}\right)\mathbf{M},$$

$$\mathbf{X}:=(\mathbf{B}-{\mathbf{B}}_m){\mathbf{M}}^{-1}.$$

$$\mathbf{\Psi }$$

$${\hat{\mathbf{\Psi }}}_{1\mathrm{st}}={\mathbf{M}}^{-1}-{\mathbf{M}}^{-1}\left(\hat{\mathbf{B}}-{\mathbf{B}}_m\right){\mathbf{M}}^{-1}.$$

$${\hat{\mathbf{\Psi }}}_{2\mathrm{nd}}$$

$$\mathbf{A}={\mathbf{C}}^{nn}+{\mathbf{C}}^{sn}$$

$$\mathbf{B}={\mathbf{C}}^{ss,\mathrm{halo}}.$$

$${\mathbf{B}}_m=\alpha {\mathbf{C}}^{ss,\mathrm{Gauss}}+\beta \left({\mathbf{C}}^{ss,\mathrm{halo}}-{\mathbf{C}}^{ss,\mathrm{Gauss}}\right)$$

$$({\Omega }_m,{\sigma }_8,w_0,w_a)=$$

$$(0.3156,0.831,-1,0).$$

$$p(\mathbf{\pi }|\hat{\mathbf{\xi }})\sim {\left[1+\frac{{\left(\hat{\mathbf{\xi }}-\mathbf{\xi }[\mathbf{\pi }]\right)}^T{\hat{\mathbf{C}}}^{-1}\left(\hat{\mathbf{\xi }}-\mathbf{\xi }[\mathbf{\pi }]\right)}{N_s-1}\right]}^{-N_s/2}$$

$${\mathbf{\xi }}_{\mathrm{simple}}[\mathbf{\pi }]={\mathbf{\xi }}_{\mathrm{exact}}[{\mathbf{\pi }}_0]+\sum _{i=1}^{N_p}({\pi }_i-{\pi }_{0,i})\frac{\partial {\mathbf{\xi }}_{\mathrm{exact}}}{\partial {\pi }_i}[{\mathbf{\pi }}_0].$$

$$(b_1,b_2,b_3)=(1.35,1.50,1.65)$$

$${\hat{\mathbf{\pi }}}_{\mathrm{ML}}=\underset{\mathbf{\pi }}{\min }\left\{{\left(\hat{\mathbf{\xi }}-\mathbf{\xi }[\mathbf{\pi }]\right)}^T\mathbf{\Psi }\left(\hat{\mathbf{\xi }}-\mathbf{\xi }[\mathbf{\pi }]\right)\right\}.$$

$${\mathbf{C}}_{{\hat{\mathbf{\pi }}}_{\mathrm{ML}}}={\mathbf{F}}^{-1}\left(1+\frac{(N_d-N_p)(N_s-N_d-2)}{(N_s-N_d-1)(N_s-N_d-4)}\right),$$

$$⟨{\hat{\mathbf{C}}}^2⟩\ne {\mathbf{C}}^2.$$

$$\widehat{\left({\mathbf{C}}^2\right)}=\frac{{\nu }^2{\hat{\mathbf{C}}}^2-\nu \hat{\mathbf{C}}\mathrm{tr}\hat{\mathbf{C}}}{{\nu }^2+\nu -2},$$

$${\mathbf{M}}^{1/2}{\mathbf{M}}^{1/2}=\mathbf{M}.$$

$$\mathbf{X}={\mathbf{M}}^{-1/2}\left(\mathbf{B}-{\mathbf{B}}_m\right){\mathbf{M}}^{-1/2}$$

$$\mathbf{C}={\mathbf{M}}^{1/2}\left(1+\mathbf{X}\right){\mathbf{M}}^{1/2}$$

$$\mathbf{\Psi }$$

$$|{\lambda }_i|<1,i=1,...,N_d.$$

$$\mathbf{M}=\alpha \mathbf{C}.$$

$$\mathbf{X}$$

Full text

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.

1 Introduction

Wide area surveys such as the currently running Dark Energy Survey (DES, Flaugher, 2005) or the upcoming Large Synoptic Survey Telecope (LSST, Ivezic et al., 2008) will collect vast amounts of data about the large scale structure on the universe. In cosmological analyses this data can e.g. be compressed into measurements of 2-point correlation functions of galaxy clustering or cosmic shear. In a redshift-tomographic analysis this will easily accumulate to data vectors with several hundreds of data points. Testing cosmological models from a measurement of such a large data vector requires precise knowledge of the inverse covariance matrix of the noise in this data vector. There has been extensive research on the impact of errors associated with covariance estimation on the constraints derived on cosmological parameters. Hartlap et al. (2007) discussed the fact that the inverse of an unbiased covariance estimator is not an unbiased estimator for the inverse covariance matrix (the precision matrix). They also described a way to correct for this when assuming that the covariance estimate follows a Wishart distribution (see also Kaufman 1967 and Anderson 2003). The noise properties of this corrected precision matrix estimator and its impact on the constraints derived on cosmological parameters was e.g. investigated by Taylor et al. (2013); Dodelson & Schneider (2013); Taylor & Joachimi (2014).

Sellentin & Heavens (2016a, hereafter SH16a) have presented a different approach: given a covariance estimate they marginalize over the posterior distribution of the true precision matrix to compute the likelihood in parameter space. Assuming that the covariance estimate follows a Wishart distribution they have derived a simple, closed-form expression for the resulting likelihood function. In Sellentin & Heavens (2016b) they have extended these results to derive the information loss in parameter space due to noisy covariance estimates. A fully non-Gaussian treatment of the effects discussed in Dodelson & Schneider (2013, hereafter DS13) is however still missing.

Prior knowledge on the sparsity of the covariance matrix and the precision matrix was used by Paz & Sánchez (2015) and Padmanabhan et al. (2015) to improve estimates of the precision matrix from few simulations. Pope & Szapudi (2008) investigated shrinkage estimators of the covariance, i.e. a mixing of estimated and modelled covariance matrices. This however raises the task of finding an equivalent to the Kaufman-Hartlap correction for such a mixture of estimated and analytic matrices. More recently, Joachimi (2017) describes a non-linear extension of that estimator which combines covariance estimates from two sets of independent data vector realisations and hence does not require a covariance model.

In this paper we describe a way to expand the true precision matrix around a covariance model as a power series in the deviation between model and true covariance. Assuming a Wishart realisation for the true covariance and using the results on invariant moments of the Wishart distribution by Letac & Massam (2004) we derive an unbiased estimator for the up to second order expansion of the true precision matrix. This becomes especially powerful if parts of the covariance matrix that are well understood analytically can be turned off in simulations in order to yield a direct estimate of the remaining covariance parts. In Sect. 3 we recap the main problems of estimating parameter constraints from noisy covariance estimates and present our method of "Precision Matrix Expansion" (PME). In Sect. 4 we perform numerical experiments that mimic data from the Dark Energy Survey (DES) and the Large Synoptic Survey Telecope (LSST) likelihood analyses to test the performance of our idea. Sect. 5 concludes with a discussion of our results.

2 Parameter constraints from noisy covariance estimates

We begin by outlining the main task of this paper. Let $\smash{\bm{\hat{\xi}}}$ be a vector of $N_{d}$ data points measured from observational data and let $\bm{\xi}[\bm{\pi}]$ be a model for this data vector that depends on a vector of $N_{p}$ parameters $\bm{\pi}$. If $\mathbf{C}$ is the covariance matrix of $\smash{\bm{\hat{\xi}}}$ then a standard way to constrain the parameters $\bm{\pi}$ is to assign a posterior distribution $p(\bm{\pi}|\bm{\hat{\xi}})$ to them as

[TABLE]

with

[TABLE]

and $p(\bm{\pi})$ being a prior density incorporating apriori knowledge or assumptions on $\bm{\pi}$. These expressions in fact ignore that $\mathbf{C}$ also can be depedent on $\bm{\pi}$. We will do this throughout this paper and refer the reader to Eifler et al. (2009) who investigated the impact of cosmology dependent covariance matrices on cosmic shear likelihood analyses. Another assumption that goes into Eq. 1 is that the measured data vector $\smash{\bm{\hat{\xi}}}$ is drawn from a multi-variate Gaussian distribution. In wide area surveys this is justified in the limit where one can consider the survey to consist of many independent sub-regions, such that the measurements in those regions add up to a Gaussian data vector by means of the central limit theorem.

If the covariance matrix $\mathbf{C}$ is not exactly known, it can e.g. be estimated from N-body simulations. If $\smash{\bm{\hat{\xi}}_{i}}$, $i=1...N_{s}$, are a number of independent measurements of $\bm{\xi}$ in simulations then an unbiased estimate of $\mathbf{C}$ is given by

[TABLE]

where $\nu=N_{s}-1$ and $\bm{\bar{\xi}}$ is the sample mean of the $\smash{\bm{\hat{\xi}}_{i}}$. We will assume $\smash{\mathbf{\hat{C}}}$ to have a Wishart distribution with $\nu$ degrees of freedom which follows from our assumption that $\smash{\bm{\hat{\xi}}}$ and the $\smash{\bm{\hat{\xi}}_{i}}$ are Gaussian distributed (cf. Taylor et al. 2013).

To compute the likelihood in Eq. 1 we need to know the precision matrix, i.e. is the inverse covariance matrix $\smash{\mathbf{\Psi}=\mathbf{C}^{-1}}$. According to Kaufman (1967, see also Hartlap et al. 2007; Taylor et al. 2013) an unbiased estimator for $\smash{\mathbf{\Psi}}$ can be constructed from $\smash{\mathbf{\hat{C}}}$ as

[TABLE]

and we will call the factor of $(\nu-N_{d}-1)/\nu$ the Kaufman-Hartlap-correction.

Given a measurement $\smash{\bm{\hat{\xi}}}$ of the data vector one can derive the posterior density of the model parameters $p(\bm{\pi}|\bm{\hat{\xi}})$ by means of equations 1 and 2. A noisy precision matrix estimate influences this inference in two ways:

• •

it adds noise to the width of likelihood contours derived from inserting the precision matrix estimate into the figure of merit $\chi^{2}$ (Eq. 2).

• •

it adds noise to the location of likelihood contours. E.g. the maximum likelihood estimator for the parameters $\bm{\pi}$ would be

[TABLE]

When using a noisy precision matrix the uncertainties of $\bm{\hat{\pi}}_{\mathrm{ML}}$ have contributions from both the noise in $\smash{\bm{\hat{\xi}}}$ and the noise in $\smash{\mathbf{\hat{\Psi}}}$.

The astro-statistics literature has so far focused on the first effect, i.e. on the uncertainties on contour width due to noise in the estimate $\smash{\mathbf{\hat{\Psi}}}$ (Taylor et al., 2013; Taylor & Joachimi, 2014; Sellentin & Heavens, 2016a, b). Sellentin & Heavens (2016b) provide the most complete demonstration that $\smash{\mathbf{\hat{\Psi}}}$ yields a good estimate of the width of the posterior contours as long as $N_{s}-N_{d}\gg N_{p}$.

The more critical effect however is the additional noise of $\bm{\hat{\pi}}_{\mathrm{ML}}$. DS13 (also see appendix A) showed that the uncertainty on the position of likelihood contours from noise in $\smash{\mathbf{\hat{\Psi}}}$ is only negligible if $N_{s}-N_{d}\gg N_{d}-N_{p}$ which is a much more demanding criterion for current cosmological data vectors. We demonstrate this in the left-hand panel of Fig. 1, where we show 100 randomly drawn realisations of a DES-like weak lensing data vector with $N_{d}=450$ and a halo model covariance matrix (see Sec. 4 for further details). For each of the 100 data vectors we have also generated Wishart realisations of our covariance matrix corresponding to an estimate from $N_{s}=650$ simulations. Using either the true covariance or the estimated one, we then determine the best fitting parameters $\Omega_{m}$ and $\sigma_{8}$ (after marginalizing over equation-of-state parameters of dark energy, $w_{0}$ and $w_{a}$). The best-fits obtained from a noisy covariance (green points) clearly display a much larger scatter than those obtained from the true covariance (red points). Also shown are the best fits obtained by precision matrix expansion (PME, blue points) which we are going to introduce in the next section. Here we assumed that only $N_{s}=200$ simulations are available to estimate the PME, which gives best fit values that are significantly closer to the ones obtained when knowing the true covariance matrix.

When reconstructing $p(\bm{\pi}|\bm{\hat{\xi}})$ (e.g. from a Monte-Carlo-Markov-Chain) this can lead to significant offsets between likelihood contours inferred from the true covariance matrix and likelihood contours inferred from a covariance estimate – even if the overall width of the likelihood contours is captured well by the covariance estimate. We demonstrate this in the right-hand panel of figure 1. DS13 have derived a factor (see appendix A) by which parameter contours obtained from a Wishart realisation of the covariance should be expanded in order to account for this additional scatter. However, their derivation relies on the assumption of a Gaussian parameter likelihood and is only applicable to the extent that a Fisher analysis is accurate. The current state of the art for dealing with noisy covariance estimates is hence a combination of SH16a and DS13: expanding the contours derived from the SH16a likelihood by the DS13 factor. We implement this idea for the cyan contours in Fig. 1 and show that this brings the contours derived from a standard covariance estimate into consistency with those derived from the true covariance.

Downsides of this approach are a large increase of the uncertainties on cosmological parameters and the fact that one still needs at least as many realisations as data points in the data vector to even derive a precision matrix estimate. We now want to introduce an alternative method to estimate the precision matrix which is able to drastically decrease the offset of contours seen for the standard precision matrix estimator.

3 Precision matrix expansion

Let us split the covariance matrix $\mathbf{C}$ into two contributions

[TABLE]

where for matrix $\mathbf{A}$ we have an accurate model (e.g. the shape-noise contributions to the covariance of cosmic shear correlation functions) and for $\mathbf{B}$ we have a model $\mathbf{B}_{m}$ which we know to be imperfect. We want to include this prior knowledge of the covariance matrix when estimating the precision matrix. Starting from

[TABLE]

where $\mathbf{M}=\mathbf{A}+\mathbf{B}_{m}$ is our model for the complete covariance matrix, we rewrite

[TABLE]

where $\mathbb{1}$ is the identity matrix and we have defined

[TABLE]

The precision matrix $\mathbf{\Psi}=\mathbf{C}^{-1}$ can then be expressed as the following power series in $\mathbf{X}$:

[TABLE]

We will call this series the precision matrix expansion (PME). In appendix C we show that it converges under a wide range of conditions. There we also demonstrate that the series yields at each order a symmetric approximation of $\mathbf{\Psi}$ and that at second order it is always positive definite (at each order if the series converges).

3.1 Estimating the expansion of $\mathbf{\Psi}$

Suppose we have an estimate $\smash{\mathbf{\hat{B}}}$ of the matrix $\smash{\mathbf{B}}$ from a number of N-body simulations. This especially assumes that all covariance contributions included in $\smash{\mathbf{A}}$ can be turned off in the simulations (i.e. for cosmic shear covariances $\smash{\mathbf{A}}$ could consist of shape-noise contributions which can be set to zero in simulations). We want to use $\smash{\mathbf{\hat{B}}}$ to construct unbiased estimators for the first order and second order term of the series in Eq. 10.

Our assumptions state that $\smash{\mathbf{\hat{B}}}$ is drawn from a Wishart distribution with expectation value $\mathbf{B}$. In this case also $\smash{\mathbf{M}^{-1}\mathbf{\hat{B}}\mathbf{M}^{-1}}$ is Wishart distributed but with the expectation value $\smash{\mathbf{M}^{-1}\mathbf{B}\mathbf{M}^{-1}}$. Hence an unbiased estimator for the first order PME is given by

[TABLE]

Note that this does not involve the inversion of an estimated matrix. According to Taylor et al. (2013) the standard deviation of diagonal elements of an inverse-Wishart distributed matrix is proportional to $1/\sqrt{N_{s}-N_{d}-4}$ while for Wishart distributed matrices it is only proportional to $1/\sqrt{N_{s}-1}$. Hence, avoiding the occurence of an inverted matrix estimate greatly reduces the estimation noise.

The second order term involves squares of Wishart matrices. Using the results of Letac & Massam (2004) on invariant moments of the Wishart distribution (cf. appendix B) it is still possible to construct an unbiased estimator for the second order PME as

[TABLE]

The estimator in Eq. 12 is the key result of our paper. It has two advantages over the Anderson-Hartlap corrected standard estimator. First, it only requires matrix multiplications. As a consequence, it can even be used if $N_{s}\leq N_{d}$. Second, it only needs an estimate of $\mathbf{B}$ instead of the whole covariance $\mathbf{C}$, i.e. it allows to incorporate apriori knowledge on the covariance in the form of $\mathbf{M}$ (and $\mathbf{A}$).

In the next section we demonstrate that this significantly eases the requirement of $N_{s}-N_{d}\gg N_{d}-N_{p}$. Hence, in a likelihood analysis the noise in $\smash{\mathbf{\hat{\Psi}}_{2\mathrm{nd}}}$ becomes negligible for a much smaller number of N-body simulations than required by the standard precision matrix estimator. In appendix C we also show that the bias in parameter constraints which arises from cutting the power series in Eq. 10 after a finite number of terms is negligible even for very strong deviations of our covariance model $\mathbf{M}$ from the N-body covariance $\mathbf{C}$.

4 Examples: parameter errors for LSST weak lensing and DES weak lensing and multi-probe analyses

We investigate the performance of our method in the context of ongoing and future surveys using DES and LSST as specific examples. These surveys differ in terms of survey area, galaxy number density, and redshift distribution and have different demands on the precision matrix. For DES we consider summary statistics in real space, i.e. auto- and cross-correlation functions of galaxy shear and position, for LSST we consider the corresponding Fourier quantities of a shear-shear only data vector. A summary of the scenarios considered is given in Table 1 and a more detailed description of the considered data vectors is given in appendix D.

In order to test the performance of PME we set up mock experiments where we assume the true covariance matrix of each survey to be the analytic halo-model covariance described in Krause & Eifler (2016). This model divides the covariance into three contributions: a noise-only part that consists of shape- and shot-noise contributions, $\mathbf{C}^{nn}$, a contribution from the cosmic variance of the signal, $\mathbf{C}^{ss,\mathrm{halo}}$, and a mixed term including noise and signal contributions, $\mathbf{C}^{sn}$. For shear-shear only covariances we set

[TABLE]

and

[TABLE]

The shape-noise contributions to the covariance can be modelled reliably since the ellipticity dispersion can be measured from the data itself and since the mixed term $\mathbf{C}^{sn}$ involves only the modeling of two-point statistics of the shear field. The $\mathbf{B}$ term comprises the more complex 4-point statistics of the shear field, which can be estimated from simulations by turning off shape-noise. This is more complicated for galaxy clustering where shot-noise is included in the covariance matrix (cf. 4.1.2).

In order to simulate a situation where our covariance model $\mathbf{M}=\mathbf{A}+\mathbf{B}_{m}$ deviates from the true covariance we degrade it as

[TABLE]

where $\smash{\mathbf{C}^{ss,\mathrm{Gauss}}}$ contains only the parts of the cosmic variance that are also present in a Gaussian covariance model. Hence, we allow the Gaussian and non-Gaussian cosmic variance parts to be over- or underestimated by a constant multiplicative factor. If not stated differently in this section we will use $\alpha=1.0$ and $\beta=0.5$. In appendix C we explore a wider range of rescalings and also consider more complex deformations of our fiducial covariance to show that the PME remains robust under more complicated deviations of $\mathbf{M}$ from the true covariance matrix. All simulated likelihood analyses in this paper are computed using the CosmoLike cosmology package (Krause & Eifler, 2016; Eifler et al., 2014).

4.1 Performance for DES weak lensing data vector

We now carry out mock likelihood analyses for DES and LSST weak lensing data vectors, varying the parameters $\Omega_{m}$, $\sigma_{8}$, $w_{0}$ and $w_{a}$. Our fiducial values for these parameters are

[TABLE]

We start by drawing random Gaussian realisations of our fiducial data vectors according to a covariance given by the halo model. For each realisation we also draw new Wishart realisations $\mathbf{\hat{B}}$ and $\mathbf{\hat{C}}$ of cosmic variance and total covariance to compute the PME estimate $\smash{\mathbf{\hat{\Psi}}_{\mathrm{2nd}}}$ and the standard estimator $\smash{\mathbf{\hat{\Psi}}}$. In practice, this is done by drawing additional realisations $\smash{\bm{\hat{\xi}}_{i}}$ , $i=1\ \dots\ N_{s}$, of our fiducial data vector from a multivariate Gaussian distribution whose covariance is $\mathbf{B}$ respectively $\mathbf{C}$. These realisations represent measurements from N-body simulations and inserting them into Eq. 3 generates the desired Wishart realisations $\mathbf{\hat{B}}$ and $\mathbf{\hat{C}}$ of the two matrices.

Using CosmoLike we then run likelihood chains to infer a posterior distribution for our parameters using $\smash{\mathbf{\hat{\Psi}}_{\mathrm{2nd}}}$, $\smash{\mathbf{\hat{\Psi}}}$ and the true precision matrix $\smash{\mathbf{C}^{-1}}$. When computing the likelihood from $\smash{\mathbf{C}^{-1}}$ and $\smash{\mathbf{\hat{\Psi}}_{\mathrm{2nd}}}$ we simply use standard ansatz given in Eq. 1. When deriving contours from the Wishart realisation $\smash{\mathbf{\hat{C}}}$ we furthermore compute the parameter likelihood as

[TABLE]

which SH16a have shown to be a more accurate than using the Kaufman-Hartlap correction and the standard Gaussian likelihood. We however found only small differences to using the standard likelihood ansatz, which is due to the fact that in all cases considered in this paper $N_{s}-N_{d}\gg N_{p}$.

In Fig. 2 and 3 we show see the resulting 1$\sigma$ and 2$\sigma$ contours in the $\Omega_{m}$-$\sigma_{8}$ plane (after marginalizing over the other parameters) for 3 different random draws of data vector and Wishart matrices. For each realisation of the DES data vector we assumed that $N_{s}=200$ simulations are available to estimate the PME and $N_{s}=N_{d}+200=650$ simulations for the standard estimator. For each realisation of the LSST data vector we assumed $N_{s}=400$ simulations for the PME and $N_{s}=N_{d}+400=2600$ simulations for the standard estimator.

Even though in each case we assumed many more simulations for the standard estimator than for the PME, the PME is significantly better in reconstructing the contours from the true precision matrix. In particular we find that deviations from the true contours are much smaller than the corresponding 1$\sigma$ and 2$\sigma$ uncertainties of the parameters.

Next we generalize the findings in Figs. 2 and 3. We generate $1000$ Wishart realisations of the matrices $\smash{\mathbf{\hat{C}}}$ and $\smash{\mathbf{\hat{B}}}$ for different assumptions on the number of available N-body simulations $N_{s}$. For each of the $1000$ sets of matrices we also generate $10$ realizations $\smash{\bm{\hat{\xi}}}$ of our fiducial data vector (i.e. overall $10\,000$ different realizations $\smash{\bm{\hat{\xi}}}$). Hence for each type of precision matrix estimate we perform overall $10\,000$ likelihood analyses. In each analysis we determine the best fit parameters $\bm{\hat{\pi}}_{\mathrm{ML}}$ and check whether our fiducial cosmology is outside the $68.3\%$ confidence contour around these parameters. In order to make this computationally feasible, we are now linearly approximating the calculations of CosmoLike around our fiducial cosmology $\bm{\pi}_{0}$, i.e. we use

[TABLE]

This allows us to analytically determine the maximum likelihood parameters and the $68.3\%$ confidence contours in each likelihood analysis. It is also the situation where a Fisher-matrix formalism and hence the derivations of DS13 hold exactly.

We define $\smash{F_{>1\sigma}}$ as the fraction of times that our fiducial cosmology is outside of the $68.3\%$ confidence contour around the best fit parameters and we use it as a metric for comparing the different precision matrix estimators. In Fig. 4 we show this fraction for all different types of precision matrices introduced before. The solid, dashed, and dotted lines show the fractions achieved when using the noise-less matrices $\smash{\mathbf{C}^{-1}}$, $\smash{\mathbf{M}^{-1}}$ and $\smash{\mathbf{\Psi}_{\mathrm{2nd}}}$. Especially, the noise-less matrix $\smash{\mathbf{\Psi}_{\mathrm{2nd}}}$ would be the PME-estimator in the limit of infinitely many simulations and $\smash{\mathbf{C}^{-1}}$ would be the standard estimator in the same limit. The red and blue dots show the fraction achieved when using the noisy precision matrix estimates $\smash{\mathbf{\hat{\Psi}}}$ and $\smash{\mathbf{\hat{\Psi}}_{\mathrm{2nd}}}$.

As expected, $\smash{F_{>1\sigma}}$ is very close to $32\%$ when using the true covariance $\mathbf{C}$ in the likelihood analyses. For the deformed halo model covariance $\mathbf{M}$ we assumed the two cases $\alpha=0.7$, $\beta=0.5$ (left panel) and $\alpha=1.0$ and $\beta=0.5$ (right panel). For $\alpha=0.7$ and $\beta=0.5$ our fiducial cosmology is regarded as outside the $68.3\%$ contour in more than $40\%$ of the cases. For both choices of $\mathbf{M}$ the noise-free PME significantly corrects that fraction towards the optimal value of $\sim 32\%$. Especially promising is that the PME estimate performs very similar to the noise-free PME. If 200 simulations are available to estimate the PME, it essentially converges to its best possible performance. And even if only $100$ simulations are available to estimate the PME, its value of $\smash{F_{>1\sigma}}$ comes closer to $32\%$ than when using $\mathbf{M}$ to derive the contours.

When inferring the likelihood from the standard precision matrix estimator $\smash{F_{>1\sigma}}$ is greater than $50\%$ even if we allow $N_{s}=N_{d}+800$ simulations for the covariance estimation, which corresponds to $1250$ simulations. This is due to the additional variance of $\bm{\hat{\pi}}_{\mathrm{ML}}$ caused by the noise of the precision matrix (cf. Eq. 5). Using the results of DS13 we can derive predictions for this effect (cf. appendix A). As can be seen from the red dashed lines in Fig. 4 these predictions agree well with what we find in our simulated likelihood analyses. Extrapolating the results of DS13 to higher values of $N_{s}$ we can also estimate, how many simulations would be required for the standard precision matrix estimator in order to achieve the same value of $\smash{F_{>1\sigma}}$ as the second order PME. For the left panel of Fig. 4 we find that it would take $\sim 8000$ simulations for the standard estimator to get as close to $\smash{F_{>1\sigma}=32\%}$ as the PME with only $200$ simulations. This statement however depends on the model covariance $\mathbf{M}$ since it determines how well the PME has converged after its second order.

An $\mathbf{M}$-independent way of comparing standard estimator and PME estimator is to see how many simulations it takes each to have $F_{>1\sigma}$ within $1\%$ of their best possible performance. It would take the standard estimator $\sim 24\,000$ simulations to be within $1\%$ of $F_{>1\sigma}=32\%$. The PME estimator is well within $1\%$ of its best possible performance for only 200 simulations.

Note that with the results of DS13 one can in principle correct a likelihood analysis for the additional variance caused by the standard precision matrix estimator. This would result in a decreased constraining power of the analysis and it would hence be the main benefit of the PME to prevent this loss.

4.1.1 Larger covariance matrices: LSST weak lensing data vector

We repeat the above analysis for the LSST-like weak lensing data vector. Fig. 5 shows the fractions $F_{>1\sigma}$ obtained from PME and standard precision matrix estimator. The PME estimator now requires $\sim 2400$ simulations to be less than $1\%$ away from its best possible performance. As before, this statement does not include any additional biases between PME and true precision matrix that might arise from the biased model matrix $\mathbf{M}$ used to carry out the matrix expansion. The standard precision matrix estimator would need $N_{s}>115\,000$ simulations to be less than $1\%$ away from its best possible performance.

4.1.2 Defining $\mathbf{A}$ and $\mathbf{B}$ for multi-probe covariances

We now repeat the analysis of Fig. 4 for a DES-like multi-probe data vector. This vector includes contributions from galaxy clustering and galaxy-galaxy lensing, which introduces shot-noise terms to the covariance. These shot-noise contributions are in principle well understood theoretically and include, similar to the cosmic shear case, at most two-point statistics of the cosmic density field. Hence, one could absorb them into the matrix $\mathbf{A}$ (cf. Eq. 6) and use N-body simulations only for the remaining part of the covariance - i.e. to define $\mathbf{B}$ as only the cosmic variance. This is however difficult since most N-body simulations provide only simulated galaxy catalogues that are affected by shot noise themselves, which makes it impossible to independently estimate the cosmic variance. If however all shot-noise contributions are included in $\mathbf{B}$ when defining and estimating the PME, then the estimator $\smash{\mathbf{\hat{\Psi}}_{2\mathrm{nd}}}$ will have a higher variance in many of its elements. Hence, the additional scatter of best fit parameters due to a noisy precision matrix might not be negligible anymore.

In Fig. 6 we compare the fractions $F_{>1\sigma}$ obtained from different estimates of the precision matrix in our simulated likelihood analyses - this time for the DES multi-probe data vector. In each likelihood analysis we now vary $7$ parameters, since for each lens bin we include a galaxy bias parameter in our model. The fiducial bias values are

[TABLE]

in order of increasing redshift. Fig. 6 shows the results obtained for each of the mentioned options of defining $\mathbf{B}$. It is clear that the noisy PME approaches its best possible performance already for a smaller number of simulations if the cosmic variance can be estimated directly. In practice this would however require density maps in thin redshift slices for each simulation in order to measure the correlation functions of the projected density fields without shot-noise.

Assuming one can directly measure the cosmic variance from simulations we again want to asses how many simulations are required for the standard precision matrix estimator and the PME estimator to be within $1\%$ of their best possible performance. Extrapolating the results of DS13 we find that it would take the standard estimator $\sim 44\,000$ simulations to be within within $1\%$ of $F_{>1\sigma}=32\%$. The PME estimator is within $1\%$ of its best possible performance for 1600 simulations. For $N_{s}=2000$ the performance of the PME becomes almost solely restricted by the deviation between $\mathbf{M}$ and $\mathbf{C}$ in our mock experiment. However, below $N_{s}=1600$ there seems to be significant additional scatter of the best fitting parameters due to the noise of the PME estimate. We demonstrate this in Fig. 7 for $N_{s}=400$. Regardless of how $\mathbf{B}$ is defined, we can nevertheless conclude that also for multi-probe covariances the PME poses a vast improvement over the standard precision matrix estimator.

5 Conclusions

It was the starting point of our analysis to find a method for using a priori knowlegde about the covariance matrix when estimating the precision matrix from simulations. This requires finding an equivalent of the Kaufman-Hartlap correction when only parts of the covariance are estimated. Using the results of Letac & Massam (2004) we partly solved this task by calculating an expansion of the precision matrix and showing how the leading terms of this expansion can be estimated from simulations. Our method enables the use of preexisting knowledge on the covariance structure to improve the convergence of the PME and to reduce the noise in its estimation. It also has the advantage that the relative uncertainties of the elements of the PME estimate scale with the number of available simulations $N_{s}$ as $\sim 1/\sqrt{N_{s}-1}$, which is typically much smaller than the uncertainties of the standard precision matrix estimator. The latter also depend on the number of data points $N_{d}$ and scale as $\sim 1/\sqrt{N_{s}-N_{d}-4}$.

We demonstrated that the PME converges even for drastic deviations between the model covariance and the N-body covariance and we also showed that it provides a much less noisy estimate of the parameter likelihood compared to estimating the precision matrix in the standard way. For a DES weak lensing data vector $N_{s}\gtrsim 8000$ simulations would be required for the standard estimator to reconstruct the likelihood similarly well as the PME with only $N_{s}=200$ - even if the model covariance heavily underestimates Gaussian and non-Gaussian covariance parts. If we assume more realistic deviations between model and N-body covariance, up to $24\,000$ simulations would be needed for the standard estimator to reconstruct the 1$\sigma$ quantile of the parameter distribution at the same precision as the PME with only 200 simulations. For an LSST-like weak lensing data vector with $N_{d}=2200$ we found that up to $115\,000$ simulations would be required for the standard estimator to reconstruct the 1$\sigma$ quantile as well as the PME with only 2400 simulations. It should however be stressed that these statements depend on the quality of the model covariance $\mathbf{M}$ that was used to compute the PME.

Additional complications arise when galaxy clustering correlation functions are included in the data vector. A performance similar to the weak lensing case can still be achieved if one manages to estimate the cosmic variance of the correlation functions directly, i.e. without shot-noise. For this case, we find that a DES-like multi-probe data vector requires up to $44\,000$ simulations for the standard precision matrix estimator to reconstruct the 1$\sigma$ quantile of the parameter distribution as well as the PME with 1600 simulations.

One aspect that should be addressed in future work, is to find a priori criteria for the convergence of the PME. In appendix C we demonstrate that it converges for very strong deformations of the halo-model covariance, but one can not be certain whether and how fast it will converge for all possible data vectors and covariance models. As we show in appendix C, situations where the PME does not converge can at least be identified a posteriori by a comparison of the first order and second order expansion. A strong oscillation of likelihood contours derived from the first and second order PME indicate a significant deviation of model and N-body covariance. This way, the PME provides a clear criterion for testing covariance models with simulations – even when the number of available simulations is small.

Acknowledgments

This work was supported by SFB-Transregio 33 ‘The Dark Universe’ by the Deutsche Forschungsgemeinschaft (DFG). We also acknowledge the support by the DFG Cluster of Excellence "Origin and Structure of the Universe". The simulations have been carried out on the computing facilities of the Computational Center for Particle and Astrophysics (C2PAP). This research is partially supported by NASA ROSES ATP 16-ATP16-0084 grant. Part of the research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration.

This paper has gone through internal review by the DES collaboration. We want to thank Eric Baxter, Gary Bernstein, Scott Dodelson, Franz Elsner, Juan Garcia-Bellido, John Peacock, Stella Seitz and the DES multi-probe pipelines group for very helpful comments and discussions on this project.

Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey. The DES data management system is supported by the National Science Foundation under Grant Number AST-1138766.

The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Enérgeticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l’Espai (IEEC/CSIC), the Institut de Física d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University.

The DES participants from Spanish institutions are partially supported by MINECO under grants AYA2012-39559, ESP2013-48274, FPA2013-47986, and Centro de Excelencia Severo Ochoa SEV-2012-0234. Research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) including ERC grant agreements 240672, 291329, and 306478.

Appendix A Influence of noisy covariance estimates on the scatter of best fitting cosmological parameters

Using a noisy precision matrix estimate to determine the best fitting cosmological parameters

[TABLE]

leads to an additional scatter in these parameters. Especially, this additional noise is not accounted for by the width of contours generated from the precision matrix estimate. This effect has e.g. been described by Dodelson & Schneider (2013) who also derived a prediction for the additional noise assuming a Gaussian parameter likelihood. They find that the actual parameter covariance when using an inverse-Wishart realisation of the precision matrix is given by

[TABLE]

where $N_{p}$ is the number of considered parameters and $\mathbf{F}$ is the Fisher matrix computed from the true precision matrix. Hence, in the case of a Gaussian parameter likelihood, best fit parameters $\bm{\hat{\pi}}_{\mathrm{ML}}$ that are computed from a Wishart realisation of the covariance have also a Gaussian distribution but with a rescaled parameter covariance.

Appendix B Unbiased estimator of the square of a Wishart matrix

Let $\smash{\mathbf{\hat{C}}}$ be distributed according to a Wishart distribution with $\nu$ degrees of freedom and expectation value $\mathbf{C}$. Then

[TABLE]

However, using the results of Letac & Massam (2004) it is possible to devise an unbiased estimator of $\mathbf{C}^{2}$. It is given by

[TABLE]

where $\smash{\mathrm{tr}\mathbf{\hat{C}}}$ denotes the trace of $\smash{\mathbf{\hat{C}}}$. Using this formula, it is straight foreward to derive the estimator of the second order PME given in Eq. 12.

Appendix C General properties and convergence of the power series

C.1 General properties

In order to derive some general properties of the PME series, let us slightly change the notation of Sec. 3. First, let $\mathbf{M}^{1/2}$ be the unique symmetric and positive definite matrix such that

[TABLE]

This matrix exists as long as our covariance model $\mathbf{M}$ is positive definite. Let us then re-define

[TABLE]

where $\mathbf{M}^{-1/2}$ is the inverse of $\mathbf{M}^{1/2}$, and $\mathbf{B}$ and $\mathbf{B}_{m}$ are the same as in Sec. 3. The complete covariance can then be written as

[TABLE]

and the precision matrix expansion now reads

[TABLE]

Since both $\mathbf{M}^{-1/2}$ and $\mathbf{X}$ are symmetric matrices, it is immediately clear that this gives a symmetric approximation of $\mathbf{\Psi}$ at each order of the power series. The series converges if and only if all eigenvalues of $\mathbf{X}$ fulfill

[TABLE]

In each eigendimension of $\mathbf{X}$ the series $\smash{\left(\mathbb{1}-\mathbf{X}+\mathbf{X}^{2}+\mathcal{O}\left[\mathbf{X}^{3}\right]\right)}$ is simply the geometric series. For $|\lambda_{i}|<1$ the value of this series is $>0$ at each finite order. At second order, the value of this series is $>0$ regardless of the values of $\lambda_{i}$. Hence, the second order PME is always positive definite.

C.2 Special cases

C.2.1 Rescaling of the covariance

Let us investigate the convergence properties of the power series in Eq. 10 in a couple of special cases. We start by assuming that our model for the covariance matrix, $\mathbf{M}$, under- or overestimates the true covariance matrix by a constant factor $\alpha$, i.e.

[TABLE]

In this case we have

[TABLE]

Hence, all eigenvalues of $\mathbf{X}$ are given by $\lambda=\frac{1-\alpha}{\alpha}$. This has absolute value smaller than $1$ for all $\alpha>0.5$. This especially means that the series used to define the PME converges even if the model covariance overestimates the true covariance by an arbitrarily high overall factor. Since we cut Eq. 10 after the second order we must however look at how well the series is converged after that order. The relative error on each element of the precision matrix is given by

[TABLE]

This is $<10\%$ for $\alpha\in[0.69,1.86]$ and $<1\%$ for $\alpha\in[0.83,1.27]$.

C.2.2 Partial rescaling of the covariance

Now let us assume that $\mathbf{C}$ falls into two contributions $\mathbf{A}$ and $\mathbf{B}$ and that only $\mathbf{B}$ is mischaracterized by a constant factor in our model,

[TABLE]

Let us furthermore assume that $\mathbf{B}$ has a dominant eigenvalue $\lambda$ and that $\mathbf{v}$ is an eigenvector to it. If

[TABLE]

then the matrix $\mathbf{C}$ and – for values of $\alpha$ that are not too small – also the matrix $\mathbf{M}$ will have an eigendimension close the that of $\mathbf{B}$ with eigenvalues $\lambda_{C}\approx\lambda$ and $\lambda_{M}\approx\alpha\lambda$. As a consequence, the matrix $\mathbf{X}$ will have an eigendimension with eigenvalue close to $\lambda_{X}\approx\frac{1-\alpha}{\alpha}$ which allows the same conclusion in C.2.1.

In section 4 we considered a deformation of the halo model covariance of the form $\mathbf{M}=\mathbf{A}+\mathbf{B}_{m}$ with

[TABLE]

This is similar to the situation described above. To illustrate how the rescaling factors $\alpha$ and $\beta$ impact the convergence of the PME we can e.g. compare the Fisher contours derived from $\mathbf{C}^{-1}$, $\mathbf{M}^{-1}$, $\mathbf{\Psi}^{1st}$ and $\mathbf{\Psi}^{2nd}$. In Fig. 8 we show the $1\sigma$ and $2\sigma$ Fisher-contours for the parameter pair $\Omega_{m}$-$\sigma_{8}$ derived for the DES multi-probe data vector using different values of $\alpha$ and $\beta$. The figure shows that the PME manages to correct the bias between contours derived from $\mathbf{C}^{-1}$ and contours derived from $\mathbf{M}^{-1}$ even for rather drastic choices of the rescaling factors. Especially for $\alpha,\ \beta>1.0$ the convergence is very robust. As predicted by our considerations above, it however breaks down for $\alpha,\ \beta<0.5$ where one can see strong oscillations between $\mathbf{\Psi}_{1\mathrm{st}}$ and $\mathbf{\Psi}_{2\mathrm{nd}}$. The convergence of the contours in Fig. 8 is very similar when other parameter combinations are considered or when the contours are derived for the other data vectors considered in this paper.

C.2.3 Log-normal motivated approximation to the halo-model covariance

Motivated by the work of Hilbert et al. (2011) on approximating the shear-shear covariance matrix with a log-normal approach (cf. their equation 26) we approximate the non-Gaussian parts of the covariance of shear correlation functions as

[TABLE]

where $\theta_{i}$ labels the different angular bins, $A$ and $B$ label the different auto- and cross-correlation functions and $R_{AB}$ is just a constant factor (depending only on the pair $A,B$ and not on whether $\xi_{+}$ or $\xi_{-}$ are involved). We fix the values of $R_{AB}$ by demanding that our approximation coincides with the halo-model for $\langle\Delta\xi_{+}^{A}(\theta)\Delta\xi_{+}^{B}(\theta)\rangle_{\mathrm{non}\ \mathrm{Gauss.}}$ where $\theta$ is a certain angular scale which we chose to be either our smallest angular bin ($\theta\approx 3^{\prime}$) or a slightly larger scale ($\theta\approx 20^{\prime}$). Note that this is a very crude approximation – even to the log-normal model by Hilbert et al. (2011) since they have not even considered cross-correlations between redshift bins.

We nevertheless use the above matrix as our model covariance $\mathbf{M}$ for the DES shear-shear data vector and compare it to the halo-model covariance $\mathbf{C}$ and the PME. All eigenvalues of matrix $\mathbf{X}$ have in that case $|\lambda_{i}|<1$. The three most dominant eigenvalues are

[TABLE]

in the case where we match the amplitudes of $\mathbf{M}$ and $\mathbf{C}$ at $\theta\approx 20^{\prime}$ and

[TABLE]

when we match the amplitudes at $\theta\approx 3^{\prime}$. In both of these cases the PME in principle converges. However, in the second case at least one eigenvalue comes dangerously close to $1$. In Fig. 9 we show that in terms of the Fisher contours in the $\Omega_{m}$-$\sigma_{8}$-plane the PME nevertheless converges and significantly corrects for the deviations between halo-model and log-normal motivated covariance. We have also checked other parameter combinations and find similar results. The reason that a matching at larger scales gives smaller eigenvalues (i.e. better agreement between halo-model and log-normal motivated covariance) is probably that the scaling of Eq. 35 fails at small scales.

C.2.4 Scale dependent rescaling of the cosmic variance of the multi-probe data vector

Another alternative way to deform the halo-model covariance is to apply different rescaling factors $\alpha$ and $\beta$ for the Gaussian and non-Gaussian cosmic variance parts for different angular scales (cf. Eq. 15). If e.g. the finite area of a survey is not correctly accounted for in a covariance model, the results of Friedrich et al. (2015) indicate that this leads to a scale dependent miss-characterization of the Gaussian cosmic variance and to an almost scale independent over- or underestimation of the non-Gaussian parts. Covariance parts involving shape- or shot-noise on the other hand are less sensitive to the survey area (only to the product of area and galaxy density which is the total number of galaxies).

Motivated by this we replace Eq. 15 by

[TABLE]

where we choose $\beta=0.5$ and $\alpha_{ij}=\sqrt{a_{i}a_{j}}$ setting $a_{i}$ to $1.0$ at the smallest scales and to $0.5$ at the largest scales of the data vector and linearly interpolating for intermediate bins (interpolating in terms of the bin-index).

The most dominant eigenvalues of the matrix $\mathbf{X}$ for this choice of the matrix $\mathbf{M}$ are

[TABLE]

i.e. the PME converges. The bottom panel of Fig. 9 also shows that the Fisher contours derived from the 2nd order PME around this model almost coincide with the ones derived from the halo-model covariance again.

C.3 Convergence in the General Case

Let us now consider the general case. We want invert the equation

[TABLE]

where

[TABLE]

Since both $\mathbf{M}$ and $\mathbf{C}$ are positive definite matrices we can immediately infer that also the matrix $\mathbb{1}+\mathbf{X}$ must be positive definite, i.e. all its eigenvalues must be greater that [math]. As a consequence, all eigenvalues $\lambda_{i}$ of $\mathbf{X}$ must fulfil

[TABLE]

In order to invert $\mathbb{1}+\mathbf{X}$ let us change into the eigenbasis of $\mathbf{X}$ by means of an orthogonal matrix $\mathbf{U}$, i.e.

[TABLE]

It is not a priori clear whether we can invert this by means of the geometric series, since we do not know a priori that $|\lambda_{i}|<1$. As discussed in section 5 in the case that $|\lambda_{i}|>1$ the PME can at least help to identify differences between a covariance model and covariance from (possibly very few) simulations since in that case the 1st order and 2nd order PME will display a divergent behaviour. However, since we know a priori that $\lambda_{i}>-1$ we can in principle apply a trick to let the PME series converge in any case. This trick is to expand $1/(1+\lambda)$ not around $\lambda_{0}=0$ but around some other point $\lambda_{0}=a>0$:

[TABLE]

In terms of the PME series this is in fact equivalent to replacing the model covariance $\mathbf{M}$ by $(1+a)\mathbf{M}$. This way, one can in principle always ensure convergence of the series. This however comes at the expense of the series converging very slowly for eigenvalues of $\mathbf{X}$ that are already close to or smaller than [math]. Since in a real case scenario $\mathbf{M}$ is assumed to be our best guess for the true covariance we hence recommend to stay with $a=0$ and interpret a divergent PME as a significant difference between model and N-body covariance.

Appendix D Data vectors

D.1 Weak lensing data vectors

The redshift distribution and tomographic binning used for our LSST-like weak lensing data vector was chosen to be exactly that of Krause & Eifler (2016, see section 3). This means we assumed an overall source density of $26\mathrm{arcmin}^{-2}$ and a source distribution with a median redshift of $\approx 0.7$ that extends out to $z\gtrsim 3.0$. The tomographic bins were defined by first splitting the redshift distribution into $10$ non-overlapping bins of equal source density and then assuming a Gaussian photoz uncertainty of $\sigma_{z}=0.05$. The intrinsic ellipticity dispersion of the sources was assumed to be $\sigma_{\epsilon}=0.26$ per ellipticity component.

The redshift distribution for the DES-like data vector was chosen to be shallower as for the LSST case reflecting the smaller depth of DES. Here our source distribution has a median redshift of $\approx 0.5$ and extends out to $z=2.0$. The overall source density was taken to be $10/\mathrm{arcmin}^{2}$ and the $5$ tomographic bins where defined assuming a photoz uncertaintz of $\sigma_{z}=0.08$. The intrinsic ellipticity dispersion was chosen to be the same as for the LSST-like case.

D.2 Lens galaxies

For the DES multi-probe data vector we also considered galaxy clustering and galaxy-galaxy lensing correlation functions. For this we were assuming a sample of foreground galaxies with a constant comoving density motivated by the DES redMaGiC sample (Rozo et al., 2016) divided into $3$ tomographic bins whose redshift ranges are $(0.20,0.35)$, $(0.35,0.50)$ and $(0.50,0.65)$. For these galaxies we assumed zero redshift uncertainties motivated by the fact that the redMaGiC redshift errors are small compared to the values for our source samples. The overall density of forground galaxies was taken to be $0.15/\mathrm{arcmin}^{2}$.

D.3 Binning and scales

The real space data vectors use 15 logarithmic angular bin from $\theta=2.5^{\prime}$ to $\theta=250^{\prime}$ for each correlation function and the Fourier space data vector uses 40 logarithmic bin from $\ell=20$ to $\ell=5000$ for each power spectrum. Data vector I contains the correlation functions $\xi_{+}$ and $\xi_{-}$ for each possible combination of source bins. Data vector Ia also contains the auto-correlation of the lens bins and all possible combinations of lens-source correlations (i.e. only those combinations where the sources are at higher redshifts than the lenses). Data vector II contains the auto- and cross-power spectra of all possible combinations of source bins.