{'fraction_non_alphanumeric': 0.029549884813042708, 'fraction_numerical': 0.009170653907496013, 'mean_word_length': 4.713110604910149, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cosmological Parameter Estimation and Inference using Deep Summaries

The ability to obtain reliable point estimates of model parameters is of crucial importance in many fields of physics. This is often a difficult task given that the observed data can have a very high number of dimensions. In order to address this problem, we propose a novel approach to construct parameter estimators with a quantifiable bias using an order expansion of highly compressed deep summary statistics of the observed data. These summary statistics are learned automatically using an information maximising loss. Given an observation, we further show how one can use the constructed estimators to obtain approximate Bayes computation (ABC) posterior estimates and their corresponding uncertainties that can be used for parameter inference using Gaussian process regression even if the likelihood is not tractable. We validate our method with an application to the problem of cosmological parameter inference of weak lensing mass maps. We show in that case that the constructed estimators are unbiased and have an almost optimal variance, while the posterior distribution obtained with the Gaussian process regression is close to the true posterior and performs better or equally well than comparable methods.


I. INTRODUCTION
Parameter estimation and inference are common problems in many fields of physics. In parameter estimation, one tries to construct unbiased point estimators for model parameters and in parameter inference one is interested in the posterior distribution of the underlying model parameters, given a set of observations. These two problems are often difficult to solve in real-world scenarios for multiple reasons. There are usually only very few relevant parameters in physical models that are used to describe highly complex systems whose interactions are governed by physical laws. Due to the nature of these complex interactions it is rarely possible to directly observe the underlying model parameters, but rather one has to extract them from high dimensional data. A prime example for this is the field of cosmology or, more specifically, weak graviational lensing (see e.g. [54] for a review). The currently most accepted ΛCDM (dark energy, cold dark matter) cosmological model contains only five relevant parameters that have to be extracted from the shape measurements of millions of galaxies. The relationship between the model parameters and the observed data is highly non-linear and not analytically tractable. A practical way to extract the desired parameters is via complex models and large cosmological simulations that can require millions of CPU and GPU hours to run [51]. Given the large amount of data produced by these simulations, a standing problem is the compression of the observed data, such that one can optimally constrain the cosmological parameters. Frequently used methods for the data compression include MOPED and (nuisance hardened) score compression [2,3,21]. However, the recent advances in machine learning offer new possibilities in the field of cosmology. Especially deep neural networks, have gained a lot of attention [10,52] and have been used in particular for parameter estimation [22] and inference [17,26].
An important problem in the construction of parameter estimators is the quantification of their bias. For deep neural networks there has been extensive work concerning uncertainty estimation [18]. Common approaches include drop out [19] and Bayes by Backprop [9]. These methods help to quantify the intrinsic model uncertainties, however, they do not quantify potential biases of the models. In [12] they are able to construct unbiased estimators, but the presented method is not applicable for general non-linear relationships between model parameters and observations that are considered in this work.
In parameter inference one tries to evaluate the posterior distribution and consequently confidence intervals of the underlying model parameters given a set of observations. If the posterior distribution can be evaluated analytically, one can resort to classic methods such as Markov chain Monte Carlo (MCMC) to obtain the desired confidence intervals. However, even if the posterior distribution is known, it might be unfeasible to perform a MCMC run if the likelihood requires the evaluation of expensive simulators. Approaches using Gaussian process regression have been developed [48,58] in order to reduce the number of required likelihood evaluations and will be discussed in more detail below. Such schemes can potentially reduce the number of necessary likelihood evaluations by orders of magnitude. However, they are not applicable if the likelihood is either unknown or can't be evaluated. In these cases one can resort to likelihood free inference (LFI) methods such as approximate Bayesian computation (ABC) [15]. ABC methods make use of a discrepancy measure between observed and simulated data to estimate the posterior distribution. Given a sufficient amount simulations it is theoretically possible to obtain the true posterior via this measure and a kernel function. However, ABC methods, while accurate and stable, again require a large number of simulations that might be unfeasible to run in a reasonable time.
Combinations of Gaussian process regression and ABC have also been proposed, e.g. to emulate the discrepancy measure between simulations [25]. This emulation can greatly reduce the required number of simulations, as the posterior can be estimated directly through the emulated discrepancy measure. Further, [24] analysed Gaussian process regression of noisy log-likelihood estimates. This was primarily done in a synthetic likelihood setting, where the unknown likelihood is approximated by a Gaussian distribution. Using simulations, it is possible to estimate the Gaussian approximation of the likelihood, as well as, its uncertainty. This approximation can then be emulated with Gaussian processes to obtain the posterior. In [24] they then focus on the uncertainty estimation of the emulated posterior and acquisition strategies for new simulations that optimally improve improve the posterior estimate.
There are also combinations of ABC and deep neural networks, for example learning summary statistics for ABC inference has been proposed in [16,27,61], e.g. by training the networks to predict the model parameters directly from simulated data. Further, a method to generate globally sufficient summaries of arbitrary data using mutual information has been developed in [11]. However, these approaches usually require on-the-fly simulations and do not use the generated summaries for parameter estimation.
Likelihood-free inference using neural networks has been used in [1,32,46,47]. These approaches mostly aim to learn the distribution of the (compressed) data that can be used to obtain the posterior distribution of the underlying model parameters and are not considering a fixed grid of pre-computed simulations. An example of such an approach using neural density estimation to cosmological data can be found in [4]. Lastly, in [33] they examine the compression of weak lensing data similiar to the approach in this work, the key difference being that they have access to on-the-fly simulations and use neural density estimation for the inference instead of our approach that is explained below.
In this work, we consider a scenario that is especially common in weak graviational lensing, where we only have access to a fixed grid of simulations generated prior to the analysis, additional simulations are therefore not obtainable, making it impossible to use acquisition strategies like the ones discussed in [24]. The problem we address is to construct estimators of the underlying model parameters and to calculate confidence intervals of these parameters given a set of observations without having access to the likelihood function. An overview of the general approach is shown in Figure 1 where we distinguish between two parts of the pipeline: i) parameter estimation, where we construct an estimator of the model parameters and ii) parameter inference where we find confidence regions for the model parameters using the estimator constructed Overview of the presented approach that consists of two parts. In the first part we construct parameter estimatorŝ θ using the highly compressed summariesŝ(x) of the simulated data x. And in the second part we use these estimators to calculate the ABC posterior pABC which is then emulated with a Gaussian process regression to find the confidence regions of the model parameters.
in i). To achieve this, we start by relaxing the prior assumptions about the distribution of the compressed data made in [10] on their loss function. Next, we develop a method to construct estimators of the underlying model parameters with quantifiable biases using the compressed simulations. Finally, we demonstrate how one can utilize the constructed estimators to calculate ABC log-posterior estimates and their corresponding uncertainties, which can in turn be used for parameter inference using Gaussian process regression. Each of these aspects is important to run the pipeline as a whole in a realistic setting and we validate our findings using a cosmological model that is comparable to modern weak lensing surveys (e.g. [62]). Finally, the code used in this work is available on github [Note1].

II. PROBLEM DESCRIPTION
We start with a detailed description of the two problems of parameter estimation and inference. In parameter estimation, the goal is to estimate the unobserved, usually low dimensional, model parameter θ ∈ Θ ⊆ R d θ from observed, high dimensional data x ∈ X ⊆ R dx . The constructed estimatorθ(x) should be unbiased, i.e.
and minimize the variance around the true value, i.e.
Such an estimator does not always exist. For example, if we have a family of distributions which fulfills p(x|θ) = p(x|θ ) for θ = θ , it is impossible to construct an unbiased estimatorθ(x) for all θ ∈ Θ. We therefore only consider distributions p(x|θ) that allow the construction of such an estimator. In parameter inference the goal is to calculate the posterior distribution p(θ|x) of the model parameter given an observation x and subsequently to find confidence intervals for the model parameter. In general this is difficult, because the likelihood p(x|θ) is often unknown. However, we assume that there exists a computationally expensive simulator, which can produce samples of the observed space x ∼ p(x|θ) for any possible model parameter θ. Throughout this work, we will consider an experiment design where the computationally expensive simulator was run on a set of n θ different model parameters and n x samples were drawn from the distribution p(x|θ α ) for each parameter θ α in the set. We will use Greek indices to indicate parameters contained in a set (e.g. θ α ∈ R d θ ) and Latin indices (e.g. θ i ∈ R) to refer to individual elements of vectors and matrices.

III. PARAMETER ESTIMATION
We start with a description of the first part of the pipeline shown in Figure 1, where our goal is to construct an estimator of the underlying model parameters. A simple way to construct such an estimatorθ(x) using the generated simulations is by doing regression. One can train a neural network h(x, w), parameterized by weights w ∈ R dw to predict the underlying model parameter θ in a supervised fashion using, for example, the mean squared error (MSE) over mini-batches of size N b , where θ True α is the true parameter of the randomly sampled simulation x α . However, using this approach, there is no guarantee that the trained network h(x, w) is an unbiased estimator of the model parameter. In fact, the resulting estimatorθ(x) = h(x, w * ), where w * are the trained weights obtained by minimizing the loss function of equation (3), will approximate the expected value E p(θ|x) [θ] [27] and can be arbitrarily biased. In our setting with a fixed grid of simulations an additional source of bias can be introduced by the choice of simulated parameters θ α . By choosing a fixed grid of parameters θ α one implicitly places a prior over the parameters the network tries to estimate. In [52] and [60] it was shown that the network can learn this prior, potentially leading to large biases in the predictions of the network in areas where there is a lower parameter sampling density.
To alleviate this issues, we propose to split the problem into two parts. We start by training a network to produce informative, low-dimensional summary statisticsŝ(x) = h(x, w) ∈ R ds , with d θ ≤ d s d x . In a second step, we use these summary statistics to construct an estimator θ(x) ≡θ(ŝ(x)).

A. Obtaining Informative Summaries
As mentioned earlier, we face the problem of producing a compact representation of the data, which we name summary (corresponding to step 1 ○ in Figure 1). To achieve this, [10] proposed to optimize the Cramér-Rao bound that expresses a lower bound on the variance of an estimator as a function of the Fisher information, starting with the assumption that the summary statistics generated by a neural network follow a Gaussian distribution. We will next show that this assumption is not necessary. Formally, the Fisher Information Matrix (FIM) is defined as where θ i denotes the i-th entry of the underlying model parameter θ ∈ R d θ . The Fisher information matrix is used in the Cramér-Rao bound, which bounds the covariance matrix of any summaryŝ(x) as where the derivative with respect to θ is understood as a Jacobian matrix, A B indicates that A − B is positive definite and we defined Note, that if the summaryŝ(x) is itself an unbiased estimator of θ, the Jacobian becomes the identity and the covariance matrix is bounded directly by the Fisher information matrix. A summary that achieves equality is called sufficient, if its Jacobian with respect to the underlying model parameter has full rank [53]. This type of sufficiency also implies Bayesian sufficiency, i.e.
which means that the summary retains all relevant information about the model parameters. Similarly to the work of [8] we will only consider summary statistics having the same dimensionality as the model parameter d s = d θ . However, in the appendix A we explain how the same approach can be used for summaries with d s > d θ . Assuming that the Jacobian is invertible, one can transform equation (5) into The right side of this inequality can be used as a loss function to construct informative summaries. Note that this loss is equivalent to the log-determinant of the inverse Fisher matrix used in [10], without relying on the assumption that the summary follows a Gaussian likelihood, showing the generality of the approach. Estimating the covariance matrix of the summaries in equation (8) is straightforward. The simulator, however, is generally not differentiable, preventing the exact calculation of the Jacobian. As in [10], we calculate the Jacobian via finite differences. One starts by generating a set of simulations using a fiducial parameter [Note2] θ fid and small perturbations of it θ ± fid = θ fid ± ∆θ ± , resulting in 2d + 1 parameters. Using this set of parameters one can approximate the Jacobian as where the expectation Ψ θ (ŝ) is empirically calculated over multiple samples from p(x|θ). This reduces the variance at the price of an increased computational effort spent on a single step during the training. However, training the network only around a single (or a few) fiducial points leads to the advantage that the remaining datapoints can be used for testing and validation. It also reduces the overall training time compared to simple regression networks as one does not have to cycle through the entire data set of simulations.

Regularization
The minimization of the righthand side of inequality (8) could be numerically unstable because the loss allows the entries of the Jacobian or the summary statistics itself to become arbitrary large or small. The reason for this is that there are infinitely many summary statistics that have the same Fisher information. In fact, any one-toone transformation of summary statistics preserves the Fisher information, showing the degeneracy of the possible solutions. To cope with this issue we add a regularization term to the loss where λ is a hyper-parameter, I d represents the ddimensional identity matrix and we use the Frobenius norm as matrix norm. Additionally one can add a regularization term proportional to the length of the summary vectors ŝ(x) to avoid arbitrarily large summaries. However, we never observed any sort of over-or underflow in the summary statistics during our experiments and did not use this type of regularization.

Adding Expert Knowledge
We briefly mention that the approach described above allows for an easy incorporation of expert knowledge into the summaries. Consider the scenario where a set of summaries (ŝ 1 (x), . . . ,ŝ n (x)) generated using expert knowledge is available. One can simply incorporate this additional information into the model by concatenating the summaries to the fully connected layers. Another alternative is to perform another round of compression by training a different neural network to compress the list of summaries containing the previously trained summaries and the expert summaries (ŝ(x),ŝ 1 (x), . . . ,ŝ n (x)).

B. Constructing Estimators
Next, we use the generated summaries to construct estimators of the underlying parameters (corresponding to step 2 ○ in Figure 1). To do so, we first perform the following series expansion: where we used the Einstein summation convention, meaning that we sum over all indices that appear as lower and upper index in a term. One can choose the coefficients a i , b ij , etc. such that they minimize the bias of the resulting estimator around the chosen fiducial point θ fid . In this work, we will only consider first-order expansions, however, the construction of higher-order terms follows the same principle. We start by expanding the expectation value of the summaries around the fiducial parameter where we defined ∆θ = θ − θ fid . Plugging this into the first-order expansion of equation (11), one can find the coefficients In practice one has to estimate the coefficient a i and b ij from the sampled simulations x of the fiducial parameter, leading to an additional small error from this approximation. We briefly note here, that the construction of higher-order estimators does not only require the estimation of the expected values Ψ θ fid (ŝ) and its derivative, but also from higher-order moments. We show how to construct a second order estimator, having an error of order O(∆θ 3 ), in appendix B.
Since the construction of the estimators requires only simulations of the fiducial parameter and its perturbations, the simulations from the other parameters can be used to measure the performance of the constructed estimators.
Finally, although the approach we described uses a single estimation point, it can simply be extended to combine several estimation points. If one has access to multiple fiducial points (and perturbations of them), one can construct approximately unbiased estimators for each point. The different estimators can either come from the same underlying network or one can train a different network for every fiducial point individually. The constructed estimators can then be combined such that the bias and variance of the final estimators is minimized. For instance, [31] proposed to use a linear combination.

IV. PARAMETER INFERENCE
In this section, we show how one can use the previously described estimators to constrain the underlying model parameters using parameter inference (second part of the pipeline in Figure 1). The goal of parameter inference is to find confidence regions of the underlying model parameter θ, given an observation x obs (corresponding to step 3 ○ in Figure 1). This is no trivial task in practice, even if the posterior distribution p(θ|x obs ) is known, because it requires its integration to derive marginal distributions.
A. Common Approaches

Markov Chain Monte Carlo
A popular way of performing parameter inference with a known likelihood function is to use Markov chain Monte Carlo (MCMC). However, even if the posterior distribution p(θ|x obs ) is known, running a MCMC can be unfeasible if the evaluation of the likelihood is expensive. A possible solution to this issue was proposed by [48]. They used Gaussian process regression to emulate the expensive posterior p(θ|x obs ) and performed a MCMC on the emulated posterior to obtain the parameter constraints. The obtained one and two dimensional marginal distributions from the Gaussian Process emulation were in excellent agreement with the ground truth. However, this approach is not applicable if the analytic form of the likelihood, or good approximations of it, are not available.

Neural Density Estimators
In the case of an unknown or intractable posterior distribution p(θ|x obs ) one has to use likelihood free inference (LFI). A promising LFI technique is based on neural density estimators, as for instance proposed in [1], where a combination of mixture density models and auto-regressive flows is used to learn the distribution of the summaries p(ŝ|θ). The posterior distribution p(θ|s obs ) can then be obtained by using Bayes theorem. We use the implementation of this method, called PyDelfi [Note3], as baseline for our experiments in section V. The advantage of PyDelfi is that it can be applied to our settings as it does not require any on-the-fly simulations. Further, the constructed estimators, described in section III, are usually close to the actual parameter making it much simpler to learn the distribution p(θ(x)|θ), as opposed to summaries whose distribution is completely unknown or have a higher dimensionality.

Approximate Bayesian Computation
Another common way of performing LFI is called approximate Bayesian computation (ABC). In ABC one usually replaces the high dimensional data x with a low dimensional summaryŝ(x). The ABC approximation of the posterior (e.g. see [15]) is then given by (12) where s obs =ŝ(x obs ) is the summary of the observation, K h (x) = K(x/h)/h is a kernel density function with scale parameter h, · defines a norm and p(θ) denotes the prior distribution of the model parameters. It is straightforward to see that the ABC posterior converges to p(θ|s obs ) for h → 0. If the summaryŝ(x) is sufficient, we even have lim h→0 p ABC (θ|s obs ) = p(θ|s obs ) = p(θ|x obs ).
The simplest ABC algorithm is called rejection ABC which works similarly to MCMC. It first generates a random sample θ from p(θ) and then generates a sample x from p(x|θ). The new parameter θ is accepted if ŝ(x)− s obs < , where the threshold is a hyper-parameter corresponding to the scale h. Using this technique, one can generate samples from the ABC posterior and use them to calculate the marginal distributions and eventually the confidence regions. The downside of this algorithm is that it requires a large number of simulations x, even with a carefully selected threshold. More advanced ABC algorithms such as population Monte Carlo ABC (PMC-ABC) [23] require fewer on-the-fly simulations, but can not be applied to a fixed grid of simulations.

B. Gaussian Process ABC
A combination of Gaussian process regression and ABC has been proposed by [25], using the Gaussian process to emulate the discrepancy measure of the ABC. Later, [24] used Gaussian process regression to emulate (noisy) likelihood estimates. However, they mainly consider a synthetic likelihood setting where the unknown likelihood is approximated by a Gaussian distribution, focusing on the error estimation of the posterior and acquisition strategies for new simulations. In this work we compute ABC posterior estimates and their corresponding uncertainties using the estimator presented in Section III. Afterwards, we use Gaussian process regression to emulate the posterior similarly to [24]. This does not require on-the-fly simulations as in standard ABC methods but still allows likelihood free inference. Further, it makes it possible to directly learn the posterior p(θ|θ obs ) conditioned on a set observationθ obs and does not require to learn the conditional distribution p(θ|θ) jointly in θ andθ like neural density estimators.
We use our available n x draws from the distribution p(x|θ) and utilize the constructed estimator of the previous section. We can then obtain an unbiased, consistent and non-negative Monte Carlo estimate of the ABC posterior described by equation (12) (corresponding to step 4 ○ in Figure 1) viâ where we used our estimatesθ α ∼ p(θ|θ) as summaries and the normalization is irrelevant for our purpose. Note that the scaling parameter h can be used to smooth the ABC posterior if the fixed grid is too coarse. A good choice for the norm inside the kernel density function is the Fisher distance [10] around the fidicual parameter θ α −θ obs 2 = (θ α −θ obs ) TÎ fid (θ α −θ obs ), where we defined the estimated Fisher matrix For a first order estimator this is equivalent to the Mahalanobis distance of the summary and scales the entries of the summary according to their variances. As [48] pointed out, performing Gaussian Process regression on the posterior p(θ|θ obs ) itself is difficult. One reason for this, is that the posterior is usually close to zero in a large subset of the parameter space Θ. This can be alleviated by performing the regression on the log-posterior log(p(θ|θ obs )) instead. Using equation (14) we can construct an estimator of the ABC log-posterior where we left out a constant offset coming from the normalization constant of the ABC posterior estimate. One should note that this estimator is still consistent, but not unbiased. The bias follows directly from Jensen's inequality and becomes negligible for sufficiently large sample sizes n x . Another issue with normal Gaussian process regression is that the variance of the estimators depends on the underlying parameter θ and is not constant over the parameter space Θ. Using gpflow[Note4] [14] we adopt the following model for the Gaussian process regression (corresponding to step 5 ○ in Figure 1) where k(·, ·) is a kernel function and σ 2 α is the estimated variance of the ABC log-posterior estimate. To obtain the variances σ 2 α we use a Taylor series expansion of the moments. If g : U ⊆ R → R is a smooth function and X ∈ U a continuous random variable with finite mean and variance µ X , σ 2 X < ∞, one can approximate the variance of the transformed random variable as We find this to be a good approximation [Note5] for g(x) = log(x) if X is approximately normal distributed and µ X /σ X 2.5. Therefore, we start by estimating the variance of our ABC posterior estimates by taking the variance on both side of equation (14) Var and then propagate this variance to the log-posterior via the Taylor series expansion The evaluation of the Gaussian process with optimized kernel parameters is computationally cheap, making it straight forward to obtain all relevant properties of the posterior distribution p ABC (θ|θ obs ) with standard methods like MCMC.

V. EXPERIMENTS
To test our method, we apply it to the inference of cosmological parameters from weak gravitational lensing maps. For this purpose, we consider simplified Gaussian maps for which the sufficient summary statistics is known. In weak lensing, one measures the shapes of millions of galaxies across the sky that have been distorted by the density inhomogeneities in the universe on the line-ofsight, according to general relativity. The correlations between the measured shapes can then be used to construct convergence maps, which are maps of the mass overdensity projected on the sky. Convergence maps can be generated using large, computationally expensive, cosmological simulations. These simulations are usually so expensive to run that it is impossible to run them onthe-fly during the data analysis. Comprehensive reviews about weak gravitational lensing can be found in [7,54].
Common ways to analyze weak lensing convergence maps usually rely on two-point correlation functions (e.g. [6,55]). However, it is known that the two-point correlations function is not a sufficient summary of the data and there is an on going search for so called non-Gaussian statistics (e.g. [62]). This problem is going to become more significant as the quality and quantity of the observed data is going to increase in future surveys like Eulcid [Note6] [30]. Replacing the human made summaries with automatically learned summaries that can be used for LFI can therefore be extremely beneficent to achieve tight and robust cosmological parameter constraints.

A. Model Parameter
Our experiments are based on the flat wCDM cosmological model that extends the standard ΛCDM model with the dark energy equation of state parameter w 0 . The remaining parameters of the ΛCDM are the total matter density fraction Ω M , the baryonic matter density fraction Ω b , the Hubble parameter H 0 ≡ 100h, the spectral index n s and the fluctuation amplitude σ 8 . Here, flat is referring to the total energy density parameter being equal to one Ω = Ω CDM + Ω b + Ω Λ ≡ 1, where Ω CDM = Ω M − Ω b is the dark matter density fraction and the dark energy density fraction Ω Λ is determined by the constraint Ω = 1, which corresponds to a geometrically flat universe. We consider two cases, i) a two-dimensional model (hereafter 2D) where we only vary the parameters weak lensing is most sensitive to, namely Ω M and σ 8 ; and ii) a six-dimensional model (hereafter 6D) that allows all six parameters to vary.
We chose the same fiducial parameters and perturbations as [59] for the training, which are listed in table I. Note that we place an additional prior on the degeneracy parameter S 8 = σ 8 Ω/0.3, corresponding to the width of the degenerate contours (as shown in Figure 7).

B. Mock-Survey Settings
We consider a stage-III like survey, meaning that the convergence maps span ∼ 5000 deg 2 on the sky (see Figure 2). Further, instead of using convergence maps generated from cosmological simulations, we use Gaussian random fields (GRF). The advantage of GRF is, that while being similar to realistic convergence maps, they are easy to model statistically. A spherical GRF is completely determined by its spherical harmonic power spectrum (see appendix C for more details), which means that the power spectrum is a sufficient statistic and that one can use power spectra to generate GRF. The power spectrum of convergence maps is itself a summary statistic often used for cosmological parameter inference [6,55] and advanced codes exist that allow the fast generation of spectra (∼ 1 second) for any set of cosmological parameters. In our model we use the publicly available code CCL[Note7] (The Core Cosmology Library) described in [13] to predict the spectra.

Redshift Distribution
The generation of convergence power spectra does not only require the cosmological parameters, but also observational parameters. The most important observational parameter is the galaxy distribution along the line-of-sight of the observer, i.e. how likely it is to observe a galaxy at a certain distance (redshift) from the observer, called the redshift distribution. We use the same total redshift distribution as [62], which has an analytic formula n(z) ∝ z 1.5 exp − z 0.31 1.1 , and is a realistic redshift distribution for a stage-III cosmological survey. The normalization of the distribution depends on the redshift boundaries of the survey (maximum observed distance). For our model we chose z min = 0.0 (position of the observer) and z max = 3.5. Redshift distributions are normalized like any probability distribution by the condition zmax zmin n(z)dz = 1.
Further, we consider a survey setting where the observed galaxies are assigned into four tomographic bins according to their redshift. This allows the creation of four convergence maps, all having a different redshift distribution. Note however, that their features are correlated in a certain way, since they stem from the same universe. These correlation are defined by their cross-spectra, which can also be predicted using CCL. Each of the tomographic bins should have approximately the same number of galaxies.
Optimally, one could measure the redshift of a galaxy exactly, such that redshift distributions of the individual bins follow the same redshift distribution as the survey truncated with different z min and z max . However, in reality redshift measurements are prone to errors such that the individual redshift distributions broaden. We model this measurement error following [5] assuming that the standard error of a redshift measurement is given by σ z = δ(1 + z), with δ = 0.01. The resulting redshift distributions of the individual bins are shown in Figure 3 and do not have a closed form expression.

Mock-Data Generation
Generating a GRF with a single power spectrum is fairly straight forward. However, generating different GRF with given cross spectra requires a more advanced generation procedure. The necessary algorithm is described in appendix C of [57]. To actually generate the spherical GRF we use the Hierarchical Equal Area iso-Latitude Pixelization tool (HEALPix [Note8]) described in [20], using CCL to calculate the input power-and crossspectra of convergence maps given a set of cosmological parameters and the necessary redshift distributions. Our cosmological model used CCL to predict spectra up to max = 1000. In Healpix the resolution of a map is given by the nside parameter, which we fixed for all experiments to nside = 512, leading to ∼ 400, 000 relevant pixel (see Figure 2) for each of the four redshift bins.
We generated 10, 000 simulations with the fiducial parameter and each of 12 individual perturbations. Additionally, we generated n x = 50 simulations for n θ = 1024 different parameters sampled as a Sobol sequence inside the reported priors for the 2D model and n θ = 2048 for the 6D model.

Observational Noise
The measured data in realistic cosmological surveys contains a substantial amount of noise. The observational noise can be modeled as Gaussian random noise [50] that mainly depends on the number of observed galaxies. The distribution is centered at zero and its standard deviation is given by where σ e = 0.3 is a typically measured ellipticity dispersion, A is the area of a given pixel and n is the galaxy number density. We add observational noise to all our simulations with a chosen galaxy density of n = 12 galaxies/arcmin 2 which translates to 3 galaxies/arcmin 2 per redshift bin. Further, it can be shown that the power spectrum of a pure (pixelized) noise map is constant where f sky corresponds to the fraction of the sky that is covered by the survey and n pix to the number of pixels required to cover the whole sphere.

C. GCNN Network
The chosen survey configuration of the model makes it impossible to treat the data as flat images. There exist multiple approaches to deal with spherical data. We chose to use DeepSphere [Note9], described in [49], that treats the pixels of the the generated GRF as nodes in a graph. The distances between the nodes are given by the spherical geometry of the maps. The graph convolutional filters in DeepSphere are expressed in terms of Chebyshev polynomials and applied to the data using the graph Laplacian. Additionally it uses the hierarchical pixelization approach of HealPix to downsample the data while maintaing is spherical geometry. This downsampling procedure is refered to as pseudo convolutions and displayed in Figure 4. However, to use this feature one has to add a padding to the input of the GCNN such that  II. Architecture of the used GCNN. We report layer type, output shape (N b being the batch size) and number of trainable parameters. The residual layer (see Figure 5) is repeated ten times. The output shape depends on the model. For the 2D model we have Nout = 2 and Nout = 6 for the 6D model.
The architecture of the used network is shown in table II. All presented layer use the ReLU activation function. The architecture of the residual layers is shown as a flowchart in Figure 5. All graph convolutions use Chebyshev polynomials with degree K = 5. Further details and a concrete implementation can be found on our github [Note10] page.

Training
All graph convolutional neural networks presented in this work were trained for 100,000 steps with batch size of 128 distributed among 32 GPUs on the cluster computer Piz Daint [Note11]. This batch size has to be multiplied with 2d θ + 1, because of the necessary perturbations to estimate the derivatives with respect to the underlying parameters. Beside the loss defined in equation (8) we added the regularization term of equation (10) with parameter λ = 10. The observational noise (see section V B 3) modeled as Gaussian random noise was sampled anew for each training step, but the same seed was used for each perturbation. The weights were optimized using the Adam optimizer [28] with initial learning rate of 0.0001 and moments β 1 = 0.9 and β 2 = 0.999. Additionally, to avoid large weight updates, we applied global gradient clipping, reducing the global norm of the gradients in a training step to 5.0 if necessary.

D. Ground Truth
The advantage of this model is that a sufficient summary statistics is known from the beginning. All maps are generated from the power spectra calculated with CCL. Let C(θ) be such a spectrum, already corrected for the survey geometry (see appendix C 1). Given a spectrum C obs measured from an observation x obs , the posterior distribution is given by [56] p(θ|C(θ), C obs , S) ∝ det(S) − 1 where Q is defined as the covariance matrix S is estimated from the N = 10, 000 simulations of the fiducial parameter and C noise is the noise contribution, which does not depend on the cosmological parameters. We make the common assumption that the covariance matrix S does not depend on the cosmological parameters θ. For parameter inference, one can directly integrate the posterior distribution with a MCMC. For parameter estimation one can use the likelihood to calculate the Fisher information matrix, which in turn gives the smallest possible variance of any estimator. It can be shown [29], that the Fisher information matrix of the likelihood defined in equation (19) is given by Leading to a lower bound on the variance of any unbiased estimator Var θ i ≥ (I −1 ) ii . We calculate the lower bound for the fiducial parameters using the perturbed spectra to numerically approximate the derivatives.

E. Parameter Estimation
We trained a GCNN (see section V C) for 100, 000 training steps using a batch size of 128. Afterwards, we used the 10, 000 simulations of the fiducial parameter and its perturbations to construct a first order estimator.
Additionally, inspired by [16,27], we trained the same network with the mean squared error loss (MSE) defined in equation (3)  size of 128. The predictions of this model where used as a benchmark for the parameter estimation task. We compare the true parameters and the mean predictions of the constructed estimator of the 2D model and the predictions of the MSE model for the n x = 50 maps generated for the n θ = 1024 test parameters in Figure 6. Even though the 2D model network was only trained on the fiducial parameter and its perturbations, the mean predictions are not significantly biased throughout the considered prior range. Further, the variance of the constructed estimator is close to the calculated lower bound of the known sufficent summary. A plot of the residues, i.e. the biases, can be found in appendix D. The model trained with the MSE loss exhibits a larger bias despite being trained on the whole parameter space. This bias also increased towards the edges of our chosen priors as expected and shown in [52] and [60].

F. Parameter Inference
Using the trained network we perform a GP ABC. We evaluated the ABC posterior estimates as described in section IV using a sigmoid kernel with scale parameter h = 0.35 for the 2D model and h = 0.45 for the 6D model. We investigate the impact of the scale parameter on the resulting posterior in appendix E. For the 2D model we use the first 128 parameters of the generated Sobol sequence and fit a Gaussian process using a standard Matèrn 5/2 covariance function with automatic relevance determination (ARD). This was done with the gpflow package. We also used the parameter normalization suggested in [48], rotating and scaling the input parameters θ α such that they are uncorrelated with zero mean and unit variance. The posterior distribution obtained by running a MCMC with the fitted GP is shown in Figure 7. For the 2D model we performed an additional GP ABC using the predictions of the model trained with the MSE loss using the same parameters, the only difference being that we used the Mahalanobis distance, using the fiducial covariance matrix, as norm. The results are also shown in Figure 7.
For the 6D model, we use a Sobol sequence of n θ = 2048 parameters. In this model, most of the variance is caused by the observational noise. This can be used to effectively increase the sample size n x = 50. For each simulation we create 50 version of observational noise and use the 50×50 mock surveys to estimate the ABC posterior. Afterwards, the same Gaussian process regression as for the 2D model is performed. We provide a plot of the two-dimensional marginal posterior distributions in appendix F.
We use exactly the same data to train a neural density estimator. We use the same model as presented in section 6 of [1] only adapting the dimensions. The results can be seen in Figure 7. We also trained a neural density estimator of the summary statisticsŝ as opposed to the estimatesθ for the 2D model. However, the network did not converge to an acceptable result, showing that the distribution of the constructed estimators is easier to learn.
We assess the performance of our method in two ways: by comparing the mean of the posterior and by calculating the Jensen-Shannon divergence (JSD) with respect to the MCMC run on the true posterior. A high precision is extremely important for comsological parameter inference. All results are listed in table III. Both methods, GP ABC and PyDelfi perform very well on both models. The GP ABC, however, yields mean parameter values that are on average closer to the ground truth, as well as a lower JSD. The difference of the JSD is larger for the 2D model because PyDelfi produces a posterior distribution that is longer in the S 8 = σ 8 Ω M /0.3 direction (i.e. thickness of the degenerate contours).

VI. CONCLUSION
In this work, we presented a novel approach to construct parameter estimators from highly informative summaries. The constructed estimators have a quantifiable bias and are suitable for parameter inference using the presented GP ABC approach or the neural density estimators presented in [1]. We validated our approach with a realistic TABLE III. Overview of the parameter inference results. The number of evaluations corresponds to the number of different parameter θ used to obtain the posterior distribution. The constraints report the mean values and the 68% confidence regions. The JSD was calculated with respect to the MCMC run, which we consider as our ground-truth (lower values are better).  It is also useful to define the cross-spectra between two GRF T 1 and T 2 that might stem from different power spectra, which is usually done in terms if the complex spherical harmonic coefficients C 12 = E R(a 1 m )R(a 2 m ) + I(a 1 m )I(a 2 m ) . The cross spectrum can be estimated in the same way as the normal power spectrum and it easy to see that if C 1 = C 2 , then, the cross spectrum reduces to the normal power spectrum C 12 = C 1 = C 2 .

Survey Geometry
The C estimation described by equation (C1) is only unbiased if the spherical harmonic coefficients were calculated from a full sky map. Realistic surveys, however, do not observe the entire night sky and often have complicated survey masks (see Figure 2). Such a survey mask leads to a bias in the estimation of the power spectrum called mode-coupling. It is possible to remove this bias with the Wigner matrix formalism. It turns out that the expected value of the estimated masked sky spectrum C ,masked is connected to the underlying power spectrum via a linear transformation where W is the Wigner matrix. The Wigner matrix is completely defined by the survey mask and its calculation is described in appendix D of [57]. For continuous GRF it is sufficient to have a non-vanishing area to guarantee that the Wigner matrix is invertible. This holds true most of the time even for pixelized GRF. However, due to numerical stability it is common to apply the Wigner matrix to the underlying power spectrum (theoretical prediction) as oposed to reconstruct an unbiased estimator from the masked power spectrum (observation).