# Parameter inference for weak lensing using Gaussian Processes and MOPED

Arrykrishna Mootoovaloo, Alan Heavens, Andrew Jaffe, Florent Leclercq

► To cite this version:

Arrykrishna Mootoovaloo, Alan Heavens, Andrew Jaffe, Florent Leclercq. Parameter inference for weak lensing using Gaussian Processes and MOPED. Monthly Notices of the Royal Astronomical Society, 2020, 497 (2), pp.2213-2226. 10.1093/mnras/staa2102. hal-05347496

HAL Id: hal-05347496

https://hal.science/hal-05347496v1

Submitted on 6 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Parameter inference for weak lensing using Gaussian Processes and MOPED

Arrykrishna Mootoovaloo,*, Alan F. Heavens, Andrew H. Jaffe and Florent Leclercq

ICIC, Astrophysics, Imperial College, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK

Accepted 2020 July 10. Received 2020 July 3; in original form 2020 May 19

## ABSTRACT

In this paper, we propose a Gaussian Process (GP) emulator for the calculation both of tomographic weak lensing band powers, and of coefficients of summary data massively compressed with the MOPED algorithm. In the former case cosmological parameter inference is accelerated by a factor of ~10–30 compared with Boltzmann solver *CLASS* applied to KiDS-450 weak lensing data. Much larger gains of order 10³ will come with future data, and MOPED with GPs will be fast enough to permit the Limber approximation to be dropped, with acceleration in this case of ~10⁵. A potential advantage of GPs is that an error on the emulated function can be computed and this uncertainty incorporated into the likelihood. However, it is known that the GP error can be unreliable when applied to deterministic functions, and we find, using the Kullback–Leibler divergence between the emulator and *CLASS* likelihoods, and from the uncertainties on the parameters, that agreement is better when the GP uncertainty is not used. In future, weak lensing surveys such as Euclid, and the Legacy Survey of Space and Time, will have up to ~10⁴ summary statistics, and inference will be correspondingly more challenging. However, since the speed of MOPED is determined not the number of summary data, but by the number of parameters, MOPED analysis scales almost perfectly, provided that a fast way to compute the theoretical MOPED coefficients is available. The GP provides such a fast mechanism.

**Key words:** gravitational lensing: weak–methods: data analysis–methods: statistical–cosmological parameters–large-scale structure of Universe.

## 1 INTRODUCTION

With the continuous advancement of technology, our understanding of the Universe has progressively improved, starting with the radiation of the cosmic microwave background (CMB) from experiments such as *COBE* (Smoot et al. 1992; Jaffe et al. 2001), *WMAP* (Spergel et al. 2003), and *Planck* (Planck Collaboration et al. 2014, 2016) to ongoing observations of supernovae as standard candles (Betoule et al. 2014) and large-scale structure probes (Anderson et al. 2014). While these experiments place the Λ-cold dark matter (ΛCDM) model on a firm footing, there is still the need to understand better dark matter and the evolution of dark energy. A powerful probe of the geometry of the Universe is cosmic shear, the weak lensing effect observed as a result of bending of light between the observer and background galaxies due to intervening large-scale structure. Cosmic shear is proving to be a powerful method for measuring the temporal and spatial properties of dark species statistically over a large sample of sources (Köhlinger et al. 2017).

However, one possible bottleneck in the era of massive surveys lies in the forward computation of two-point statistics from cosmological parameters, either correlation functions or power spectra (Kilbinger 2015). The complexity of the problem is exacerbated when performing weak lensing analysis in *n* tomographic redshift slices, for which *n*(n + 1)/2 auto- and cross-correlations are needed for each angular scale probed. If intrinsic alignments (Heavens, Refregier & Heymans

2000a; Hirata & Seljak 2004) are included, we require an additional calculation of *n*(*n* + 1) power spectra. In general, running this full forward model is computationally expensive. The problem is even worse in very large simulation settings. As an example, Heitmann et al. (2009) argued that a naive analysis to obtain cosmological constraints from future surveys such as Legacy Survey of Space and Time (LSST, Seppala 2002; LSST Science Collaboration et al. 2017) from weak lensing shear spectra will take a 2048 processor machine 20 yr. An alternative approach is to replace the simulator (the full forward model) by an approximate mathematical function, often referred to as a metamodel, surrogate model, or emulator.

Various types of emulators have been designed for different purposes in cosmology. As an example, Parameters for the Impatient Cosmologist (PICO) (Fendt & Wandelt 2007b) was developed to accelerate parameter estimation for the CMB. The underlying idea is to use an order-p polynomial and a clustering method to interpolate between power spectra generated at specific points in parameter space. On the other hand, neural networks have also been used as emulators for generating power spectra. Auld et al. (2007) used a three-layer multilayer perceptron (MLP) algorithm to construct an emulator to accelerate the calculation of the CMB power spectra. Neural networks learn the non-linearities between the input space and output space by optimizing the weights associated with each neuron. In the same spirit, Schmit & Pritchard (2018) developed a two-layer MLP algorithm to emulate the 21-cm power spectrum. Similarly, Agarwal et al. (2012, 2014) designed a neural network emulator for non-linear matter power spectrum interpolation based on six cosmological parameters. Charnock, Lavaux & Wandelt (2018) also

<footnote>* E-mail: a.mootoovaloo17@imperial.ac.uk</footnote>
used neural networks to find functionals of the data that maximize the Fisher information matrix. Manrique-Yus & Sellentin (2020) devised a neural network scheme, coupled to a Markov Chain Monte Carlo (MCMC) algorithm, to accelerate cosmological parameter inference. In an alternative approach, Habib et al. (2007) developed a statistical framework using Gaussian Processes (GPs) to emulate matter power spectrum, and Schneider, Holm & Knox (2011) used a similar method to emulate CMB temperature power spectrum. An analogous approach was used and extended in the Coyote Universe collaboration (Heitmann et al. 2009, 2010; Lawrence et al. 2010; Heitmann et al. 2014). Recently, Bird et al. (2019) and Rogers et al. (2019) developed a GP emulator for Lyman-$\alpha$ forest simulation. Kern et al. (2017) designed a GP emulator for the 21-cm power spectrum.

An alternative to predicting the signal using machine learning techniques is to predict the likelihood. Instead of building emulators at the level of the theory/data, another option is to construct a likelihood regressor. Fendt & Wandelt (2007a) extended the PICO formalism to fit a likelihood function. Leclercq (2018) developed the Bayesian Optimisation for Likelihood-Free Inference (BOLFI) algorithm (Gutmann & Corander 2015) to construct a likelihood regressor which fits two cosmological parameters $w_0$ and $\Omega_m$ for the Joint Light-curve Analysis (JLA) data set (Betoule et al. 2014). In particular, in the BOLFI approach, the uncertainty from the GP appears in the acquisition function, which is used to choose the next design point where the simulation will be run. However, the surrogate model itself is not a perfect reconstruction of the true likelihood function. Instead, if the emulator were built at the level of the theoretical model, the GP uncertainty can be propagated in the full forward model. In the same spectrum, Alsing, Wandelt & Feeney (2018) developed a density-estimation likelihood-free method and argued that it requires less tuning compared to traditional approaches such as Approximate Bayesian Computation (ABC). Leclercq et al. (2019) also developed a likelihood-free approach to infer power spectrum and cosmological parameters from forward simulations only.

Each type of emulator (polynomial, neural network, and GPs) has its own advantages and disadvantages. In the case of polynomial regression, one has to specify the basis functions prior to performing interpolation. However, it provides a full predictive distribution of the function at a test point. On the other hand, neural network regression requires much empirical work, for example, choosing an appropriate optimiser, setting the learning rate of the optimiser and choosing the number of layers and neurons. The training stage is computationally very challenging if we have tens or hundreds of functions to interpolate, as in our case, but it can be done (Manrique-Yus & Sellentin 2020). The outputs from a neural network are point estimates, unless we consider Bayesian neural networks which attempt to get a full posterior distribution on the weights of the neural network using variational inference methods, but this is also challenging (Kendall & Gal 2017).

In this paper, we develop a GP emulator for speeding cosmological parameter inference for weak lensing cosmology whilst keeping the number of calls to the full solver to the order of a few tens of thousands only, to avoid an expensive and potentially unstable Cholesky decomposition. In other contexts, where the generation of training point data may be done through expensive full numerical simulations, the number may need to be further restricted. Despite this restriction, an important and interesting result is that the posterior densities between the emulator and the full solver are comparable. On the other hand, in the case of full N-body simulations, one might be restricted to only a few hundred forward simulations,

and the mechanism (MOPED + GP) proposed in this work (see discussion below) will enable us to develop promising inference schemes for this scenario. A useful advantage of GPs is that we also get a theoretical predictive uncertainty. However, note that the latter depends on the GP model, including the kernel. Building such a framework is a challenging task. It is well known that non-parametric regression methods such as GPs suffer from the so-called curse of dimensionality (Geenens et al. 2011) as a result of the training data being sparse in high-dimensional spaces. The only way to obtain accurate interpolation is to add information to our model and this can be achieved by adding more and more training points. However, this involves a penalty, not only in terms of computations but also storage. Moreover, from the cosmology perspective, drawing inferences from current weak lensing data is a demanding task, for we have to correctly account for sources of systematic error which add new effective parameters to our model. In addition to this, given the current status of the weak lensing field with relatively few data points at low signal-to-noise ratio, cosmological parameters are not very well-constrained compared to (e.g.) the CMB. As a result, our training points have to be distributed over a larger volume in parameter space.

Fortunately, sophisticated techniques such as Latin Hypercube (LH) sampling (McKay, Beckman & Conover 1979), with the appropriate transformations at the level of both the input covariates and the response, can simplify the problem. Re-casting the full problem as a hierarchical model enables us to leverage the predictive uncertainty of the reconstructed probabilistic functions, thus propagating the uncertainty consistently in the statistical framework. This is in principle a major advantage of GPs, but the GP error is not always a reliable measure when applied to deterministic functions (Karvonen et al. 2020; Wang 2020), and we find that the posteriors are more accurate if it is not included (see Section 6 for further details).

A powerful application of the GP emulator is to combine it with the extreme data compression algorithm MOPED (Heavens, Jimenez & Lahav 2000b) to further speed up parameter inference. MOPED is essentially a lossless data compression technique which, irrespective of the size of the original data, compresses the latter to the number of parameters in the model. One can simply replace the theoretical prediction by the surrogate model in the MOPED-likelihood function and use the MOPED vectors to compute the MOPED coefficients. However, in order to achieve the full acceleration that MOPED allows, a fast mechanism to generate the expected coefficients themselves is needed, and GPs can do this.

In this study, as a first application, the GP emulator is built first for tomographic weak lensing band power spectra, and secondly for MOPED coefficients. With a well-designed emulator, it is possible to obtain reliable marginalized posterior densities of cosmological and nuisance parameters with a few hundred forward simulations only. We find that for a small number of training points, the emulator yields posterior distributions close to the true posterior whilst still maintaining a decent speed for each likelihood evaluation. When combined with MOPED data compression, the GP allows for very large acceleration of analysis of large data sets, which may otherwise be prohibitively expensive.

The paper is organized as follows: in Section 2, we provide a brief overview on weak lensing theory, and in Section 3, we touch upon GPs before systematically going through the main steps taken to build the emulator in Section 4. In Section 5, we discuss how the MOPED algorithm can be used alongside the emulator and we present our results in Section 6 before concluding in Section 7. Throughout, we will assume a flat universe.
## 2 WEAK LENSING

Gravitational lensing is the bending of light as it propagates through the inhomogeneous Universe, which leads to a coherent distortion of galaxy images. In particular, weak lensing can only be studied in a statistical sense since the distortion is small and requires averaging over a large sample of galaxies as a result of broad distribution of intrinsic ellipticities. Weak lensing has the advantage of probing the distribution of matter and not of biased tracers such as galaxies which is hard to predict. In essence, we would ideally want to perform a full 3D statistical analysis but much of the cosmological information can still be retained by tomography, where objects are separated by redshift. We refer the reader to common literature on weak gravitational lensing (Bartelmann & Schneider 2001; Kilbinger 2015) for further details on these techniques.

In a weak lensing analysis, observables include galaxy positions, photometric redshifts, and shapes. The latter is given in terms of the ellipticity components, $\epsilon_1$ and $\epsilon_2$. In particular, this observed ellipticity, $\epsilon = \epsilon_1 + i\epsilon_2$ is related to the unlensed intrinsic ellipticity, $\epsilon_0$ via the shear field such that

$$ \epsilon \approx \epsilon_0 + \gamma. \quad (1) $$

The observed ellipticities are binned into pixels and redshift bin, $i$. An estimate of the shear field, $\hat{\mathbf{r}}_i$ is obtained by averaging the ellipticities in each pixel. This (complex) shear field can be expanded in terms of spin-weighted spherical harmonics, $sY_{\ell m}$ (Hu 2000; Castro, Heavens & Kitching 2005)

$$ \gamma_1(\mathbf{r}) \pm i\gamma_2(\mathbf{r}) = \sqrt{\frac{2}{\pi}} \sum_{\ell m} \int dk \, k^2 \gamma_{\ell m}(k)_{\pm 2} Y_{\ell m}(\hat{\mathbf{n}}) j_{\ell}(kr), \quad (2) $$

where $j_\ell$ is a spherical Bessel function, $k$ is a radial wavenumber, and $\ell$ is a positive integer while $m = -\ell, \dots, \ell$. The coefficients $\gamma_{\ell m}$ are related to the transform of the lensing potential, $\phi(\mathbf{r})$ by

$$ \gamma_{\ell m} = \frac{1}{2} \sqrt{\frac{(\ell + 2)!}{(\ell - 2)!}} \phi_{\ell m}(k). \quad (3) $$

Similarly, the expansion coefficients for the convergence field, $\kappa$, is

$$ \kappa_{\ell m} = -\frac{\ell(\ell + 1)}{2} \phi_{\ell m}(k). \quad (4) $$

The shear field can be decomposed into E and B modes corresponding to the curl- and divergence-free components. In particular, the convergence field, $\kappa^E$ contains most of the cosmological information since $\kappa^B$ is negligible in the absence of systematics (Castro et al. 2005). Under this condition, the E-mode lensing power spectrum between tomographic bins $i$ and $j$ is equal to the convergence power spectrum, that is, $C_{\ell, ij}^{EE} = C_{\ell, ij}^{KK}$ and is given, in the Limber approximation (Limber 1953; Loverde & Afshordi 2008) by

$$ C_{\ell, ij}^{EE} = \int_{0}^{\chi_H} d\chi \frac{w_i(\chi)w_j(\chi)}{\chi^2} P_{\delta} \left( k = \frac{\ell + 1/2}{\chi} ; \chi \right), \quad (5) $$

where $\chi$ is the comoving radial distance and $\chi_H$ is the comoving distance to the horizon. Without the Limber approximation, the integrals can be slow to compute, although faster methods are being developed (Fang et al. 2019). Crucially, the tomographic convergence power spectrum is sensitive to the background geometry and the growth of structure. It depends on the the 3D matter power spectrum, $P_\delta(k; \chi)$ which is a function of redshift (Weinberg et al. 2013). The weight function $w_i$ for a flat universe is

$$ w_i(\chi) = \frac{3\Omega_m H_0^2}{2c^2} \chi(1+z) \int_{\chi}^{\chi_H} d\chi' n_i(\chi') \left( \frac{\chi' - \chi}{\chi'} \right), \quad (6) $$

which depends on the lensing kernel. $\Omega_m$ is the present matter density, $H_0$ is the Hubble constant, and $c$ is the speed of light. An important quantity is the redshift distribution, $n_i(z)dz = n_i(\chi)d\chi$ which is normalized such that

$$ \int n_i(\chi) d\chi = 1. \quad (7) $$

For a weak lensing survey, the data vector consists of the measured shear per pixel for each redshift bin. At this point, in order to extract the shear power spectrum, one can either take a quadratic estimator approach using maximum-likelihood technique (Bond, Jaffe & Knox 1998) or employ a pseudo-$C(\ell)$ approach (Hinshaw et al. 2007). Alternatively, one can also build a full Bayesian hierarchical model, to infer the full shear power spectrum (Alsing et al. 2016; Alsing, Heavens & Jaffe 2017). Here, we focus on the tomographic band power spectra, as determined by Köhlinger et al. (2017).

## 2.1 Astrophysical systematics

Coupled to the E-mode power spectrum are various systematics which we should consider. For example, baryon feedback results in altering the power in high $k$. Although feedback is not fully understood, it is often parametrized through the bias function, $b^2(k, z)$ such that the modified power spectrum is

$$ P_{\delta}^{\text{mod}}(k, z) = b^2(k, z) P_{\delta}(k, z). \quad (8) $$

As an example, for the KiDS-450 analysis, the following fitting formula from van Daalen et al. (2011) was used

$$ b^2(k, z) = 1 - A_{\text{bary}} \left[ A_z e^{(B_z x - C_z)^3} - D_z x e^{E_z x} \right], \quad (9) $$

where $x = \log_{10}(k/1 \text{ Mpc}^{-1})$ and the other parameters $A_z$, $B_z$, $C_z$, $D_z$, and $E_z$ depend on the scale factor $a$. Moreover, we must account for intrinsic alignment effects which give rise to a preferred ellipticity orientation. The total lensing power spectrum between two redshift slices is a linear combination of the gravitational lensing (EE), intrinsic alignment (II), and interference (GI) power spectra. Specifically, the II effect is due to correlation of ellipticities in the local environment and contributes positively towards the total lensing spectrum. The second effect, GI, is due to correlation between tidally stretched foreground galaxies and the shear of background galaxies. The GI term subtracts from the total lensing spectrum. We model the power spectrum, following Köhlinger et al. (2017), as

$$ C_{\ell, ij}^{\text{tot}} = C_{\ell, ij}^{\text{EE}} + A_{\text{IA}}^2 C_{\ell, ij}^{\text{II}} - A_{\text{IA}} C_{\ell, ij}^{\text{GI}}, \quad (10) $$

where the II power spectrum, $C_{\ell, ij}^{\text{II}}$ and the GI power spectrum, $C_{\ell, ij}^{\text{GI}}$ respectively are

$$ C_{\ell, ij}^{\text{II}} = \int_{0}^{\chi_H} d\chi \frac{w_i(\chi)w_j(\chi)}{\chi^2} P_{\delta} \left( k = \frac{\ell + 0.5}{\chi} ; \chi \right) F^2(\chi), \quad (11) $$

$$ C_{\ell, ij}^{\text{GI}} = \int_{0}^{\chi_H} d\chi \frac{w_i(\chi)n_j(\chi) + w_j(\chi)n_i(\chi)}{\chi^2} P_{\delta} \left( k = \frac{\ell + 0.5}{\chi} ; \chi \right) F(\chi), \quad (12) $$

$$ F(\chi) = C_1 \rho_{\text{crit}} \frac{\Omega_m}{D_+(\chi)}, \quad (13) $$

and $A_{\text{IA}}$ is a free parameter to be inferred during sampling. This allows for the flexibility of rescaling the otherwise fixed normalization value, $C_1 = 5 \times 10^{-14} h^{-2} M_{\odot}^{-1} \text{Mpc}^3$. $\rho_{\text{crit}}$ is the critical density of the Universe, while $D_+(\chi)$ refers to the linear growth factor normalized to unity today.
## 3 GAUSSIAN PROCESSES

In this section, we provide a general outline of the basic concepts behind the GP, which is at the core of our cosmological parameter inference scheme. We refer the reader to Rasmussen & Williams (2006) for further details.

Our goal is to design an emulator, which takes as input a set of cosmological parameters, $\theta$, and outputs summary statistics, $y$, for example, power spectrum or band powers. We consider our regression problem to be in the form,

$$y = f + \epsilon, \quad (14)$$

where $f$ is a multivariate Gaussian when evaluated at an arbitrary number of points, and $\epsilon$ is an error. The distribution of $f$ is informed by a set of $N_{train}$ training points in parameter space, and has a Gaussian prior distribution, $f \sim N(0, k_{pq})$ with a covariance $k_{pq} \equiv k(\theta_p, \theta_q)$.

The kernel function $k$ which encapsulates the correlation between points in the parameter space is used to model the smoothness of the function. We choose the automatic relevance determination (ARD) kernel, also referred to as the squared exponential kernel, radial basis function or simply the Gaussian kernel:

$$k(\theta_p, \theta_q) = A^2 \exp \left[ -\frac{1}{2} (\theta_p - \theta_q)^T \Lambda^{-1} (\theta_p - \theta_q) \right], \quad (15)$$

where $\Lambda = \text{diag}(\lambda_1, \lambda_2, \dots, \lambda_d)$. A and $\lambda$ are referred to as the amplitude and characteristic squared length-scale respectively. In particular, the former determines the average distance of the function from the mean while $\lambda$ controls the smoothness of the function. The set of hyperparameters for this particular kernel is $\{A, \lambda\}$. This kernel has the nice property that it is fully differentiable and positive definite.

Unlike parametric regression where we define priors over parameters, we will now define a prior covariance over the functions directly and using Bayes' theorem, the joint posterior distribution of the functions is

$$p(f | \theta, y) = \frac{p(y | \theta, f) p(f | \theta)}{p(y | \theta)}, \quad (16)$$

where the likelihood $p(y | \theta, f) \sim N(f, \Sigma)$ and the prior is $p(f | \theta) \sim N(0, K)$. $\Sigma$ is the noise covariance covariance matrix. $K = k(\theta, \theta')$ is the kernel matrix by calculating equation (15) for every pair of points in $\theta$. Prior to making predictive inference, the GP is trained by finding the set of hyperparameters $\{A, \lambda\}$ which maximizes the Bayesian Evidence (marginal likelihood)

$$\log p(y | \theta) = -\frac{1}{2} y^T K_y^{-1} y - \frac{1}{2} \log |K_y| + \text{constant}, \quad (17)$$

where $K_y = K + \Sigma$. The first term in the marginal likelihood controls the fit to the data while the second term controls the model complexity. For numerical stability, we first compute the Cholesky factor, $L$, of $K_y = LL^T$, solve for $u$ in the linear system $Lu = y$ followed by solving for $\alpha$ in $L^T\alpha = u$. The marginal likelihood is then given by

$$\log p(y | \theta) = -\frac{1}{2} y^T \alpha - \sum_i \log L_{ii} + \text{constant}. \quad (18)$$

Moreover, the partial derivatives of equation (17) with respect to the kernel hyperparameters, $\eta = \{A, \lambda\}$ can be computed in closed form

$$\frac{\partial}{\partial \eta_i} \log p(y | \theta) = \frac{1}{2} \text{tr} \left[ (\alpha \alpha^T - K^{-1}) \frac{\partial K}{\partial \eta_i} \right] \quad (19)$$

and $\alpha = K_y^{-1} y$. The gradients are useful when maximizing the marginal likelihood when using gradient-based optimization. An-

<figure><img src="image_3.png" /><figcaption>**Figure 1.** A diagrammatic form of the core principle in this work. We substitute the most expensive part of the pipeline by surrogate models (GPs) built at the level of the band powers. The other blocks in the inference procedure, for example, for the computations related to the nuisance parameters, are unaltered.</figcaption></figure>

other option, which we will not use for this particular application, is to marginalize over the kernel hyperparameters (given an appropriate set of priors) in a fully Bayesian formalism. The reason for not taking this route is that the number of latent variables in the data model becomes large ($\sim 10^2$) and recall that the marginal likelihood for a GP is an expensive calculation.

For a given test point $\theta_*$, the posterior distribution of the function, $f_* = f(\theta_*)$ is a Gaussian distribution with mean and variance given by

$$\bar{f}_* = k_*^T \alpha \\ \text{var}(f_*) = k_{**} - k_*^T K_y^{-1} k_* \quad (20)$$

The GP approach is quite appealing as it predicts both the mean and variance, at the expense of an $O(N^2)$ operation, of the function for a particular test point. Moreover, as seen from equation (20), the mean of a GP is a linear predictor. Once $\alpha$ is calculated, the mean function can be quickly and easily calculated for any set of test points since it involves only $O(N)$ operations. In addition, if required, the analytical gradient with respect to the mean function at a particular test point can also be calculated. However, a GP can be prohibitively expensive for large data sets since it involves the computation of the Cholesky factor which has computational complexity of $O(N^3)$ during the training phase and requires $O(N^2)$ operations for the predictive variance. Besides, a GP has $O(N^2)$ memory requirements for storing the Cholesky factor (if the computation of the predictive error is required) and the vector $\alpha$.

In the next section, we look into how a GP emulator can be useful for future weak lensing surveys, for which we naively expect around 10 tomographic redshift bins and thousands of summary statistics (Euclid Collaboration et al. 2019). In particular, one can either emulate the band powers directly or the MOPED coefficients. We investigate both and discuss the advantages of using the MOPED coefficients in the following sections.

## 4 EMULATOR

In this section, we use the formalism presented above to build the emulator. In brief, the latter involves four main stages, (1) generating a set of design points, (2) running the full forward simulator at these points, (3) training the emulator, and (4) making predictions at test locations in the parameter space. Once this is done, the emulator is connected to an MCMC sampler to obtain the marginalized posterior distributions of the parameters in our model. A simple flow of the core idea is shown in Fig. 1. In the following, we touch briefly on the data we have used for our analysis before systematically going through the steps we have taken to build the emulator.
## 4.1 Data

We use the publicly available weak lensing data from Köhlinger et al. (2017) to test the performance of our emulator. We use three tomographic redshift bins, namely, $0.10 < z < 0.30$, $0.30 < z < 0.60$, and $0.60 < z < 0.90$ and the convergence power spectrum is computed in the range $10 < ℓ < 4000$. Moreover, we follow Köhlinger et al. (2017) and drop the first, second-to-last and last band powers in our analysis, that is, we use only the band powers corresponding to the following $ℓ$-ranges: $76 ≤ ℓ < 220$, $221 ≤ ℓ < 420$, $421 ≤ ℓ < 670$, and $671 ≤ ℓ < 1310$. For a three-bin tomographic analysis, we have six auto- and cross-tomographic power spectra to calculate. The data and covariance matrix for this problem are shown in Figs 2 and 3, respectively.

The emulator can be built at the level of the power spectra or the band powers. Here we choose to build a GP for each band power, giving 24 GPs. Alternatively, for likelihood-free inference methods, one can also emulate the likelihood directly using the GPs (see Fendt & Wandelt 2007a; Leclercq 2018). For power spectrum reconstruction, one can use the PICO method or an alternative, but constrictive, stance is to adopt the approach taken by Habib et al. (2007) to first learn a set of basis functions via singular value decomposition and model the resulting weights by a GP. However, building an emulator for weak lensing analysis needs to account for systematic effects, but some of these can be included analytically without emulation, resulting in an 8D GP, rather than 12 (six cosmological and six systematic parameters) if we were to emulate the likelihood.

## 4.2 Training points

The generation of the training points is a key ingredient for the emulator to perform well. Accurate high-dimensional regression is not easy, mainly due to the curse of dimensionality. With the formalism presented in this work and depending on the complexity of the function, one can reconstruct the function precisely and accurately in low dimensions, hence leading to an accurate likelihood as would be the case if we were to use the full simulator, CLASS (Lesgourgues 2011) in this case. As the dimensionality of the problem increases, we need an exponentially increasing number of training points to emulate the true function accurately.

In PICO, the training points were generated uniformly from a box whose sides were centred on the mean of a converged MCMC chain (consisting of ~60 000 cosmological models) and width $3\sigma$ along each direction. In the second release of PICO, they selected training points which lie within 25 log-likelihoods of the WMAP peak (Fendt & Wandelt 2007a). On the other hand, Auld et al. (2007) first drew 2000 training points from the same box defined in PICO and also added an extra 5000 training points drawn from a Gaussian distribution, whose covariance was twice the expected covariance matrix, centred on the maximum likelihood. These techniques perform quite well for two reasons: (1) by restricting the prior volume of the training points to the high likelihood regions allows the sampler explicitly to explore this specific region in parameter space, and (2) creating a data set with thousands of training points will also improve any regression method. A shortcoming of using these approaches is that the algorithm will not perform well in regions where there is no training point nearby (see appendix A in Habib et al. 2007 for a comparison of their method with PICO). This is a typical manifestation of almost any machine learning algorithm. They are good at making reliable predictions within a pre-defined prior, provided they are trained with enough data points. Building

machine learning algorithms in the small data regime is still in its infancy, hence an active area of research (Barz & Denzler 2019).

Moreover, if the training points are naively generated randomly from our pre-defined priors, we might not obtain a suitable coverage of the parameter space. A possible solution to this, is to use a grid but then the number of training points grows exponentially as the dimensionality of the problem increases. As an example, say, we have a 7D problem and we choose to have 10 points per parameter, then our training set will have 10 million points.

Alternatively, we can use LH sampling (McKay et al. 1979) which is a method for generating random samples from a multidimensional distribution in a controlled (quasi-random) way. A point is assigned such that it uniquely occupies its row and column respectively. This procedure generalizes to higher dimensional designs. In Fig. 4, we show the projection of 5 LH samples, which have been generated from a box in 8D and scaled by the pre-defined priors in Section 4.3, in 2D. In particular, we show the projection for three parameters only but the same applies for the other parameters, where each point uniquely occupies its corresponding row and its column. The LH method is now a ubiquitous tool for performing emulation in large simulation scenarios (Habib et al. 2007; Schneider et al. 2011; Schmit & Pritchard 2018) and is seen to be quite efficient, not only in producing a fair interpolation, but also provides reasonable posterior densities.

In this work, we adopt the LH approach to generate our training set. The LH samples are generated using the maximinLHS function from the lhesR package (Carnell 2012). This particular design relies on distance criterion (Johnson, Moore & Ylvisaker 1990) and the final design is a result of maximizing the minimum distance between points.

## 4.3 Priors

In our baseline emulator, we generate 1000 LH samples from a box, between 0 and 1. We first linearly transform these samples to the range of the pre-defined prior box for the six cosmological and two systematics parameters,

$$ \Omega_{cdm}h^2 \sim U[0.01, 0.40] $$

$$ \Omega_b h^2 \sim U[0.019, 0.026] $$

$$ \ln(10^{10} A_s) \sim U[1.70, 5.00] $$

$$ n_s \sim U[0.70, 1.30] $$

$$ h \sim U[0.64, 0.82] $$

$$ A_{\text{bary}} \sim U[0.0, 2.0] $$

$$ A_{\text{IA}} \sim U[-6.0, 6.0] $$

$$ \Sigma m_v \sim U[0.06, 1.0] $$

followed by running the full simulator at these points to obtain the total band powers. $U[a, b]$ denotes a uniform distribution with lower and upper limits $a$ and $b$ respectively. We apply a more restrictive prior than the original KiDS-450 prior [0.01,0.99] for $\Omega_{cdm}h^2$ since otherwise a large fraction of the LH samples we generate lie outside the region of parameter space constrained by the current weak lensing analysis. Moreover, having a smaller volume of parameter space also improves the performance of the emulator. The prior for the $A_{\text{bary}}$ is set to an upper limit of 2 (instead of 10 in Köhlinger et al. (2017)) because we found that, values of $A_{\text{bary}} \gtrsim 3$ lead to negative $b^2$, which implies an unphysical negative power spectrum. In the same spirit, large values of $A_{\text{bary}}$ lead to negative auto-correlated band powers and in some cases, the band power matrix (equation 21) was not positive definite. We also found that large values of neutrino masses, $\Sigma m_v \gtrsim$
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The E-mode band powers (data) used in our inference scheme, similar to the KiDS-450 analysis (Köhlinger et al. 2017). The $ℓ$-ranges are as follows: 76 ≤ $ℓ$ < 220, 221 ≤ $ℓ$ < 420, 421 ≤ $ℓ$ < 670, and 671 ≤ $ℓ$ < 1310. In particular, the auto-correlation band powers are along the main diagonal ($z1$ × $z1$, $z2$ × $z2$, and $z3$ × $z3$) for the three redshift bins 0.10 < $z_1$ ≤ 0.30, 0.30 < $z_2$ ≤ 0.60, and 0.60 < $z_3$ ≤ 0.90. The off-diagonal blocks show the unique cross-correlation band powers. The blue shaded regions indicate the $1\sigma$ level errors from the covariance matrix.</figcaption></figure>

1 eV result in almost half of the CLASS band powers in our training set to be nan. We therefore set an upper limit for $\Sigma m_v$ to 1 eV.

## 4.4 Transformations

Training the GPs with the LH samples from above might be sub-optimal, the reason being that the volume occupied by a hypercube grows exponentially with increasing dimensions. On the other hand, a sphered training set (hypersphere) has a smaller volume compared to its corresponding hypercube but with the same scaling with dimension. This transformation step is analogous to the one used by Fendt & Wandelt (2007b). Schneider et al. (2011) assessed in detail the effect of various transformations prior to building an emulator for the CMB power spectrum. They found that de-correlating the input space leads to significant improvements compared to working with the original form of the input parameters. The interpolation can further be improved if one uses a known Fisher information matrix specific to the problem.

The transformation matrix can be calculated as follows: we first compute the sample covariance, $C_θ$ of the 1000 input parameters, $θ$ to the emulator (see Table 1), which we diagonalize, $C_θ = UDU^T$. U is a $d \times d$ orthonormal matrix and D is a diagonal $d \times d$ matrix consisting of the (necessarily positive) eigenvalues. The transformation matrix which whitens $θ$ is then $UDU^{-1}$, such that the transformed input covariates are $X = UDU^{-1}θ$, and the covariance of X is the identity matrix. Also, having a pre-whitened basis also justifies the use of a diagonal kernel matrix such as the ARD kernel in equation (15), for which it is often blindly assumed (without transforming the inputs) that the correlation among the input parameters vanishes.

Next, we consider the transformation of the band powers. The distribution of the original band powers in our training set is left-skewed. For a fixed $ℓ$ in our three-bin tomographic analysis, the resulting $3 \times 3$ matrix,

$$ \mathbf{B}_{\ell} = \begin{pmatrix} B_{\ell, 00} & B_{\ell, 01} & B_{\ell, 02} \\ B_{\ell, 10} & B_{\ell, 11} & B_{\ell, 12} \\ B_{\ell, 20} & B_{\ell, 21} & B_{\ell, 22} \end{pmatrix} \qquad (21) $$
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The data correlation matrix for the KiDS-450 analysis. We have ordered the covariance matrix in order of the tomographic labelling *ij*. Note that we have four band powers per tomographic bin, hence 6 × 4 blocks in the covariance matrix.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Five LH samples (using the maxi–min method) projected in 2D. In particular, we generate five LH samples in 8D and we scale them according to our pre-defined priors. In the figure, we show the projection in 2D for three parameters and as expected, each point occupies it corresponding row and column.</figcaption></figure>

must be positive-definite and emulating the matrix elements individually will not guarantee this.

To ensure that the 3 × 3 band power matrix remains positive-definite during the prediction phase when using the emulator, we instead build the latter on each element of the logarithm $\mathbf{B}_i$ (lower or upper triangular part, essentially all the unique elements),

$$ \mathbf{V}_{\ell} = \mathbf{R}\tilde{\mathbf{A}}\mathbf{R}^{\mathrm{T}} = \log(\mathbf{B}_{\ell}), \qquad (22) $$

where $\mathbf{B}_{\ell} = \mathbf{R}\tilde{\mathbf{A}}\mathbf{R}^{\mathrm{T}}$, $\tilde{\mathbf{A}}_{\nu\nu} = \log(\mathbf{\Lambda}_{\nu\nu})$ and $\mathbf{\Lambda}$ and $\tilde{\mathbf{\Lambda}}$ are diagonal. Moreover, since we normally assume a GP with mean zero and kernel, $\mathbf{K}$, we do an additional linear scaling such that the mean of the band powers in our training set is zero and has a standard

<table><caption>**Table 1.** Set of cosmological and systematic parameters which are used as inputs in the emulator. $\theta$ is the set of parameters (first eight parameters in the table below) used for the emulation scheme and $\beta$ is the set of the remaining four parameters which are also marginalized over.</caption><thead><tr><th>Definition</th><th>Symbol</th></tr></thead><tbody><tr><td>CDM density</td><td>Ω<sub>cdm</sub>h<sup>2</sup></td></tr><tr><td>Baryon density</td><td>Ω<sub>b</sub>h<sup>2</sup></td></tr><tr><td>Scalar spectrum amplitude</td><td>ln(10<sup>10</sup>A<sub>s</sub>)</td></tr><tr><td>Scalar spectral index</td><td>n<sub>s</sub></td></tr><tr><td>Hubble parameter</td><td>h</td></tr><tr><td>Free amplitude baryon feedback parameter</td><td>A<sub>bary</sub></td></tr><tr><td>Intrinsic alignment parameter</td><td>A<sub>IA</sub></td></tr><tr><td>Neutrino mass (eV)</td><td>Σm<sub>v</sub></td></tr><tr><td>Free amplitude (bin 1)</td><td>A<sub>1</sub></td></tr><tr><td>Free amplitude (bin 2)</td><td>A<sub>2</sub></td></tr><tr><td>Free amplitude (bin 3)</td><td>A<sub>3</sub></td></tr><tr><td>Multiplicative bias</td><td>m</td></tr></tbody></table>

deviation of one, for example, for the *i*th transformed band power,

$$ \mathbf{v}'_i \rightarrow \frac{\mathbf{v}_i - \bar{\mathbf{v}}_i}{\sigma_i} \qquad (23) $$

and the predictive mean and variance are

$$ \begin{gathered} \mathbb{E}[v_{i(*)}] = \sigma_i \mathbb{E}[v'_{i(*)}] + \bar{v}_i \\ \mathrm{var}[v_{i(*)}] = \sigma_i^2 \mathrm{var}[v'_{i(*)}]. \end{gathered} \qquad (24) $$

## 4.5 Training the emulator

We now have our training set {$\mathbf{X}, \mathbf{V}_{\ell,ij}$}. Therefore we have a set of 24 GPs due to each element of the transformed band powers. Prior to building the emulator, a crucial step is to choose a kernel function for the GP. Here we use the ARD kernel, defined in equation (15).

To ensure a good performance, we have to find the set of hyper-parameters which maximizes the marginal likelihood, as discussed in Section 3. An important ingredient is the analytical gradient of the marginal likelihood with respect to the kernel hyperparameters to guarantee convergence to the global minimum. The gradients are

$$ \begin{gathered} \frac{\partial k_{pq}}{\partial A} = \frac{2}{A} k_{pq}^{\mathrm{ARD}} \\ \frac{\partial k_{pq}}{\partial \ell_i} = k_{pq}^{\mathrm{ARD}} \frac{(\theta_{p(i)} - \theta_{q(i)})^2}{\ell_i^3}, \end{gathered} \qquad (25) $$

where *i* indicates the *i*th dimension of the problem. We use the Limited memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS-B) algorithm (Zhu et al. 1997; Press et al. 2007) along with the gradients defined above to optimize for these hyperparameters by minimizing the negative log-marginal likelihood, in equation (17), via gradient descent. However, it is a known fact that training a GP is not an easy task because the marginal likelihood has various local maxima (Rasmussen & Williams 2006). We adopt the standard approach of restarting our optimizer at different positions and we find that $N_{restart} = 5$ was sufficient in practice to ensure that we find the set of hyperparameters corresponding to the global optimum (see Fig. 5). Although this is not guaranteed, we also want to emphasize that the use of the gradients was required to find the global optimum. Once the GP is trained, the kernel parameters are fixed at the optimized values of the hyperparameters and then use equations (20) to make predictions.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** The figure shows the marginal likelihood of the GP (with 3000 training points) for the fourth band power matrix and $i = j = 2$. Note the local minimum for the third run of the optimizer. The other bars have almost the same value, hence showing that $N_{\text{restart}} = 5$ is a good choice for training the GP.</figcaption></figure>

## 4.6 The GP uncertainty

In this section, we look into propagating the GP uncertainty through the full forward model when we use the emulator. To be more specific, we seek the posterior distributions of the cosmological parameters and the two nuisance parameters ($A_{\text{IA}}$, $A_{\text{bary}}$), that is,

$$ \theta = [\Omega_{\text{cdm}}h^2, \Omega_b h^2, \ln(10^{10}A_s), n_s, h, A_{\text{bary}}, A_{\text{IA}}, \Sigma m_v] $$

and the other four nuisance parameters,

$$ \beta = [A_1, A_2, A_3, m] $$

marginalized over the probabilistic band powers. $A_1$, $A_2$, $A_3$ correspond to free parameters which determine excess noise in the autocorrelation power spectrum, while $m$ is the shear multiplicative bias parameter (Köhlinger et al. 2017). Using equation (10) and defining $v$ as the total band powers, we can write the joint posterior, $p(\theta, \beta | d)$ as

$$ \begin{aligned} p(\theta, \beta | d) &= \int p(\theta, \beta, v | d) \, dv \\ &= \int p(d | v, \beta) p(v | \theta) \, dv p(\theta) p(\beta). \end{aligned} \quad (26) $$

If $p(v|\theta)$ were a Gaussian distribution of the band power from the GP, the above integration would be a convolution of two Gaussian distributions and the likelihood part would be Gaussian.

However, in our analysis, the predictive distribution is Gaussian in each element of the logarithm of the band power matrix. For example, in Fig. 6, we show the GP mean and variance for two elements across a slice in parameter space. As previously discussed, if the GP predictions were Gaussian in the band powers, we could marginalize analytically over the theoretical uncertainty. Since they are Gaussian in each element of $V_{\ell,ij}$, we marginalize by drawing samples of the cosmological and nuisance parameters (see Fig. 7) and perform a Monte Carlo integration, which is relatively fast and approximating the joint posterior as

$$ p(\theta, \beta | d) \propto \frac{p(\theta)p(\beta)}{N_s} \sum_{i=1}^{N_s} p(d | V_{\ell,ij}). \quad (27) $$

$N_s$ is the number of random band powers drawn after computing the predictive mean and variance. We use $N_s = 20$ at every step in the MCMC to take into account the uncertainty from the GP. Recall that each band power is being modelled independently by a GP and hence the Monte Carlo integral in equation (27) requires few draws of the probabilistic band powers.

## 5 DATA COMPRESSION

The next era of weak lensing surveys such as Euclid and LSST will have ~10 tomographic bins, and with multiple band powers or correlation functions per bin, the number of summary statistics will be large, ~$10^3-10^4$. As an example, Euclid Collaboration et al. (2019) considered 100 bandpowers per bin, and 10 tomographic bins, which gives a minimum of 1000 summaries, and 5500 if cross-band powers are included. The setup, in the previous section, is not a scalable approach for these future surveys. In particular, emulating each band power is not an entirely feasible approach because one will have to train and store thousands of separate GPs and this process in itself can be quite expensive.

In this section, we show that the emulator can be used with the MOPED algorithm (Heavens et al. 2000b) which reduces the number of data points from $N$ to just $p$ numbers. $N$ is the number of data points and $p$ is the number of parameters in our model. For current weak lensing analysis, the gain is not significant (since we are working with only 24 band powers) but the method proposed in this work is expected to yield fast parameter inference in the regime of a large number of band powers, $N \sim 10^4$, with only $p \sim 10$ parameters of interest.

Here, we briefly cover the MOPED algorithm. The latter essentially finds some weighing vector, $b$ which encapsulates as much information as possible for a specific model parameter $\theta_\alpha$. This vector is then used to find linear combination of the data, $d$ such that the compressed data is

$$ y_\alpha \equiv b_\alpha^T d. \quad (28) $$

The first and subsequent MOPED vectors are given respectively by

$$ b_1 = \frac{\mathbf{C}^{-1} \mu_{,1}}{\sqrt{\mu_{,1}^T \mathbf{C}^{-1} \mu_{,1}}} \quad (29) $$

and

$$ b_\alpha = \frac{\mathbf{C}^{-1} \mu_{,\alpha} - \sum_{\beta=1}^{\alpha-1} (\mu_{,\alpha}^T \mathbf{b}_\beta) \mathbf{b}_\beta}{\sqrt{\mu_{,\alpha}^T \mathbf{C}^{-1} \mu_{,\alpha} - \sum_{\beta=1}^{\alpha-1} (\mu_{,\alpha}^T \mathbf{b}_\beta)^2}} \quad (\alpha > 1), \quad (30) $$

where $\mathbf{C}$ is the data covariance matrix and $\mu_{,\alpha}$ is the vector obtained by calculating the gradient of our theoretical model at a fiducial parameter set. In the previous applications of the MOPED algorithm, it was assumed that the covariance matrix is fixed. In our case, Köhlinger et al. (2017) constructed a covariance matrix which depends on the $m$ parameter, the multiplicative bias. In this work, we fix $\mathbf{C}$ at the average fiducial value provided<fnref n="1" /> in the data. Data compression with parameter-dependent covariance matrix has been explored by Heavens et al. (2017). If $\mathbf{B} \in \mathbb{R}^{N \times p}$ is the matrix whose columns consist of the MOPED vectors, the compressed data vector is just

$$ \mathbf{y} = \mathbf{B}^T \mathbf{x}. \quad (31) $$

By construction, the MOPED vectors $b_\alpha$ are orthogonal to each other, that is, $b_\alpha^T C b_\beta = \delta_{\alpha\beta}$. Therefore, the covariance matrix of $\mathbf{y}$,

<footnote n="1">¹Cosmological parameter inference depends mildly on the parameter $m$.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The left-hand plot shows the predicted band power across a slice in parameter space. In other words, we choose a point within the prior box and compute the GP mean, variance, and the actual band power for $\Omega_{cdm}h^2 \in [0.05, 0.40]$. The same procedure is repeated in the right-hnd plot, but we instead choose a point from a training set, to illustrate the fact that the predicted GP uncertainty tends to zero near the training point and the predictive variance increases towards the edge of the prior box.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** The full forward model can be understood as follows: at each step in the inference procedure, a random set of samples of the cosmological, $\theta$ and nuisance, $\beta$ is drawn from the prior, followed by a random realization of the probabilistic band powers, centred on its mean and variance before computing the likelihood. Note that the kernel hyperparameters are fixed to their optimized values.</figcaption></figure>

$\mathbf{B}^T\mathbf{CB} = \mathbf{0}$, the identity matrix, of size $p \times p$. As a result of this orthogonality condition, elements from the compressed data vector are uncorrelated. Hence, the log-likelihood is straightforwardly

$$ \log \mathcal{L} = -\frac{1}{2} \sum_{\alpha=1}^{p} (y_{\alpha} - \mathbf{b}_{\alpha}^{T} \mu)^2 + \text{constant}, \quad (32) $$

where $\mathbf{b}_{\alpha}^{T} \mu$ is computed using the emulator. The fact that the likelihood of the compressed data involves only $O(p)$ operation makes parameter inference very fast since the $O(N^3)$ operation in the standard likelihood is completely eliminated, provided $\mathbf{b}_{\alpha}^{T} \mu$ can be rapidly computed.

By emulating the MOPED coefficients directly with separate GPs, we have a very powerful tool. The GPs are still functions of just eight parameters (six cosmological and two systematics) and we now have only 11 separate GPs. Crucially, this setup is interesting because

increasing the number of band powers (e.g. in forthcoming lensing surveys) will not affect the MOPED timings at all.

# 6 RESULTS

Fig. 6 shows two band powers, evaluated across the $\Omega_{cdm}h^2$ slice in parameter space. In particular, the function in black corresponds to the accurate solver, CLASS while the broken red function corresponds to the GP mean, with the tan shading giving the $3\sigma$ credible interval of the GP. Note also that the right-hand panel shows the GP prediction through a given training point and as expected, the GP uncertainty tends to zero. As seen in Fig. 6, the GP is able to predict the band powers quite well.

Since the predictive function is a Gaussian distribution, we can build a simple emulator by just using the mean, or propagate the uncertainty from the GP through the model. Either method gives reasonable posterior densities as shown in Fig. 8. On a high-end desktop computer, the evaluation is quite fast. Computing one likelihood with the mean of the GP takes 0.03 s compared to 0.09 s if we include the GP uncertainty with 20 Monte Carlo samples to marginalize over the GP uncertainty. On the other hand, CLASS takes 0.65 s for a likelihood evaluation. If we use 1000 training points, this yields an overall speed-up by a factor of ~12–30 depending on whether we use the mean or the GP variance. In our case, we generate 360 000 MCMC samples using EMCEE (Goodman & Weare 2010; Foreman-Mackey et al. 2013) for which the full simulator takes about 44 h, while the GP emulator, using the mean, takes about ~1.5 h. On the other hand, when we emulate the MOPED coefficients using 1000 training points, each likelihood computation takes ~0.03 s with either the mean or the variance of the GPs. As an example, with the MOPED compression, CLASS takes ~44 h to generate 330 000 MCMC samples (note that there is no significant improvement in speedup because we have just 24 band powers and each likelihood computation with or without the MOPED compression is almost the same). However, with the emulator, we obtain the same number of MCMC samples in ~1.5 h with either the mean or variance of the GPs. All experiments with EMCEE were run on a single core. An interesting additional feature for the emulation scheme would be to exploit parallelization to speedup inference further.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** The full 1D and 2D marginalized posterior distributions obtained using three different methods – the one in tan colour corresponds to posterior distributions with the full simulator (CLASS), while the solid brown one corresponds to the GP emulator when random functions of the band powers are drawn, hence marginalizing over the GP uncertainty. The posterior in blue shows the distributions obtained when only the mean of the GP was used in the inference routine. The contours denote the 68 per cent and 95 per cent credible interval, respectively. Note that some parameters are dominated by their respective priors and are not constrained at all. A similar conclusion was drawn by (Köhlinger et al. 2017). However, the important point here is that the posterior from the GP is close to that obtained with CLASS.</figcaption></figure>

The distribution of the log-posterior (up to a normalization constant) of the MCMC samples obtained by using CLASS (in pale blue) is shown in Fig. 9. In the same plot, the red and green histograms show the distribution of the log-posterior when using the mean and error from the GP respectively. In the same figure in panel (b), we show the log-posterior of the samples obtained after compressing the data using the MOPED formalism. Note that, the distribution of the log-posterior of the different MCMC samples gives an indication of how faithful the function reconstruction with the GP is. With a small number of training points, there is a small shift of the log-posterior distribution of the GP emulator (either with the mean or the uncertainty) relative to the CLASS distribution.

To compare the two distributions, we compute the Kullback–Leibler (KL) divergence between the CLASS distribution and the GP distribution, that is,

$$D_{\text{KL}}(p \| q) = \sum p(\theta, \beta | d) \log \left[ \frac{p(\theta, \beta | d)}{q(\theta, \beta | d)} \right] \quad (33)$$

where $p(\theta, \beta | d)$ and $q(\theta, \beta | d)$ are the posterior probabilities computed using CLASS and the GP at the same points in parameter space. Since the posterior probability is cheap to compute with the GP, we use all the MCMC samples obtained using CLASS to compute $q(\theta, \beta)$. The KL-divergence in nats, as a function of the number of training points, is shown in Table 2. In general, the reconstruction of the band powers is almost perfect as the number of training points
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Samples of the log-posterior obtained with the three methods investigated – in panel (a), the pale blue histogram refers to the log-posterior samples from CLASS, while the red and green step histogram correspond to the mean and error of the GP respectively. A similar plot is shown in panel (b) but after applying the MOPED compression step.</figcaption></figure>

<table><caption>**Table 2.** Computational cost comparison between CLASS and the GP emulator.</caption><thead><tr><th>N<sub>train</sub></th><th>Training</th><th>MCMC (mean)</th><th>MCMC (uncertainty)</th><th>D<sub>KL</sub> (mean)</th><th>D<sub>KL</sub> (uncertainty)</th></tr></thead><tbody><tr><td>1000</td><td>20</td><td>84</td><td>216</td><td>0.84</td><td>1.00</td></tr><tr><td>1500</td><td>48</td><td>85</td><td>290</td><td>0.63</td><td>0.89</td></tr><tr><td>2000</td><td>92</td><td>86</td><td>396</td><td>0.60</td><td>0.81</td></tr><tr><td>2500</td><td>139</td><td>88</td><td>524</td><td>0.47</td><td>0.68</td></tr><tr><td>3000</td><td>209</td><td>90</td><td>692</td><td>0.09</td><td>0.65</td></tr></tbody></table>

*Notes*: the training and sampling time (columns 2–4) are given in minutes and the KL divergence is computed in units of nats (scaled by a constant); the largest $D_{KL}$ is with the 1000 training points when we include the GP uncertainty, in which case, $D_{KL} = 2.03 \times 10^{-12}$.

increases. This can also be deduced from the fifth column in Table 2 where the KL divergence decreases with increasing training points. If one could afford additional simulations, one option would be to just use the mean of the GP to sample the posterior distribution since it is not only faster compared to the case where the GP uncertainty is included, but is also closer to the actual true posterior distribution.

To assess the convergence of our MCMC chains, we also compute the Gelman–Rubin statistics (Gelman & Rubin 1992) for different scenarios. The latter is simply defined as $\hat{R} = V/W$, where V is the between-chain variance and W is the within chain variance. $\hat{R}$ is calculated for different cases, for example, for a fixed number of training points, we use the MCMC samples using the GP (mean) and the MCMC samples obtained using CLASS. This is repeated with the MCMC samples where the GP uncertainty is included. In all cases, we apply a threshold of 1.05 to ensure that the chains satisfy the ergodicity condition.

We are also interested in the $S_8 = \sigma_8 \sqrt{\Omega_m/0.3}$ cosmological parameter constraint. Recall that the GPs for sampling the posterior are built using the eight parameters (six cosmological and two systematics) and they do not allow us to predict $\sigma_8$ directly. However, the latter is a function of just the six cosmological parameters, since it involves an integration over the power spectrum. Therefore, as we compute the band powers with the 1000 training points, we also record the $\sigma_8$ values, as generated by CLASS. We then construct a training set with inputs

$$[\Omega_{\text{cdm}}h^2, \Omega_{\text{b}}h^2, \ln(10^{10} A_s), n_s, h, \Sigma m_v]$$

which is then used to build an additional GP for $\sigma_8$. This then allows us to predict $\sigma_8$ at any point in the parameter space within the prior box. We find that it takes only 1 min to predict $\sigma_8$ for 360 000 MCMC samples.

In Fig. 10, we show the 2D marginalized posterior distribution of the derived parameters $S_8$ and $\Omega_m$ using three different methods is shown. In particular, we compare the distribution obtained from CLASS with the mean and uncertainty of the GP and we conclude that we are able to recover comparable posterior densities for these two quantities, $S_8$ and $\Omega_m$.

In high dimensions, the GP uncertainty inflates between any two points. It is expected that adding more training points will improve the performance of the emulator (either with the mean or GP uncertainty) since the reconstruction of the emulated function will converge to the original function. In general, with increasing number of training points, the GP uncertainty will also decrease. The effect of the number of training points is indicated by the values of the KL-divergence in the last two columns of Table 2. However, we empirically found that the KL-divergence when we use the mean of the emulator, is always better compared to the GP uncertainty.

One might expect the inclusion of the GP uncertainty to broaden the likelihoods, so the KL divergence would not be an appropriate measure of success. However, this does not appear to be the case: marginal errors are not noticeably increased. Our conclusion is that inclusion of the GP uncertainty does not improve results, but this might vary with application. The reason is probably that we are emulating a precise function, where the training points have
<figure><img src="image_12.png" /><figcaption>**Figure 10.** $S_8$ versus $\Omega_m$ plane for our analysis. The left-hand panel shows that the GP emulator, which is a function of our cosmological parameters, for computing $\sigma_8$ is accurate and precise enough compared to CLASS. The middle panel shows the constraints without MOPED compression, while the right-hand panel includes MOPED compression. The inner and outer contours correspond the 68 per cent and 95 per cent credible interval respectively.</figcaption></figure>

zero error, and in this circumstance, the GP (which makes some assumptions that do not hold in detail) provides an error that is only approximately correct (Karvonen et al. 2020; Wang 2020).

## 7 CONCLUSIONS

We have designed a principled and detailed GP emulator for constraining not only weak lensing cosmological parameters but also the nuisance parameters. In summary, for this problem, (1) the (expensive) solver is queried a few thousand times only, to generate a training set (compared to a conventional MCMC routine where the solver is queried at every likelihood computation, (2) the emulator is ~20 times faster compared to the full solver and this makes inference very quick and (3) by emulating the MOPED coefficients, the number of separate GPs is equal to the number of parameters in the model and inference, irrespective of the number of data points, and is a very powerful technique for analysing large data sets. Moreover, the posterior distributions obtained from the emulator are quite robust compared to the full run of the simulator, with and without MOPED.

We have also demonstrated that the emulator can be used to emulate the MOPED coefficients directly. Both combined are expected to accelerate cosmological parameter inference. Emulating the MOPED compressed data has two major advantages. The first is a feature of MOPED itself, that the compressed data set does not grow at all as the original data set increases in size, so scales exceptionally well to Euclid and LSST. The second is that MOPED is only fast if the theoretical values of the MOPED coefficients can be computed very quickly. The GP provides this functionality. This is the most important conclusion of this paper.

In addition, we have used the KL-divergence as a metric to assess the performance of the emulator in obtaining reliable high-dimensional posterior distributions. As evident from Table 2, the larger the number of training points, the better the reconstruction of the emulated function and hence the lower the KL-divergence between the accurate CLASS posterior distribution and the emulator posterior distribution.

We also recommend using the mean of the emulator for this application. In Table 2, the KL-divergence between the emulator posterior and the CLASS posterior shows that the mean is always better than the emulator with the GP uncertainty. From a computational perspective, this has various other advantages, for example, inference is very fast since the GP mean prediction requires $O(N)$ operations (recall that the GP mean is a linear predictor) and storage.

An exciting application of this emulator can be in the case where one requires non-Limber computation of the power spectra. This certainly applies to galaxy clustering statistics (Fang et al. 2019) and the weak lensing bispectrum (Deshpande & Kitching

2020), even if for the weak lensing power spectrum it is a good approximation (Kitching & Heavens 2017; Kilbinger et al. 2017). In general, the latter is computationally expensive to be calculated accurately, especially at large $\ell$ because of the rapid oscillations of the spherical Bessel functions (Lemos, Challinor & Efstathiou 2017). For example, if the CLASS run were to be repeated without the Limber approximation, the emulator would have been ~$10^3$ times even faster. In future surveys, because the number of tomographics bins will be large, one would require more power spectra computations. For example, 10 tomographic bins lead to 55 auto- and cross-power spectra and the emulator would be ~$10^3$ and ~$10^5$ faster with and without the Limber approximation respectively.

Emulation has various other key advantages, apart from speeding up inference. As an example, one has to choose a good proposal distribution, which often requires tuning, to run an MCMC chain with the full simulator. The emulator can be used to explore the parameter space quickly and learn a suboptimal proposal distribution which can then be used with the full simulator.

The accuracy of the reconstructed function can be improved by adding more training points as we have demonstrated. However, scaling GPs to large number of training points results in a major computational bottleneck, mainly due to $O(N^3)$ operations in training and $O(N^2)$ in predicting the uncertainty (see Fig. 11). Fortunately, here a few hundred training points suffices to give cosmological results with only a few per cent degeneracy in error bars. Moreover, in this work, the training points have been placed according to the prior range itself. However, the interpolation scheme can be improved if we have more constrained parameters where we can use better prior information such as a Fisher matrix to intelligently place the training points. Alternatively, one can also do a quick optimization to find the maximum-likelihood estimator and the Hessian matrix, both of which can be used to construct an optimal design for the training points.

An alternative option to accelerate the computation of GP uncertainty is to intelligently partition the training set by using a clustering algorithm, for example, $k$-means clustering (Hastie, Tibshirani & Friedman 2001). During the prediction step, one can then use a local expert, which has a smaller kernel size, to compute the GP uncertainty swiftly.

A quantity which is often required in optimization or Monte Carlo methods such as Hamiltonian Monte Carlo is the gradient with respect to the negative log-likelihood (cost function). Conveniently, the gradient with respect to the mean of the GP surrogate model is analytic and this opens a new avenue towards accelerating gradient computation as well.

GPs should not only be interpreted as a method for just accelerating computations. Instead, it effectively allows us to compute the
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Illustrating the performance of the emulator as a function of the number of training points. The expensive computations reside in training and predicting the GP uncertainty. Sampling the posterior with the GP mean is quick, even with the 3000 training points. The graphs do not perfectly follow the expected scaling with *N* because of various overheads.</figcaption></figure>

posterior distribution of a function by placing a prior over it. In this work, the EE band powers and MOPED coefficients are modelled independently as GPs and we have shown that we can recover robust cosmological parameters, whilst still marginalizing over the nuisance parameters.

## ACKNOWLEDGEMENTS

We thank the referee for helpful and meaningful comments. AM is supported by the Imperial College President's Scholarship. We thank Jonathan Pritchard for suggesting the use of LH samples for building the emulator and Pat Scott, Daniel Jones for useful discussion. We also thank Zafiirah Hosenie for providing useful suggestions to improve this manuscript. We also thank Prof. Marc Deisonruth and Prof. David van Dyk for insightful discussion at the beginning of this project.

## DATA AVAILABILITY

The code and data products underlying this article will be made available at https://github.com/Harry45/gp_emulator.

## REFERENCES

* Agarwal S., Abdalla F. B., Feldman H. A., Lahav O., Thomas S. A., 2012, MNRAS, 424, 1409

* Agarwal S., Abdalla F. B., Feldman H. A., Lahav O., Thomas S. A., 2014, MNRAS, 439, 2102

* Alsing J., Heavens A., Jaffe A. H., 2017, MNRAS, 466, 3272

* Alsing J., Heavens A., Jaffe A. H., Kiessling A., Wandelt B., Hoffmann T., 2016, MNRAS, 455, 4452

* Alsing J., Wandelt B., Feeney S., 2018, MNRAS, 477, 2874

* Anderson L. et al., 2014, MNRAS, 441, 24

* Auld T., Bridges M., Hobson M. P., Gull S. F., 2007, MNRAS, 376, L11

* Bartelmann M., Schneider P., 2001, Phys. Rep., 340, 291

* Barz B., Denzler J., 2019, preprint (arXiv:1901.09054)

* Betoule M. et al., 2014, A&A, 568, A22

* Bird S., Rogers K. K., Peiris H. V., Verde L., Font-Ribera A., Pontzen A., 2019, J. Cosmol. Astropart. Phys., 2019, 050

* Bond J. R., Jaffe A. H., Knox L., 1998, Phys. Rev. D, 57, 2117

* Carnell R., 2012, R package version 0.10

* Castro P. G., Heavens A. F., Kitching T. D., 2005, Phys. Rev. D, 72, 023516

* Charnock T., Lavaux G., Wandelt B. D., 2018, Phys. Rev. D, 97, 083004

* Deshpande A. C., Kitching T. D., 2020, Phys. Rev. D, 101, 103531

* Euclid Collaboration et al., 2019, preprint (arXiv:1910.09273)

* Fang X., Krause E., Eifler T., MacCrann N., 2019, J. Cosmol. Astropart. Phys., 5, 1

* Fendt W. A., Wandelt B. D., 2007a, preprint (arXiv:1910.09273)

* Fendt W. A., Wandelt B. D., 2007b, ApJ, 654, 2

* Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

* Geenens G., et al., 2011, Stat. Surv., 5, 30

* Gelman A., Rubin D. B., 1992, Stat. Sci., 7, 457

* Goodman J., Weare J., 2010, Commun. Appl. Math. Comput. Sci., 5, 65

* Gutmann M. U., Corander J., 2015, preprint (arXiv:1501.03291)

* Habib S., Heitmann K., Higdon D., Nakhleh C., Williams B., 2007, Phys. Rev. D, 76, 083503

* Hastie T., Tibshirani R., Friedman J., 2001, The Elements of Statistical Learning (Springer series in Statistics), Springer, New York Inc

* Heavens A., Refregier A., Heymans C., 2000a, MNRAS, 319, 649

* Heavens A. F., Jimenez R., Lahav O., 2000b, MNRAS, 317, 965

* Heavens A. F., Sellentin E., de Mijolla D., Vianello A., 2017, MNRAS, 472, 4244

* Heitmann K., Higdon D., White M., Habib S., Williams B. J., Lawrence E., Wagner C., 2009, ApJ, 705, 156

* Heitmann K., Lawrence E., Kwan J., Habib S., Higdon D., 2014, ApJ, 780, 111

* Heitmann K., White M., Wagner C., Habib S., Higdon D., 2010, ApJ, 715, 104

* Hinshaw G. et al., 2007, ApJS, 170, 288

* Hirata C. M., Seljak U., 2004, Phys. Rev. D, 70, 063526

* Hu W., 2000, Phys. Rev. D, 62, 043007

* Jaffe A. H. et al., 2001, Phys. Rev. Lett., 86, 3475

* Johnson M., Moore L., Ylvisaker D., 1990, J. Stat. Plan. Infer., 26, 131

* Karvonen T., Wynne G., Tronarp F., Oates C. J., Särkkä S., 2020, preprint (arXiv:2001.10965)

* Kendall A., Gal Y., 2017, preprint (arXiv:1703.04977)

* Kern N. S., Liu A., Parsons A. R., Mesinger A., Greig B., 2017, ApJ, 848, 23

* Kilbinger M., 2015, Rep. Prog. Phys., 78, 086901

* Kilbinger M. et al., 2017, MNRAS, 472, 2126

* Kitching T. D., Heavens A. F., 2017, Phys. Rev. D, 95, 063522

* Köhlinger F. et al., 2017, MNRAS, 471, 4412

* Lawrence E., Heitmann K., White M., Higdon D., Wagner C., Habib S., Williams B., 2010, ApJ, 713, 1322

* Leclercq F., 2018, Phys. Rev. D, 98, 063511

* Leclercq F., Enzi W., Jasche J., Heavens A., 2019, MNRAS, 490, 4237

* Lemos P., Challinor A., Efstathiou G., 2017, J. Cosmol. Astropart. Phys., 2017, 014

* Lesgourgues J., 2011, preprint (arXiv:1104.2932)

* Limber D. N., 1953, ApJ, 117, 134

* Loverde M., Afshordi N., 2008, Phys. Rev. D, 78, 123506

* LSST Science Collaboration et al., 2017, preprint (arXiv:1708.04058)

* Manrique-Yus A., Sellentin E., 2020, MNRAS, 491, 2655

* McKay M. D., Beckman R. J., Conover W. J., 1979, Technometrics, 21, 239

* Planck Collaboration et al., 2014, A&A, 571, A16

* Planck Collaboration et al., 2016, A&A, 594, A13

* Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 2007, Numerical Recipes 3rd Edition: The Art of Scientific Computing. Cambridge Univ. Press, Cambridge

* Rasmussen C. E., Williams C. K. I., 2006, Gaussian Processes for Machine Learning, MIT Press, Cambridge, MA, USA

* Rogers K. K., Peiris H. V., Pontzen A., Bird S., Verde L., Font-Ribera A., 2019, J. Cosmol. Astropart. Phys., 2019, 031

* Schmit C. J., Pritchard J. R., 2018, MNRAS, 475, 1213
Schneider M. D., Holm Ó., Knox L., 2011, ApJ, 728, 137

Seppala L. G., 2002, in Tyson J. A., Wolff S., eds, Improved optical design for the Large Synoptic Survey Telescope (LSST), SPIE, p. 111

Smoot G. F. et al., 1992, ApJ, 396, L1

Spergel D. N. et al., 2003, ApJS, 148, 175

van Daalen M. P., Schaye J., Booth C. M., Dalla Vecchia C., 2011, MNRAS, 415, 3649

Wang W., 2020, preprint (arXiv:2002.01381)

Weinberg D. H., Mortonson M. J., Eisenstein D. J., Hirata C., Riess A. G., Rozo E., 2013, Phys. Rep., 530, 87

Zhu C., Byrd R. H., Lu P., Nocedal J., 1997, ACM Trans. Math. Softw., 23, 550

This paper has been typeset from a TFormula file prepared by the author.