Sheer shear: weak lensing with one mode

TL;DR

This paper demonstrates that a single radial eigenmode derived from Karhunen-Loève decomposition captures nearly all cosmological information in cosmic shear data, greatly simplifying data analysis and covariance computation.

Contribution

The study shows that one radial eigenmode suffices for cosmological constraints in shear surveys, reducing data dimensionality and computational complexity significantly.

Findings

Almost all cosmological information is contained in one eigenmode.

Data compression reduces covariance computation by a factor of 30 to 800.

Eigenfunctions are nearly ell-independent, enabling simplified analysis.

Abstract

3D data compression techniques can be used to determine the natural basis of radial eigenmodes that encode the maximum amount of information in a tomographic large-scale structure survey. We explore the potential of the Karhunen-Lo\`eve decomposition in reducing the dimensionality of the data vector for cosmic shear measurements, and apply it to the final data from the \cfh survey. We find that practically all of the cosmological information can be encoded in one single radial eigenmode, from which we are able to reproduce compatible constraints with those found in the fiducial tomographic analysis (done with 7 redshift bins) with a factor of ~30 fewer datapoints. This simplifies the problem of computing the two-point function covariance matrix from mock catalogues by the same factor, or by a factor of ~800 for an analytical covariance. The resulting set of radial eigenfunctions is close to ell-independent, and therefore they can be used as redshift-dependent galaxy weights. This simplifies the application of the Karhunen-Lo\`eve decomposition to real-space and Fourier-space data, and allows one to explore the effective radial window function of the principal eigenmodes as well as the associated shear maps in order to identify potential systematics. We also apply the method to extended parameter spaces and verify that additional information may be gained by including a second mode to break parameter degeneracies. The data and analysis code are publicly available at https://github.com/emiliobellini/kl_sample.

Tables

Table 1. Table 1: Redshift bins used in this analysis along with the effective density n eff subscript 𝑛 eff n_{\rm eff} per bin and field. The last column shows the total n eff subscript 𝑛 eff n_{\rm eff} per bin.

Bin No.	$z_B$	$n_{\mathrm{eff}}$ [arcmin-2]
		W1	W2	W3	W4	TOT
1	$0.15$ – $0.29$	$1.41$	$1.38$	$1.39$	$1.37$	$1.40$
2	$0.29$ – $0.43$	$1.18$	$1.23$	$1.18$	$1.12$	$1.18$
3	$0.43$ – $0.57$	$1.16$	$1.17$	$1.30$	$1.19$	$1.21$
4	$0.57$ – $0.70$	$1.70$	$1.57$	$1.84$	$1.80$	$1.73$
5	$0.70$ – $0.90$	$2.69$	$2.58$	$2.91$	$2.87$	$2.76$
6	$0.90$ – $1.10$	$1.96$	$1.49$	$1.91$	$1.67$	$1.85$
7	$1.10$ – $1.30$	$1.07$	$0.89$	$0.94$	$0.83$	$0.98$

Table 2. Table 2: We show the ℓ ℓ \ell -bandpowers used in this analysis. The first column indicates the number of the bandpower, the second one the ℓ ℓ \ell -range that each bandpower covers and the last column indicates which bandpowers were actually used to get cosmological constraints.

No.	$\ell$-range	Used
1	30 – 80
2	80 – 260	✓
3	260 – 450	✓
4	450 – 670	✓
5	670 – 1310	✓
6	1310 – 2300	✓
7	2300 – 5100

Table 3. Table 3: Cosmological parameter priors used in this analysis. All the priors are uniform.

Parameter	Prior
$h$	$0.61$ – $0.81$
${\Omega }_b{h}^2$	$0.013$ – $0.033$
${\Omega }_c{h}^2$	$0.001$ – $0.99$
$\ln \left({10}^{10}{A}_s\right)$	$2.3$ – $5.0$
$n_s$	$0.7$ – $1.3$

Table 4. Table 4: We show the size of the data vector along with the Figure of Merit (FoM) for each case considered here. The first column indicates the KL-compression used. The second column shows the size of the data vector for each method, where s ℓ ​ θ = 5 subscript 𝑠 ℓ 𝜃 5 s_{\ell\theta}=5 in Fourier space ( 5 5 5 ℓ ℓ \ell bandpowers) and s ℓ ​ θ = 10 subscript 𝑠 ℓ 𝜃 10 s_{\ell\theta}=10 in Real space ( 6 6 6 angles θ 𝜃 \theta for ξ + subscript 𝜉 \xi_{+} and 4 4 4 for ξ − subscript 𝜉 \xi_{-} ). The other columns show the FoM for each case, calculated as the Area -1 of the 1 ​ σ 1 𝜎 1\sigma contour in Ω m subscript Ω 𝑚 \Omega_{m} - S 8 subscript 𝑆 8 S_{8} . “Fourier SI” refers to the analysis in Fourier Space with ℓ ℓ \ell -independent KL-transform, “Fourier SD” refers to the ℓ ℓ \ell -dependent KL-transform, while “Real” shows the KL-transformed results in real space.

KL	Size	Fourier SI	Fourier SD	Real
full	$28\times s_{\ell \theta }$	$109.7$	$109.8$	$83.0$
6_off	$21\times s_{\ell \theta }$	$103.0$	$104.6$	$82.1$
5_off	$15\times s_{\ell \theta }$	$107.7$	$106.2$	$85.4$
4_off	$10\times s_{\ell \theta }$	$102.7$	$109.0$	$83.1$
3_off	$6\times s_{\ell \theta }$	$104.2$	$100.1$	$85.2$
2_off	$3\times s_{\ell \theta }$	$105.8$	$103.1$	$89.5$
7_diag	$7\times s_{\ell \theta }$	$104.5$	$102.9$	$81.0$
6_diag	$6\times s_{\ell \theta }$	$104.1$	$114.0$	$73.0$
5_diag	$5\times s_{\ell \theta }$	$105.0$	$102.3$	$76.4$
4_diag	$4\times s_{\ell \theta }$	$106.6$	$115.9$	$72.4$
3_diag	$3\times s_{\ell \theta }$	$99.9$	$111.8$	$76.6$
2_diag	$2\times s_{\ell \theta }$	$98.9$	$101.1$	$81.7$
1	$1\times s_{\ell \theta }$	$98.8$	$107.0$	$83.0$

Equations

$$({\gamma }_1+i{\gamma }_2{)}_{\alpha }=-ðð\left[{\int }_0^{{\chi }_H}d\chi \Phi (\chi \hat{\mathbf{n}})q_{\alpha }(\chi )\right],$$

$$q_{\alpha }(\chi )=\chi {\int }_{z(\chi )}^{\infty }dz\left(1-\frac{\chi }{\chi (z)}\right)n_{\alpha }(z).$$

$$ð{}_{s}f\equiv -(\sin \theta {)}^s\left(\frac{\partial }{\partial \theta }+\frac{i}{\sin \theta }\frac{\partial }{\partial \phi }\right)(\sin \theta {)}^{-1}{}_{s}f.$$

$$S_{\ell }^{\alpha \beta }=\int \frac{d\chi }{{\chi }^2}{\tilde{q}}_{\alpha }(\chi ){\tilde{q}}_{\beta }(\chi )P\left(k=\frac{\ell +1/2}{\chi },z(\chi )\right)$$

$${\tilde{q}}_{\alpha }(\chi )\equiv F_{\ell }\frac{3H_0^2{\Omega }_M}{2a(\chi )}{q}_{\alpha }(\chi ),F_{\ell }\equiv \sqrt{\frac{(\ell +2)!}{(\ell -2)!}}\frac{1}{(\ell +1/2{)}^2},$$

$$C_{\ell }^{\alpha \beta }=S_{\ell }^{\alpha \beta }+N_{\ell }^{\alpha \beta }.$$

$${\xi }_{\pm }^{\alpha \beta }(\theta )=\int \frac{d\ell \ell }{2\pi }{J}_{2\mp 2}(\ell \theta )C_{\ell }^{\alpha \beta },$$

$${\partial }_{\theta }\mathsf{C}{\mathbf{e}}_p={\lambda }_p\mathsf{C}{\mathbf{e}}_p,$$

$$\mathsf{C}{\mathbf{e}}_p={\lambda }_p\mathsf{N}{\mathbf{e}}_p,$$

$$\left[{\mathsf{L}}^{-1}\mathsf{C}({\mathsf{L}}^{-1}{)}^{†}\right]{\tilde{\mathbf{e}}}_p={\lambda }_p{\tilde{\mathbf{e}}}_p.$$

$$n_{\mathrm{eff}}=\frac{{\left({\sum }_i{w}_i\right)}^2}{A_{\mathrm{eff}}\left({\sum }_i{w}_i^2\right)}.$$

$${\tilde{\gamma }}_i({\hat{\mathbf{n}}}_p)=\frac{{\sum }_{g\in p}{w}_g\left(e_{i,g}-c_{i,g}\right)}{\left(1+m_p\right){\sum }_{g\in p}{w}_g},$$

$$⟨{\tilde{C}}_{\ell }^{fg}⟩\equiv ⟨\frac{1}{2\ell +1}\sum _{m=-\ell }^{\ell }\mathrm{Re}(f_{\ell m}^v{g}_{\ell m}^{w\ast })⟩=\sum _{{\ell }^{\prime }}{M}_{\ell {\ell }^{\prime }}^{vw}{C}_{{\ell }^{\prime }}^{fg},$$

$${\tilde{\mathcal{C}}}_k^{fg}=\sum _{\ell ={\ell }_{\min }^k}^{{\ell }_{\max }^k}{\tilde{C}}_{\ell }^{fg}.$$

$${\mathcal{M}}_{kq}^{vw}\equiv \sum _{\ell ={\ell }_{\min }^k}^{{\ell }_{\max }^k}\sum _{{\ell }^{\prime }={\ell }_{\min }^q}^{{\ell }_{\max }^q}{M}_{\ell {\ell }^{\prime }}^{vw}$$

$${\mathcal{C}}_k^{fg}={\left({\mathcal{M}}^{vw}\right)}_{kq}^{-1}\sum _{\ell ={\ell }_{\min }^q}^{{\ell }_{\max }^q}({\tilde{C}}_{\ell }^{fg}-{\tilde{N}}_{\ell }^{fg}),$$

$${\mathcal{C}}_{\mathrm{tot}}={\mathsf{\Sigma }}_{\mathrm{tot}}\sum _{f=1}^4{\mathsf{\Sigma }}_f^{-1}{\mathcal{C}}_f,{\mathsf{\Sigma }}_{\mathrm{tot}}^{-1}=\sum _{f=1}^4{\mathsf{\Sigma }}_f^{-1}.$$

$$\mathsf{\Sigma }=\frac{1}{N_s-1}\sum _{i=1}^{N_s}\left({\vec{\mathbf{C}}}_i-⟨\vec{\mathbf{C}}⟩\right){\left({\vec{\mathbf{C}}}_i-⟨\vec{\mathbf{C}}⟩\right)}^T,$$

$${\bar{e}}_p^{\alpha }=\frac{{\sum }_{\ell }(2\ell +1)e_{p,\ell }^{\alpha }}{{\sum }_{\ell }(2\ell +1)},$$

Full text

1 Introduction

Imaging galaxy surveys are one of the most promising paths towards understanding some of the main open questions in cosmology (Hu, 1999; Amara and Réfrégier, 2007), such as the precise nature of the energy components dominating the late-time expansion (Linder, 2003), the level of non-Gaussianity in the cosmic inhomogeneities induced by early-time physics (Dalal et al., 2008; Alonso and Ferreira, 2015), or the mass of neutrinos (Hu et al., 1998). To do so, these experiments combine the constraining power of several cosmological probes, such as Type-Ia supernovae, strong lenses, the distribution of galaxy clusters, the clustering of galaxies, and the weak gravitational lensing effect through the correlated distortion of galaxy shapes.

Current and near-future photometric surveys, such as the Dark Energy Survey (DES Collaboration et al., 2017), the Kilo-Degree Survey (Hildebrandt et al., 2018), the Hyper Suprime-Cam (Hikage et al., 2018), Euclid (Laureijs et al., 2011) or the Large Synoptic Survey Telescope (LSST) (LSST Science Collaboration et al., 2009) face a number of observational, statistical and computational challenges in order to extract the maximum amount of cosmological information from them without incurring significant systematic biases. Arguably the most critical challenge for these experiments is the correct treatment of the photometric redshifts (Newman et al., 2015; Hildebrandt et al., 2017). An object by object analysis would be expensive and the broad nature of the lensing kernel dilute the information contained by each galaxy. In addition, we lack perfectly accurate posterior redshift probability distributions for each object. This is why a so-called tomographic approach (Hu, 1999) has become the standard analysis path. In this case, the full galaxy sample is divided into photometric redshift (photo-$z$) slices, and cosmological constraints are extracted from the different cross-correlations between the 2D galaxy positions and shapes in all pairs of redshift bins. In this approach, the problem of obtaining and interpreting accurate posterior redshift distributions is replaced by the potentially simpler problem of estimating the redshift distribution of all objects in each redshift bin. The number and distribution of the different redshift bins used in the analysis is then determined as a compromise between the needs to maximize the amount of information recovered and to minimize the size of the final data vector (in the form of a collection of all auto- cross-bin two-point functions).

In Alonso (2018), we proposed a method to determine the optimal slicing for a photometric survey that recovers the maximum amount of information in the minimum number of radial modes. The method is based on the Karhunen-Loève transform, and follows a formalism similar to other 3D data compression methods111In this paper, the use of the term “3D data compression” refers to the decomposition of a three-dimensional dataset, such as galaxy shear, into the set of modes spanned by a given set of radial basis functions. proposed in the literature (e.g. Heavens 2003; Kitching et al. 2007; McEwen and Leistedt 2013; Khalid et al. 2014; Leistedt et al. 2015; Asgari and Schneider 2015). The resulting modes are in general scale- and redshift-dependent, and in the case of spectroscopic galaxy clustering can be related to the standard harmonic-Bessel transform that leads to the standard 3D power spectrum analysis, see a.g. Alonso (2018). In the case of cosmic shear, due to the radially cumulative nature of gravitational lensing, the Karhunen-Loève transform can produce a significant reduction in the dimensionality of the final data vector, which has a direct impact on the computational complexity of estimating its covariance matrix (beyond guaranteeing an optimal recovery of cosmological information).

Two main approaches have been used in the literature to attack the problem of covariance matrices for large-scale structure and weak lensing data. On the one hand, the sample covariance matrix can be computed by running a large number of simulations resembling the data, estimating the data vector in each of them, and estimating their sample variance. In order to recover a converged covariance that can be safely inverted, this requires running a number of simulations several times larger than the size of the data vector (Dodelson and Schneider, 2013; Petri et al., 2016) ($\mathcal{O}(10)$ times larger is often a good rule of thumb). Since the size of a complete set of tomographic cross-correlations scales as $\sim N_{\rm bin}^{2}$ with the number of tomographic bins, the computational complexity of estimating the sample covariance can grow very fast if a large number of bins is required in order to recover the cosmological information. On the other hand, it is also possible to compute analytical estimates of the covariance matrix. These analytical covariances (Krause and Eifler, 2017) are often faster to compute than the simulated ones, but they can also be computationally challenging, since it is often necessary to accurately account for the mode-coupling caused by the survey mask (e.g. see García-García et al., 2019), and real-space covariances can require expensive double Hankel transforms. Since all elements of the covariance matrix must be estimated individually, the complexity for analytical covariances scales with the square of the size of the data vector, or $\sim N_{\rm bin}^{4}$. The problem of estimating covariances can therefore benefit from efficient compression methods (Heavens et al., 2017).

As a proof of concept for the potential of this method, in this paper we apply it to the analysis of cosmic shear data from the Canada-France-Hawaii Telescope Lensing Survey (hereafter refered to as CFHTLenS, Heymans et al. 2012). Even if more up-to-date cosmic shear data are available, the choice of CFHTLenS has been made mainly because other 3D data compression methods have been applied in the past (Kitching et al., 2014), which makes possible a comparison between different methods. We build upon the tomographic analysis of Joudaki et al. (2017) and present the performance of our method in terms of optimal radial eigenmodes, compression factor, final constraints and goodness of fit. The paper is structured as follows: in Section 2 we provide the theoretical background for weak lensing tomographic analyses and the Karhunen-Loève transform as a data compression method. Section 3 then introduces the data used in our analysis, as well as the different data products we derive from it. Our main results are presented in Section 4, which we discuss and summarize in Section 5.

2 Methods

2.1 Cosmic shear 2-point functions

The cosmic shear field for galaxies observed in a given tomographic redshift bin $\alpha$ characterized by a redshift distribution $n_{\alpha}(z)$ is given by

[TABLE]

where $\Phi(\chi\hat{\bf n})$ is the gravitational potential in the lightcone towards a sky position with spherical coordinates $(\theta,\varphi)$ defined by the unit vector $\hat{\bf n}$, $\chi_{H}$ is the comoving radial distance to the horizon. The lensing kernel $q_{\alpha}$ in flat space becomes

[TABLE]

Here, $z$ is redshift, and $\eth$ is the spin-raising differential operator, which acts on a spin-$s$ quantity defined on the sphere $\,{}_{s}f(\hat{\bf n})$ as222Note that the action of the first $\eth$ operator on the spin-0 quantity $\Phi$ makes it a spin-1 field, which affects the action on it of the second $\eth$ operator.

[TABLE]

The harmonic decomposition of a spin-2 quantity like $\gamma$ gives rise, in general, to both $E$- and $B$-modes (e.g. see Zaldarriaga and Seljak, 1997). Since gravitational lensing only gives rise to negligible $B$-modes, we will only concern ourselves with the $E$-mode component here (although our data analysis pipeline computes the $B$-mode component as a test of systematics).

Since the gravitational lensing that shears the observed shape of a galaxy at redshift $z$ is caused by the collective effect of all the matter inhomogeneities along the line of sight to that redshift (manifested by the cumulative nature of the lensing kernels $q_{\alpha}$ in Eq. 2), the weak lensing signal in different redshift bins is usually strongly correlated. The effective number of degrees of freedom in a tomographic analysis therefore increases very slowly with the number of redshift bins, which motivates the use of radial data compression, as described in Section 2.2, to determine the optimal radial modes that saturate the information content of a given weak lensing survey.

The two-point correlator of the weak lensing $E$-mode signal in two redshift bins $\alpha$ and $\beta$ at multipole $\ell$ can be related to the power spectrum of matter fluctuations $P(k,z)$ as (Kilbinger et al., 2017)

[TABLE]

where $S$ stands for Signal, see Eq. 5 below, and

[TABLE]

$a(\chi)$ is the scale factor at the comoving distance $\chi$, $H_{0}$ is the Hubble constant and $\Omega_{M}$ is the matter abundance at $z=0$. In deriving Eq. 4 we have used the Limber approximation333Note that the Limber approximation is accurate enough for weak lensing studies on the scales used here (Kilbinger et al. 2017). (Limber, 1953; Afshordi et al., 2004), and assumed a flat cosmology. Under the assumption that the observed shear is a combination of the true underlying signal and stochastic noise, the total power spectrum is also a sum of the signal and noise contributions:

[TABLE]

We assume the noise power spectrum to be uncorrelated across bins (i.e. $N^{\alpha\beta}_{\ell}\propto\delta_{\alpha\beta}$), and we estimate it as described in Section 3.

The shear angular power spectrum can then be related to the real-space correlation functions $\xi_{\pm}$ as

[TABLE]

where $J_{\nu}(x)$ is the cylindrical Bessel function of order $\nu$, and we have used the flat-sky approximation, which is appropriate for the small angular footprint covered by the CFHTLenS survey, i.e. $\sim 154\,{\rm deg}^{2}$.

To estimate the theoretical predictions for all two-point functions used in this work we used the Core Cosmology Library (CCL) (Chisari et al., 2018), with the matter power spectrum computed by CLASS (Lesgourgues, 2011) using the version of Halofit by Takahashi et al. (2012).

2.2 The Karhunen-Loève transform and radial data compression

Suppose we have a data vector ${\bf x}$ with dimension $N_{s}$, that we assume Gaussianly distributed with zero mean, and whose covariance ${\sf C}\equiv\langle{\bf x}{\bf x}^{\dagger}\rangle$ depends on a particular parameter $\theta$ we want to measure. It is possible to find a linear transformation ${\sf E}$ such that the elements of the transformed data vector ${\bf y}\equiv{\sf E}^{\dagger}\,{\bf x}$ have unit variance and are uncorrelated (i.e. $\langle y_{p}y_{q}^{*}\rangle=\delta_{pq}$), and such that the first $m<N_{s}$ elements contain most of the information about $\theta$ (Tegmark et al., 1997; Bond, 1995; Vogeley and Szalay, 1996). The columns of ${\sf E}$, which we label ${\bf e}_{p}$, are the eigenvectors of the following generalized eigenvalue problem:

[TABLE]

where $\partial_{\theta}$ denotes the partial derivative with respect to $\theta$. The resulting linear transformation is the so-called Karhunen-Loève (KL) transform, which can be used to simplify the analysis for one particular parameter.

In a more general situation, we would like to measure several parameters simultaneously. Although under certain circumstances it may be possible to find a set of KL eigenmodes that maximize the overall information content, as described in Tegmark et al. (1997), a simpler and popular approach is to instead maximize the information on the overall amplitude of the signal. For concreteness, let us decompose our data into signal and noise contributions ${\bf x}={\bf s}+{\bf n}$, and let us introduce a fictitious amplitude $\alpha$ with fiducial value 1 such that in general ${\bf x}=\alpha{\bf s}+{\bf n}$. Under the assumption that signal and noise are uncorrelated, and redefining $2/(2-\lambda_{p})\rightarrow\lambda_{p}$, we can rewrite the generalized eigenvalue problem in Eq. 7 as

[TABLE]

where ${\sf C}={\sf S}+{\sf N}$ is the covariance matrix of ${\bf x}$ as a sum of the signal and noise covariances. Using the Cholesky decomposition of the noise covariance ${\sf N}={\sf L}{\sf L}^{\dagger}$, and defining $\tilde{\bf e}_{p}\equiv{\sf L}^{\dagger}{\bf e}_{p}$, the previous equation can be rewritten as a standard eigenvalue problem for $\tilde{\bf e}_{p}$:

[TABLE]

Solving for $\tilde{\bf e}_{p}$, the new data vector can then be computed as $y_{p}=\tilde{\bf e}_{p}^{\dagger}\,{\sf L}^{-1}{\bf x}$. The covariance of the new data vector is given by the eigenvalues $\langle y_{p}y^{*}_{p^{\prime}}\rangle=\lambda_{p}\delta_{pp^{\prime}}$, and the transformed noise covariance is simply the identity.

As described in Alonso (2018), this formalism can be used to determine the optimal set of radial modes in a given photometric redshift survey. In this case, our data vector ${\bf x}$ is the vector of $E$-mode harmonic coefficients of the cosmic shear field in different redshift bins for a fixed $(\ell,m)$, and ${\sf S}$, ${\sf N}$ and ${\sf C}$ are the signal, noise and total power spectra described in the previous section at the corresponding value of $\ell$. The KL eigenmodes $e^{\alpha}_{p,\ell}$ are therefore scale-dependent functions in general. However, as shown in Alonso (2018), this scale dependence is very mild for cosmic shear data, and in this case the eigenmodes can be thought of as redshift-dependent galaxy weights that can be used to generate the set of most informative shear maps. We will explore this in detail in Section 4. In this case, each of the KL modes can be interpreted as a shear map associated with an effective redshift distribution given by $\tilde{n}_{p}(z)=e^{\alpha}_{p}\,n_{\alpha}(z)$.

3 Data

3.1 The CFHTLenS lensing sample

We use data from the CFHTLenS survey444http://www.cfhtlens.org/ (Heymans et al., 2012). Their analysis performed weak lensing data processing with THELI (Erben et al., 2013), shear measurement with lensfit (Miller et al., 2013), and photometric redshift measurement with PSF-matched photometry (Hildebrandt et al., 2012). A full systematic error analysis of the shear measurements in combination with the photometric redshifts is presented in Heymans et al. (2012), with additional error analyses of the photometric redshift measurements presented in Benjamin et al. (2013).

CFHTLenS is a deep multi-colour survey that spans $154\,{\rm deg}^{2}$ of the sky divided in four fields – labelled W1, W2, W3 and W4. The full CFHTLenS catalogue contains $\sim 2\times 10^{7}$ objects (precisely 22,652,728) that have been classified as “galaxy”, “star” or “bad fit object” by the lensfit algorithm (see Miller et al. 2013 for details). Furthermore, some of the objects are in masked regions of the sky due to observational systematics, or their position is in MegaCam fields that did not pass the lensing residual systematics effects, as described in Heymans et al. (2012). Thus, in any shear analysis that uses CFHTLenS data, the catalogue has to be cleaned from objects that are classified as stars or lie on unreliable sky regions. In detail, the filters we used in order to obtain a clean catalogue are:

• •

MASK$=0$: to select only unmasked objects;

• •

weight$>0$: objects that have been measured with strictly positive weight;

• •

star_flag$=0$: remove stars and bad fits, and keep only galaxies;

• •

remove objects with positions in the $1\,{\rm deg}^{2}$ camera MegaCam fields that do not pass the lensing residual systematics effects. For a list of those fields see Erben et al. (2013).

After using these filters, the clean catalogue reduces to 4,180,003 objects, which we use to estimate the lensing two-point statistics.

The ellipticities of each source, as well as their additive ($c_{i}$) and multiplicative ($m_{i}$) correction biases and an approximately optimal inverse variance weight ($w_{i}$) are obtained using the Bayesian model fitting software lensfit (Miller et al. 2013). Photometric redshifts are derived from the Bayesian photometric redshift code BPZ (Benítez 2000), which returns a full redshift probability distribution $p(z)$ for each object. As in Joudaki et al. (2017), we split the sample into 7 tomographic bins, using the peak of the redshift distribution of each galaxy ($z_{B}$) as a redshift bin label. Table 1 shows the redshift bins chosen along with the effective density of galaxies in each bin, defined as

[TABLE]

Here the effective area ($A_{\rm eff}$) per field is [41.0, 11.6, 24.7, 12.2] deg*-2* for W1-4 respectively. The redshift distribution for each tomographic bin was obtained, as in Joudaki et al. (2017), by stacking the photo-$z$ $p(z)$ of all unmasked objects in that bin. The resulting normalized redshift distributions are shown in Figure 1.

3.2 Masks and maps

In this section we describe in detail the method used to generate masks and maps of shear measurements, which are then used to compute the power spectra. As mentioned, the relative small area covered by CFHTLenS allows to use a flat-sky approximation. In this case, masks and maps are 2D arrays of pixels (one array for each field). In each pixel, the masks quantify the quality of the CFHTLenS measurements in a given sky area, while the maps store the value of the two shear components estimated by averaging over the ellipticities of all sources in the pixel.

To define our sky mask, we start from a set of individual sky masks for each of the 4 CFHTLenS fields. These have a resolution of 1 arcsecond and allow us to remove problematic areas, such as those outside the CFHTLenS footprint, around stars and stellar halos, significant object overdensities, flagged pixels, and others. The full description of the masked areas can be found in Table B2 of Erben et al. (2013). While in Erben et al. (2013) the authors mention that it is possible to consider as unmasked areas those that lie in large masks around stars and stellar haloes (mask value of $1$), here we adopt a more conservative approach, and consider only areas with mask value of [math]. We also remove by hand all pixels falling inside any of the pointings identified by Heymans et al. (2012) as being systematics dominated (corresponding to a total of $\sim 25\%$ of the available area). Once all sources in masked pixels are removed from the data, we downgrade the resolution of this binary mask to a pixel size of 2 arcminutes, while preserving the information in the finer map in terms of the fractional area masked in each larger pixel. We chose this resolution to guarantee that the majority of the unmasked pixels will be populated by, on average, four galaxies in each redshift bin. Finally, to optimally weight each pixel when computing the angular power spectra, for each field and redshift bin we define the weights map used to compute the power spectrum (see Section 3.3) to be proportional to the sum of the lensfit weights of all sources in each pixel. Figure 2 shows the mask and weights map used in this analysis for a particular field (W3). The top panel shows the 2 arcminutes mask obtained by downgrading the original 1 arcsec binary mask. This is the mask that has been used to calculate the area of each field used in Eq. 10. The bottom panel represents the weights map described above, used by the power spectrum estimator. It is worth noting that the weights map is therefore different for each redshift bin, and the associated mode-coupling matrix (described in the next section) must be computed independently for each pair of bins.

To compute power spectra we also need to produce maps of the two shear components $\gamma_{1}(\hat{\bf n})$ and $\gamma_{2}(\hat{\bf n})$. A first estimate of these maps is made as a weighted average of the ellipticity components of all sources falling within each pixel corrected for their measured additive and multiplicative shape measurement bias $c_{i,g}$ and $m_{p}$

[TABLE]

where $g$ stands for galaxy and $p$ for pixel. While the additive correction can be applied on a object by object basis (noting that $c_{1}=0$ in the CFHTLenS analysis), $m$ must be estimated by averaging over a large number of sources to avoid a large bias in the final maps (Miller et al., 2013). Thus, we make a map of $m_{p}$ for each field and redshift bin by averaging, for a given pixel $p$, over the multiplicative bias parameter of all sources lying within a square of size 10 arcmin centred around $p$. In Figure 3 we show the multiplicative correction map for the W3 field and the sixth redshift bin (the bin with the highest shear signal-to-noise ratio). It is possible to notice that this gives up to a $\sim 10\%$ correction to the shear measurements. Finally, Figure 4 shows the shear maps used in this analysis for the W3 field and the sixth redshift bin, with $\gamma_{1}$ and $\gamma_{2}$ shown in the top and bottom panels respectively.

3.3 Two-point functions

We estimate all $EE$ and $BB$ cross-power spectra between the 7 different redshift bins using a pseudo-$C_{\ell}$ estimator as implemented in NaMaster Alonso et al. (2018). We describe the estimator briefly here and refer the reader to the original paper for further details.

The pseudo-$C_{\ell}$ estimator attempts to make an unbiased measurement of the power spectrum of a masked field by analytically computing the coupling between modes induced by that mask. Mathematically, let $f^{v}(\hat{\bf n})=v(\hat{\bf n})\,f(\hat{\bf n})$ be the observation of an underlying true field $f$ through a multiplicative weights map $v$. Minimally, this weights map will be a binary mask removing regions where the field has not been observed ($v=0$) leaving only the observed pixels ($v=1$). More optimally, the weights map should be proportional to the local inverse noise variance, downweighting or upweighting different regions depending on their noise properties. In our analysis, the weights map is given by the sum of lensfit weights in each pixel.

Due to the convolution theorem, a direct estimate of the power spectrum from the masked fields $\tilde{C}^{fg}_{\ell}$ ignoring their masks, will in general yield a mode-coupled version of the underlying true power spectrum

[TABLE]

where the mode-coupling matrix $M^{vw}_{\ell\ell^{\prime}}$ depends on the masks $v$ and $w$, and can be estimated analytically. The mode-coupling matrix is usually ill-conditioned and cannot be directly inverted to yield an unbiased estimate of $C^{fg}_{\ell}$. The usual approach is then to bin the pseudo-$C_{\ell}$ into bandpowers:

[TABLE]

For this work we use the same $\ell$ bins chosen by Köhlinger et al. (2016)555Note that, in order to obtain stronger cosmological constraints, in this paper we use one additional bandpower w.r.t. Köhlinger et al. (2016), namely the 6th, i.e. $\ell=1310$ – $2300$. and summarized in Table 2. This choice was motivated by the need to capture the most relevant power spectrum features while staying within the mildly non-linear regime. For sufficiently wide bandpowers, the corresponding binned coupling matrix

[TABLE]

becomes invertible, and the final estimator takes the form

[TABLE]

where $\tilde{N}^{fg}_{\ell}$ is an estimate of the noise pseudo-power spectrum (the so-called “noise bias”). This estimator, described here for scalar fields on the sphere, can be easily generalized to spin-2 quantities and flat skies, which we use in this work. We refer the reader to Alonso et al. (2018) for further details.

Finally, as in Hikage et al. (2018), we estimate the noise bias in Eq. 15 by averaging over the power spectrum measured in 1000 noise realizations, generated by randomly rotating the ellipticities of all sources in the catalogue.

This procedure allows us to estimate the power spectra in each of the 4 CFHTLenS fields. The combined power spectra are then computed as an inverse-variance sum of the per-field power spectra:

[TABLE]

Here ${\cal C}_{f}$ is the vector of all power spectra computed in field $f$, $\mathsf{\Sigma}_{f}$ is the power spectrum covariance matrix in field $f$, described in the next section, and ${\cal C}_{\rm tot}$ and $\mathsf{\Sigma}_{\rm tot}$ are the combined power spectra and covariance matrix.

The resulting measurements of the $E$-mode power spectra estimated from the CFHTLenS data are shown in Figure 5, where the error bars are computed from the covariance matrix described in Section 3.4. We find good agreement with the best-fit cosmological model of Joudaki et al. (2017) (solid lines). It is worth noting that this is the first time that a pseudo-$C_{\ell}$ analysis of these data has been carried out. Joudaki et al. (2017) used real-space correlations measured on these same redshift bins, while Köhlinger et al. (2016) used an optimal quadratic estimator on two redshift bins that would be computationally impractical for a larger number of subsamples. After implementing the scale cuts shown in Figure 5, the final data vector has 140 elements.

We use three sanity checks to validate the power spectrum estimation. First, running our pipeline on the set of 2000 Gaussian simulations used to calculate the covariance matrix (described in Section 3.4), we verify that we are able to recover the input power spectra with no appreciable bias (this was also thoroughly verified in Alonso et al. 2018). Secondly, we estimate the $B$-mode power spectrum of all pairs of redshift bins in the CFHTLenS data. Since the cosmological signal is not expected to yield a detectable $B$-mode signal, the presence of $B$-modes can be used as a diagnostic for unknown systematics. Within the range of scales used in this analysis we do not detect significant $B$-modes, with a probability to exceed (PTE666The PTE here is the probability of the $\chi^{2}$ to be larger than the value found in the data.) associated with the measurements shown in Figure 6 of $\sim 10\%$. Lastly, with our pipeline we replicated the analysis done in Köhlinger et al. (2016), and the resulting $E$- $B$-modes match perfectly. Other $B$-mode analyses of the CFHTLenS data have been carried out in Kitching et al. (2014), Asgari et al. (2017) and Alsing et al. (2017).

Although our main analysis uses Fourier-space observables, we also present results for the measurements of the correlation function made by Joudaki et al. (2017), in order to make a direct comparison with those results. The measurements correspond to the 7 tomographic bins from $z=0.15-1.3$ and 7 angular bins from $\theta=1-120$ arcmin, which post-masking results in a data vector of 280 elements. The correlation functions $\xi_{\pm}$ were computed for each pair of redshift bins in each of the four fields using the athena software777http://www.cosmostat.org/software/athena (Kilbinger et al., 2014). These measurements include the same additive shear correction as in the Fourier space analysis, and moreover include a weighted ensemble average multiplicative shear correction following Heymans et al. (2012) and Miller et al. (2013). The patch-dependent measurements were subsequently averaged by the number of galaxy pairs returned by athena.

3.4 Covariance matrix

To estimate the covariance matrix of the power spectrum measurements we make use of 2000 Gaussian simulations of the 4 fields. Each simulation is generated in two steps:

We generate a Gaussian realization of the expected cosmic shear signal $\gamma$. To do so, we draw Gaussian realizations of the Fourier-space $E$-modes with a variance given by the $E$-mode power spectrum predicted by the best-fit cosmological parameters of Joudaki et al. (2017), and no $B$-modes. These are then transformed into real space to produce a shear map. Note that the Gaussian Fourier coefficients are not drawn independently for each redshift bin, but fully preserve the expected correlation between them. 2. 2.

We then interpolate the Gaussian shear field to the position of each source in the original catalogue. The total ellipticity of each source is then calculated by adding, to the interpolated shear value, the galaxy’s true ellipticity rotated by a random angle. This last steps allows us to include shape noise in a realistic way (under the assumption that the ellipticity of a single object is dominated by noise, and not by the lensing-induced shear).

Each simulated catalogue is then put through the same power spectrum pipeline used for the real data, including masking, power spectrum estimation and field coaddition. We use the power spectra measured in each simulation to estimate the sample covariance matrix as:

[TABLE]

where $N_{s}$ is the number of simulations used to estimate the covariance, $\vec{\bf C}_{i}$ is the vector of power spectra measured in the $i$-th simulation and $\langle\cdot\rangle$ denotes averaging over simulations. When using this covariance to compute Gaussian likelihoods, we correct the inverse covariance by the multiplicative ‘Hartlap factor’ introduced by Hartlap et al. (2007) to account for the finite number of simulations used. As shown by Sellentin and Heavens (2016), this approximation is sufficiently accurate given the large number of simulations used in our analysis. This was explicitly confirmed in Joudaki et al. (2017).

Note that, as described in the next section, we use a maximum of $N_{s}=2000$ simulations when analysing the full set of power spectra, but we scale this number down with the size of the final data vector when exploring cases with different numbers of KL modes. This is to showcase the usage of data compression to reduce the number of mock realizations needed to estimate the covariance matrix while retaining the final constraining power. In detail, the number of simulations used for each case is chosen to keep the ‘Hartlap factor’ constant.

We neglect all non-Gaussian contributions to the covariance matrix. The only relevant contribution for weak-lensing measurements, as quantified by Barreira et al. (2018) is the so-called super-sample covariance, which becomes less important for low-density samples where the measurements on high-$\ell$s are noise-dominated. This fact, together with the acceptable value of the $\chi^{2}$ we obtain both for the fiducial and best-fit cosmological parameters (see Section 4) imply that the covariance matrix estimated with this method is appropriate for our analysis. It is worth noting that our objective is not to provide state-of-the-art constraints on cosmological parameters from these data, but rather to showcase the performance of the KL decomposition in terms of data compression. The resulting $E$-mode power spectrum covariance is shown in Figure 7. For a given auto- or cross-spectrum, the covariance is strongly dominated by the diagonal.

In order to compare our results with those of Joudaki et al. (2017), we follow their estimation of the covariance matrix for the real-space correlation functions as closely as possible. In this case, the covariance matrix is estimated from a large suite of 497 distinct Scinet LIght Cone Simulations (Harnois-Déraps and van Waerbeke, 2015) based on a WMAP9$+$BAO$+$SN cosmology. Following $2\times 2$ sub-division of the simulation boxes, the resulting set of 1988 pseudo-independent mock CFHTLenS shear catalogues correspond to a ‘Hartlap factor’ of 0.86, and have the correct effective source density, shape noise, photometric redshift scatter, small-scale masking, and other properties that match the real survey. Given the different effective areas of the four CFHTLenS regions, the final covariance was obtained through area-weighted averaging of the patch-dependent inverse covariances. We refer the reader to Joudaki et al. (2017) for further details. The number of simulations used for the analysis of different sets of KL modes was scaled with the size of the data vector as described in the case of the Fourier-space measurements.

4 Results

4.1 The Karhunen-Loève transform

Before presenting the main result of this paper: i.e. the level of data compression that can be achieved through the KL transform, it is useful to review the key steps of the method, as described in Section 2.2, and how we apply them in practice to the data. Figure 8 shows the diagonal part of the angular $E$-mode power spectra, $C_{\ell}^{\alpha\alpha}$, as computed by CCL, together with the expected noise $N_{\ell}^{\alpha\alpha}$ as a function of the multipole $\ell$ for the best fit parameters of Joudaki et al. (2017), which we use here as our fiducial cosmology. The fact that the noise dominates over the signal for a wide range of scales, and that the off-diagonal signal terms (not shown here) have amplitudes that are similar to the diagonal (given the radially cumulative nature of weak lensing), indicates that different modes are tightly correlated and therefore that the effective number of signal-dominated independent radial modes should be small (i.e. the information carried out by the angular power spectrum is inefficiently spread out all over the redshift bins).

As described in Section 2.2, the two advantages of working with KL-transformed power spectra, is that ideally some combination of the redshift bins is constructed to have uncorrelated radial modes and that the information is compressed in the first few modes. This second feature can be seen in Figure 9, where we show the $\ell$-dependence of the eigenvalues of Eq. 9, which represent the KL-transformed power spectra (labelled $\Lambda^{\alpha}_{\ell}$ here). This quantity measures the amount of signal to noise present in each mode, where the KL-transformed noise power spectrum is simply unity ($\delta_{\alpha\beta}$) for all $\ell$s by construction. It is clear that the majority of the information is concentrated in the first two bins, while the information we lose by considering only the first three bins is almost negligible and the remaining bins are noise dominated.

A further simplification comes by assuming our KL transform to be $\ell$-independent, as proposed in Alonso (2018). This is necessary in real space, where $\ell$ multipoles are replaced by angles $\theta$888Note that the KL transform uses orthogonal eigenfunctions, and this is automatically achieved by decomposing the full sphere into spherical harmonics (or Fourier coefficients in the flat-sky approximation), and therefore the method lives naturally in $\ell$ space., but they can be safely used even in Fourier space. This is shown in Figure 10, where we plot the first three principal eigenvectors of Eq. 9 as a function of redshift. Each colour represents a single eigenvector, while different shades refer to different $\ell$s: lighter colours for low-$\ell$s and darker ones for high-$\ell$s. Excluding the first few $\ell$s (where cosmic variance dominates the uncertainties and little information can be recovered), all eigenvectors for the same $\alpha$ clearly have the same shape. Moreover, the blue lines, which represent the first eigenmode and dominate the information content, are the least scattered even at low-$\ell$. Then, an averaged version of the KL transform $e^{\alpha}_{p,\ell}$ can be written as

[TABLE]

where the $2\ell+1$ factors are included to roughly account for the statistical weight of each harmonic multipole. When performing the real-space analysis, we used this $\ell$-independent KL transform. In Fourier space we tested both options, and we will show that dropping the scale dependence of the eigenvectors makes little or no difference.

The measured power spectra for the first three KL eigenmodes using a scale-dependent KL transform are shown in Figure 11, together with the theoretical prediction corresponding to the fiducial cosmological model used in the analysis. This illustrates how the method concentrates the bulk of the signal-to-noise into the first few modes.

Finally, it is worth emphasizing that, since a critical use of data compression is to reduce the number of simulations needed to estimate the covariance matrix, when analyzing different compression schemes in the next section, we used a covariance matrix computed from a reduced set of simulations, scaled from the full set by the size of the data vector. This allows us to show that achieving a compression factor $f$ will allow us to carry out the analysis using $f$ times fewer simulations.

4.2 Observational constraints from the KL decomposition

In order to study the efficiency of different data compression schemes, we explore the posterior distribution of cosmological parameters using an MCMC sampler999In particular we implemented our likelihood in both emcee (Foreman-Mackey et al. 2013) and MontePython (Audren et al. 2013, Brinckmann and Lesgourgues 2018).. We assume that the two-point functions (both in real- and Fourier-space) follow a Gaussian likelihood. This is known to be formally incorrect (Sellentin and Heavens, 2018; Hamimeche and Lewis, 2008), and can give rise to parameter shifts below the $1\sigma$ level. For the purposes of this analysis (i.e. to show the performance of the KL transform in terms of data compression) this effect can be ignored, particularly given the fact that most of the constraining power is obtained from the smaller scales, where the large number of available modes makes the Gaussian approximation acceptable due to the central limit theorem. As in Joudaki et al. (2017), we use a standard $\Lambda$CDM model and vary five cosmological parameters: the Hubble parameter $h$, the physical baryon density parameter $\Omega_{b}h^{2}$, the physical cold dark matter density parameter $\Omega_{c}h^{2}$, the primordial comoving curvature power spectrum amplitude $A_{s}$ and the spectral index $n_{s}$. The priors used in this analysis are uniform and reported in Table 3. As investigated in Joudaki et al. (2017), the impact of priors is non-negligible in determining the allowed parameter space. However, since the scope of this work is just to show the goodness of the KL transform in reducing the data vector without losing any significant information, we choose to keep these priors fixed to those used in previous analyses for all our runs.

In Figure 12, we show the posterior distribution for the parameters $\Omega_{m}\equiv\Omega_{c}+\Omega_{b}$ and $S_{8}\equiv\sigma_{8}\sqrt{\Omega_{m}/0.3}$ obtained with the analysis in Fourier space and with an $\ell$-dependent KL-transform. The derived quantity $S_{8}$ is of particular interest in weak lensing analyses since cosmic shear has the most constraining power on this parameter and it has proven to be less dependent on the specific choice of priors. We show a few representative compression schemes: (i) full is the case where no KL-compression is used, (ii) in 7_diag we consider all $7$ radial KL modes but only using their auto-correlations, (iii) in 2_off we use the first 2 radial modes but we consider their auto- and cross-correlations, (iv) 2_diag is the same as 2_off but using only auto-correlations, and (v) 1 is the maximum compression case, where we consider only the first KL mode. Contours for the same parameters are shown in Figures 13 and 14 for the Fourier analysis with $\ell$-independent KL-transform and the real-space analysis respectively, where the latter also compares the real-space and Fourier-space results directly. Several conclusions can be drawn from these results:

• •

We are able to place similar constraints on the most relevant cosmological parameters in terms of final uncertainties in all cases.

• •

We observe small shifts in the maximum-likelihood parameter values for different compression schemes. These are well within the $1\sigma$ level, and are to be expected, since different compression schemes effectively use different sectors of the data vector. A further source of random shifts is the use of different sets of simulations to compute the covariance matrix for the different data compression schemes. We checked using different simulations to compute the covariance matrix that these shifts can be considered as statistical fluctuations.

• •

We obtain qualitatively equivalent results in terms of data compression efficiency using $\ell$-dependent and $\ell$-independent KL modes in a Fourier-space analysis, as well as using a real-space pipeline.

In order to get a quantitative sense of the compression efficiency of the KL transform, in Table 4 we show the Figure-of-Merit (FoM) obtained in the various cases, calculated as the inverse of the areas of the 1$\sigma$ contours in the 2D $\Omega_{m}$-$S_{8}$ plane. It is possible to notice a substantial difference between the real-space and the Fourier-space analyses. This is to be expected. We believe that most of these differences can be attributed to the different methods used to estimate the covariance matrices in both cases.101010It is interesting to note that both analyses also obtain substantially different values for the best-fit parameter $\chi^{2}$ with $\chi^{2}/{\rm ndof}=1.47$ for the real-space analysis and $\chi^{2}/{\rm ndof}=0.98$ for the Fourier-space one. In addition, it is never possible to match both sets of scale cuts, giving rise to different levels of constraining power. More interestingly, for the purposes of data compression, for a given analysis choice we observe only insignificant variations in the FoM as we reduce the number of KL modes left in the data vector. This clearly indicates that the information contained in the full case is essentially the same as in the most compressed one. Also, by comparing the results for the two different Fourier-space analyses, i.e. the $\ell$-dependent (SD) and the $\ell$-independent (SI) KL transform, we can conclude that there is no significant discrepancy between the two if we take into account cross-correlations between different radial modes. When neglecting the cross-correlations of the power spectra, we find a small ($\sim 10\%$) degradation in the FoM for the $\ell$-independent case. This indicates that, while considering an $\ell$-dependent KL transform may retain marginally more information especially when we maximally compress our data vector, the loss of information is small. It is interesting to notice that the FoM sometimes increases when we further compress the data vector. By construction, it is impossible to gain constraining power by compressing the data. This indicates that the apparent increasing of the FoM can be interpreted only as a statistical fluctuation.

In addition, since $S_{8}$ is the best constrained parameter with current weak lensing surveys, we can estimate the goodness of the KL compression by looking at its 1D posterior distribution. We compared the absolute values of the relative shift in the means ($\left|\Delta S_{8}/\sigma(S_{8})\right|$), and of the relative difference of the standard deviations ($\left|\Delta\sigma(S_{8})/\sigma(S_{8})\right|$) between the full and the 1 cases. These two have been chosen as representative of the cases with none and maximum data compression respectively. For the shift in the means we find $\sim 0.38$ for both the Fourier-space analyses (SD and SI) and $\sim 0.25$ for the real-space analysis. Moreover, the relative difference of the standard deviations is $\sim 0.04$ for the Fourier-space SD analysis, $\sim 0.08$ for the Fourier-space SI analysis and $\sim 0.01$ for the real-space analysis. This indicates not only that the amount of information lost with the KL-transform is practically negligible, but also that the shifts in the posterior distributions are compatible with expected statistical fluctuations caused by ignoring part of the signal.

In the best case, comparing all the data points used in the real-space analysis of Joudaki et al. (2017) with the maximum data compression of the Fourier-space pipeline, we are able to achieve a compression factor of $\simeq 30$, which greatly simplifies the problem of computing the covariance matrix from both mock catalogues and with analytic methods.

4.3 Relaxing the assumptions

The constraints obtained in the previous section show that even with the maximum KL compression the loss of information is negligible. However, they were obtained using a number of assumptions that can be relaxed in order to test the robustness of our conclusion.

As the KL transform has to be calculated assuming a fiducial model, the loss of information is minimal if the fiducial model is a good fit to the data. This is the main reason to choose the best fit of Joudaki et al. (2017) as our fiducial model. To explore the impact of the fiducial cosmology on the performance of this method, we have repeated our analysis using a fiducial model that is significantly far away from the best fit. In particular we use a flat $\Lambda$CDM model with $(h,\,\Omega_{b}h^{2},\,\Omega_{c}h^{2},\ln(10^{10}A_{s}),n_{s})=$ $(0.8,0.02,0.16,3.2,0.8)$, which yields a $\chi^{2}$ for the CFHTLenS data that is twice as large as that of the best fit. The result is shown in Figure 15, which shows the posterior distribution for the parameters $\Omega_{m}$ and $S_{8}$. We find contours that are very similar to the ones shown in Figure 13 and, as in the previous section, we find that the shift in the means and the relative difference of the 1-$\sigma$ errors for $S_{8}$ are $\sim 0.32$ and $\sim 0.06$ respectively. The performance of the KL transform is therefore only mildly affected by the choice of fiducial model.

Moreover, the results of our main analysis are in principle only valid for standard flat $\Lambda$CDM models described by the parameters in Table 3. An extended parameter space could in principle benefit from including additional KL modes in order to break parameter degeneracies. To explore this possibility, we have repeated our analysis for a “$w$CDM” model, in which we parametrize the evolution of the equation of state of dark energy as $w(a)=w_{0}+(1-a)w_{a}$ (Chevallier and Polarski, 2001; Linder, 2003), where $a$ is the scale factor, and $(w_{0},w_{a})$ are two new free parameters of the model111111We impose flat priors on these parameters with ranges $w_{0}\in[-3,0]$ and $w_{a}\in[-3,3]$, while leaving all other priors fixed to the values in Table 3.. The results are shown in Figure 16 both for the $\Omega_{M}$-$S_{8}$ plane (top panel) and the $w_{0}$-$w_{a}$ plane (bottom panel). As expected, the parameter contours for $\Omega_{M}$-$S_{8}$ are broader than in the $\Lambda$CDM case. The corresponding parameter shift and degradation in the final error on $S_{8}$ when comparing the full and 1 cases are $\Delta S_{8}/\sigma_{8}=0.58$, $\Delta\sigma(S_{8})/\sigma(S_{8})=0.12$. These are noticeably larger than in the $\Lambda$CDM case, reflecting the fact that additional information is encoded in the remaining KL modes. This is more evident in the $w_{0}$-$w_{a}$ plane, where the uncertainties on both parameters degrade noticeably if only a single mode is used. However, we can see that, by including both the first and second KL modes and considering their cross-correlation (dashed violet contours), we can recover the same constraints found when using the full data vector (which is also true for the $\Omega_{M}$-$S_{8}$ plane). Therefore, the compression efficiency degrades in this case from a factor $\sim 30$ to a factor $\sim 10$, which is still significant.

This result should come as no surprise, since figures 9 and 11 show that, even if the $S/N$ ratio is dominated by the first KL mode, some information can still be extracted from the second one. In the context of the $w$CDM model, this can be interpreted as follows: although the overall amplitude of the matter fluctuations can be well captured by a single, high signal-to-noise measurement of the weak lensing power spectrum, corresponding to the first KL mode, information about their growth as a function of time is vital in order to distinguish between different dark energy models. As shown by the red lines in Fig. 10, the second KL mode weights the high-redshift and low-redshift data differently, and is therefore able to partially recover the growth history. This was also shown in Alonso (2018), where it was found that an LSST-like experiment will require at least 3 KL modes in order to recover optimal constraints in an evolving dark energy scenario with massive neutrinos. It is also worth emphasizing the fact that, as shown by this exercise, even though the KL modes are uncorrelated in the fiducial cosmology, additional information can be obtained by exploring the parameter dependence of their cross-correlation.

This result shows that the performance of the KL decomposition depends upon the parameter space on which it is used, and therefore a preliminary forecast is required for any experiment before deciding upon a given compression scheme. Note that this is also a required step before defining the redshift binning to be used in a standard tomographic analysis. It is also worth noting that our particular implementation of the KL transform was aimed to maximize the overall $S/N$ of the data. The method is however more general than that, and can be tuned to optimize the final constraints on any parameter. Therefore, the results shown here could potentially be improved if we fine-tuned our KL decomposition for e.g. $w_{0}$ and/or $w_{a}$.

4.4 Comparison with other analyses

We compare here the performance of the method described in this paper to other similar analyses carried out on the CFHTLenS data.

The closest study of a radial decomposition of weak lensing data is that of Kitching et al. (2014), which uses spherical Bessel functions as a radial basis. This analysis makes use of 50 radial modes, compared to our 1 or 2 KL eigenmodes. The fact that their basis is orthogonal by construction also implies that the cross-correlations between different modes cannot be neglected (see e.g. their Figure 4). One advantage of the modes used by Kitching et al. (2014) is that, for accurate enough photometric redshifts, each of them can be associated with a physical wavenumber $k$, and therefore the effect of baryons and non-linearities in the matter power spectrum can be kept under control by using a conservative cut on $k$. It is also worth noting that Kitching et al. (2014) use the map-level modes as their data, instead of summarising them into power spectra. The total number of modes used at the likelihood level is therefore very large, but the use of map-level quantities simplifies the problem of computing 2-point function covariances.

Another interesting analysis is that of Asgari et al. (2017), in which the authors make use of a compressed set of Complete Orthogonal Sets of E$/$B-mode Integrals (COSEBIs) calculated from the CFHTLenS data. Their analysis proposes the use of a set of 20 COSEBIs for these data, which corresponds to a data reduction factor of $\sim 7$ with respect to their raw initial data vector. Although this is not checked explicitly in this analysis, the resulting compressed COSEBIs should preserve the constraining power on a basic set of 5 cosmological parameters.

Finally, Simon et al. (2015) have made use of a KL decomposition to reduce the dimensionality of third-order cosmic shear statistics in CFHTLenS . The compression factor achieved here is modest ($\sim 4$), but the use of this type of analyses for these datasets is extremely relevant, given the large number of possible three-point function combinations.

5 Discussion

The main goal of this paper has been to show the potential of the Karhunen-Loève transform in reducing the dimensionality of the data vector for cosmic shear measurements. We applied our method to weak lensing data from the CFHTLenS survey. The strong correlation between the weak lensing signal at different redshifts, makes cosmic shear particularly well suited for this type of compression. The fact that other 3D data decomposition frameworks have been applied to the CFHTLenS data in the past additionally allows us to make a comparison between different strategies.

The KL transform is designed to determine the linear combinations of the data that maximize the amount of information about a given parameter. Our main analysis of the angular two-point functions has been performed in Fourier space, both considering an $\ell$-dependent and an $\ell$-independent KL transform that maximizes the $S/N$ of the data. In this context, the resulting set of radial eigenfunctions can be likened to a Fourier decomposition along the redshift coordinate (i.e. the KL transform is a generalization of the standard 3D power spectrum analysis to weak lensing data). For comparison with previous literature, we have also applied our method to real-space correlation functions, and show how the compression works in the different cases. In Fourier space we used the CFHTLenS catalogue to generate masks and maps of the observed area of the sky. We then used these maps to estimate the weak lensing power spectra. For the real-space analysis, we directly used the correlation functions and covariance matrix constructed in Joudaki et al. (2017). In all cases, theoretical predictions for the power spectra and correlation functions given the cosmological parameters were obtained using CCL. In order to calculate the KL transform, our fiducial point was the best fit of Joudaki et al. (2017). Choosing cosmological parameters close to the best fit ensures a nearly optimal compression, while random parameters would potentially degrade the efficacy of the method. We then applied the KL transform to our two-point functions and explored the posterior parameter distribution using MCMC methods.

The main result of this paper is that, both in real and Fourier spaces, the KL compression works exceptionally well. Indeed, it is possible to compress the original data vector, with $280$ and $140$ elements in real and Fourier space respectively, down to a vector with $10$ and $5$ elements without any appreciable loss of information or parameter biases. In other words, we achieve a compression factor of $\sim 30$ for CFHTLenS, which simplifies the problem of computing covariance matrices from mock catalogues by the same factor, or by a factor $\sim 800$ using analytic methods. This is shown in Table 4, where only insignificant deviations in the FoM are found as we reduce the number of KL modes considered. Additionally, we observe that considering an $\ell$-independent KL transform does not significantly degrade the constraints, which simplifies the application of this method to real-space analyses.

As shown in Section 4.3, almost the same result holds after relaxing some of our assumptions. First, we have shown that the performance of the method is robust against errors in the fiducial cosmology chosen to define the KL eigenmodes. To show this we repeated our analysis for a fiducial cosmology that is substantially disfavoured by the CFHTLenS data, finding no significant difference in the final constraints. We also extended the parameter space by replacing the cosmological constant with a fluid with time dependent equation of state, using the standard $(w_{0},w_{a})$ parametrization. In this case, we find that, even though the final constraints are affected by the prior on these extended parameters, additional information can indeed be gained by including a second KL mode. This can be expected, given the non-zero signal contained in this mode, and can be interpreted as the second mode being able to recover growth information. The compression factor in this case degrades by a factor $3$ with respect to the $\Lambda$CDM analysis, although this could potentially be further improved by using a KL transform that optimizes the recovery of specific parameters (e.g. $w_{0}$ and/or $w_{a}$).

We must note a number of caveats in our analysis. Most prominently, we have not marginalized over possible photometric redshift systematics, including catastrophic outliers. Standard tomographic approaches have traditionally encapsulated these as additional nuisance parameters per redshift bin (e.g. overall shifts in the redshift distribution) that are marginalized over when evaluating the likelihood. This is not necessarily appropriate for the KL decomposition used here. The analogous approach would be to marginalize over the redshift distributions associated to each KL mode used ($\tilde{n}_{p}(z)=e_{p}^{\alpha}n_{\alpha}(z)$, using the notation in Section 2). Intrinsic alignments (Troxel and Ishak, 2015) is another systematic that may affect the performance of the method, given its non-local nature. Given the low significance of this effect, compared to the expected lensing signal, we expect it to have a negligible impact on the set of optimal eigenmodes, which are defined purely on the basis of signal-to-noise ratio. However, it could degrade significantly the compression efficiency in terms of final parameter constraints. We leave this study for future work. We have also ignored systematics associated to shear calibration, source clustering or baryonic effects on the matter power spectrum. Current and future imaging surveys draw most of their cosmological constraining power from the joint analysis of galaxy clustering and weak lensing (van Uitert et al. 2018, Joudaki et al. 2018, Abbott et al. 2018). It would therefore be interesting to explore the different ways in which the KL transform could be applied to this combined analysis. Although a direct application of the method would yield KL modes that mix clustering and shear, an alternative would be to combine the KL modes of both probes taken individually, in order to separate the different systematic uncertainties that affect them. We leave these investigations for future work.

Finally, it is worth noting that a number of alternative data compression methods have been proposed in the literature that can achieve significantly larger compression factors (e.g. Reichardt et al. 2001; Heavens et al. 2017; Alsing and Wandelt 2018; Alsing et al. 2018), down to a single data point per free parameter of the model. Although some aspects of these methods depend on having a sufficiently large sample of simulations, they are some of the most promising avenues to alleviate the complexity of computing covariance matrices (or even to forgo likelihood computations altogether). One possible advantage of the KL transform as presented here with respect to these methods is the availability of intermediate data products (KL mode maps and power spectra) that can be inspected to identify systematics in the analysis pipeline or the theory model. Data compression techniques of any form will be of tremendous use in future surveys such as LSST or Euclid, simplifying computationally demanding tasks, such as the computation of large covariance matrices.

Acknowledgements

We thank Chris Blake for the use of CFHTLenS mock catalogues in the creation of the covariance matrix of the shear correlation functions. We also thank Chris Blake, Pedro Ferreira, Alan Heavens, Catherine Heymans, Lance Miller and Licia Verde for useful discussions. EB and SJ are supported by ERC H2020 693024 GravityLS project, the Beecroft Trust and the Science and Technology Facilities Council (STFC). DA acknowledges support from the Beecroft Trust and from STFC through an Ernest Rutherford Fellowship, grant reference ST/P004474/1. This work is based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/DAPNIA, at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l’Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency. CFHTLenS data processing was made possible thanks to significant computing support from the NSERC Research Tools and Instruments grant program.