Primordial power spectrum and cosmology from black-box galaxy surveys

TL;DR

This paper introduces SELFI, a likelihood-free inference method using simulations to accurately recover the primordial power spectrum and cosmological parameters from complex galaxy survey models, outperforming existing techniques.

Contribution

The paper presents a novel likelihood-free inference framework, SELFI, that efficiently utilizes numerical simulations for unbiased cosmological parameter estimation from galaxy surveys.

Findings

Successfully applied to synthetic data with physical structure formation.

Achieves at least 5 times more modes than state-of-the-art methods.

Provides unbiased, high-fidelity reconstruction of baryon acoustic oscillations.

Abstract

We propose a new, likelihood-free approach to inferring the primordial matter power spectrum and cosmological parameters from arbitrarily complex forward models of galaxy surveys where all relevant statistics can be determined from numerical simulations, i.e. black-boxes. Our approach, which we call simulator expansion for likelihood-free inference (SELFI), builds upon approximate Bayesian computation using a novel effective likelihood, and upon the linearisation of black-box models around an expansion point. Consequently, we obtain simple "filter equations" for an effective posterior of the primordial power spectrum, and a straightforward scheme for cosmological parameter inference. We demonstrate that the workload is computationally tractable, fixed a priori, and perfectly parallel. As a proof of concept, we apply our framework to a realistic synthetic galaxy survey, with a data model accounting for physical structure formation and incomplete and noisy galaxy observations. In doing so, we show that the use of non-linear numerical models allows the galaxy power spectrum to be safely fitted up to at least $k_\mathrm{max} = 0.5$ $h$/Mpc, outperforming state-of-the-art backward-modelling techniques by a factor of $\sim 5$ in the number of modes used. The result is an unbiased inference of the primordial matter power spectrum across the entire range of scales considered, including a high-fidelity reconstruction of baryon acoustic oscillations. It translates into an unbiased and robust inference of cosmological parameters. Our results pave the path towards easy applications of likelihood-free simulation-based inference in cosmology. We have made our code pySELFI and our data products publicly available at http://pyselfi.florent-leclercq.eu.

Tables

Table 1. Table 1: The statistical variables appearing in section II.1 and their interpretation in the context of galaxy survey data analysis.

Symbol	Meaning	Interpretation
$𝛉\in {ℝ}^S$	Target parameters	Parametrisation of the primordial matter power spectrum
$𝛙\in {ℝ}^T$	Nuisance parameters	Random numbers involved in the initial phase realisation, instrumental noise, etc.
$\text{d}\in {ℝ}^D$	Raw data	Galaxy number counts in a three-dimensional map of the survey volume
$𝚽\in {ℝ}^P$	Summary statistics of the data (observed or simulated)	Summaries of the galaxy number count field, such as its estimated power spectrum
${𝚽}_O\in {ℝ}^P$	Summary statistics of the observations	Summaries of the observed galaxy number count field
${𝚽}_{𝛉}\in {ℝ}^P$	Summary statistics of simulated data	Summaries of a galaxy number count field, simulated with primordial matter power spectrum given by $𝛉$
$\text{s}\in {ℝ}^P$	Virtual signal	True summaries of the galaxy number count field, if they were not degraded by nuisances

Table 2. Table 2: The statistical variables appearing in selfi and their interpretation in the context of galaxy survey data analysis.

Symbol	Meaning	Interpretation
${𝛉}_0\in {ℝ}^S$	Expansion point of the simulator in parameter space and prior mean	Fiducial primordial matter power spectrum
$\text{S}\in {ℝ}^{S\times S}$	Prior covariance matrix	Prior covariance matrix of the primordial matter power spectrum
${𝚽}_O\in {ℝ}^P$	Summary statistics of the observations	Summaries of the observed galaxy number count field
${\text{f}}_0\in {ℝ}^P$	Estimated average black-box at the expansion point	Estimated average of summaries of simulated galaxy fields with fiducial primordial matter power spectrum
${\text{C}}_0\in {ℝ}^{P\times P}$	Estimated covariance matrix of the black-box at the expansion point	Estimated covariance of summaries of simulated galaxy fields with fiducial primordial matter power spectrum
$\nabla {\text{f}}_0\in {ℝ}^{P\times S}$	Estimated gradient of the average black-box at the expansion point	Estimated gradient of summaries of simulated galaxy fields around the fiducial primordial matter power spectrum
$𝛄\in {ℝ}^S$	Posterior mean	Reconstructed primordial matter power spectrum, given the observed summaries ${𝚽}_O$
$𝚪\in {ℝ}^{S\times S}$	Posterior covariance matrix	Uncertainties on the reconstruction of the primordial matter power spectrum, given the observed summaries ${𝚽}_O$

Table 3. Table 3: The cosmological parameters used in this work.

$h$	${\Omega }_b$	${\Omega }_m$	$n_S$	${\sigma }_8$
0.6674	0.0486	0.3089	0.9667	0.8159

Equations

$$\text{mathpzc}P(\mathbf{\Phi }|\text{uptheta},\text{uppsi})={\text{updelta}}_{\mathrm{D}}(\mathbf{\Phi }-\text{mathpzc}B(\text{uptheta},\text{uppsi})).$$

$$\text{mathpzc}P(\text{uptheta}|\mathbf{\Phi }{)}_{|\mathbf{\Phi }={\mathbf{\Phi }}_{\mathrm{O}}}=L(\text{uptheta})\frac{\text{mathpzc}P(\text{uptheta})}{Z_{\mathbf{\Phi }}},$$

$$L(\text{uptheta})\equiv \text{mathpzc}P(\mathbf{\Phi }|\text{uptheta}{)}_{|\mathbf{\Phi }={\mathbf{\Phi }}_{\mathrm{O}}}$$

$$L(\text{uptheta})$$

$$\text{mathpzc}P(\text{uptheta}|\mathbf{\Phi },\{{\mathbf{\Phi }}_{\text{uptheta}}^{(i)}\})$$

$$\text{mathpzc}P(\text{uptheta}|\mathbf{\Phi }{)}_{|\mathbf{\Phi }={\mathbf{\Phi }}_{\mathrm{O}}}$$

$${\hat{L}}^N(\text{uptheta})\equiv \frac{\text{mathpzc}P(\mathbf{\Phi },\{{\mathbf{\Phi }}_{\text{uptheta}}^{(i)}\}|\text{uptheta}{)}_{|\mathbf{\Phi }={\mathbf{\Phi }}_{\mathrm{O}}}}{\text{mathpzc}P(\{{\mathbf{\Phi }}_{\text{uptheta}}^{(i)}\})}.$$

$$\text{mathpzc}P(\mathbf{\Phi },\{{\mathbf{\Phi }}_{\text{uptheta}}^{(i)}\}|\text{uptheta})$$

$$\text{mathpzc}P(\mathbf{\Phi },\{{\mathbf{\Phi }}_{\text{uptheta}}^{(i)}\}|\text{uptheta})\propto \int \text{mathpzc}P(\mathbf{\Phi }|\text{s})\left[\prod _{n=1}^N\text{mathpzc}P({\mathbf{\Phi }}_{\text{uptheta}}^{(i)}|\text{s})\right]\mathrm{d}\text{s}.$$

$$-2\log \text{mathpzc}P(\mathbf{\Phi }|\text{s})=\log \left|2\pi {\mathbf{\Sigma }}_{\text{uptheta}}\right|+(\mathbf{\Phi }-\text{s}{)}^{⊺}{\mathbf{\Sigma }}_{\text{uptheta}}^{-1}(\mathbf{\Phi }-\text{s}).$$

$$-2{\hat{\ell }}^N(\text{uptheta})=\log \left|2\pi {\hat{\mathbf{\Sigma }}}_{\text{uptheta}}^{\prime }\right|+({\mathbf{\Phi }}_{\mathrm{O}}-{\hat{\mathbf{\Phi }}}_{\text{uptheta}}{)}^{⊺}{\hat{\mathbf{\Sigma }}}_{\text{uptheta}}^{\prime -1}({\mathbf{\Phi }}_{\mathrm{O}}-{\hat{\mathbf{\Phi }}}_{\text{uptheta}}),$$

$${\hat{\mathbf{\Phi }}}_{\text{uptheta}}\equiv {\mathrm{E}}^N\left[{\mathbf{\Phi }}_{\text{uptheta}}\right]=\frac{1}{N}\sum _{i=1}^N{\mathbf{\Phi }}_{\text{uptheta}}^{(i)}$$

$${\hat{\mathbf{\Sigma }}}_{\text{uptheta}}^{\prime }\equiv \frac{N+1}{N}{\hat{\mathbf{\Sigma }}}_{\text{uptheta}},{\hat{\mathbf{\Sigma }}}_{\text{uptheta}}^{\prime -1}\equiv {\left(\frac{N+1}{N}\right)}^{-1}{\hat{\mathbf{\Sigma }}}_{\text{uptheta}}^{-1},$$

$${\mathbf{\Sigma }}_{\text{uptheta}}\approx {\hat{\mathbf{\Sigma }}}_{\text{uptheta}}$$

$${\mathbf{\Sigma }}_{\text{uptheta}}^{-1}\approx {\hat{\mathbf{\Sigma }}}_{\text{uptheta}}^{-1}\equiv \alpha {\left({\hat{\mathbf{\Sigma }}}_{\text{uptheta}}\right)}^{-1},\alpha \equiv \frac{N-P-2}{N-1},$$

$$-2\tilde{\ell }(\text{uptheta})\equiv \log \left|2\pi {\mathbf{\Sigma }}_{\text{uptheta}}\right|+({\mathbf{\Phi }}_{\mathrm{O}}-\text{s}{)}^{⊺}{\mathbf{\Sigma }}_{\text{uptheta}}^{-1}({\mathbf{\Phi }}_{\mathrm{O}}-\text{s}).$$

$${\hat{\mathbf{\Phi }}}_{\text{uptheta}}\approx {\text{f}}_0+\nabla {\text{f}}_0\cdot (\text{uptheta}-{\text{uptheta}}_0)\equiv \text{f}(\text{uptheta}),$$

$${\left(\nabla {\text{f}}_0\right)}_s^{⊺}\approx \frac{\text{f}({\text{uptheta}}_s)-{\text{f}}_0}{h},\mathrm{with}{\text{uptheta}}_s\equiv {\text{uptheta}}_0+h({\text{updelta}}_{\mathrm{K}}^{s{s}^{\prime }}{)}_{s^{\prime }\in [[1,S]]},$$

$$-2{\hat{\ell }}^N(\text{uptheta})\approx \log \left|2\pi {\text{C}}_0\right|+{\left[{\mathbf{\Phi }}_{\mathrm{O}}-\text{f}(\text{uptheta})\right]}^{⊺}{\text{C}}_0^{-1}\left[{\mathbf{\Phi }}_{\mathrm{O}}-\text{f}(\text{uptheta})\right].$$

$${\left(\text{K}\right)}_{s{s}^{\prime }}\equiv \exp \left[-\frac{1}{2}{\left(\frac{k_s-k_{s^{\prime }}}{k_{\mathrm{corr}}}\right)}^2\right].$$

$$(\text{u}{)}_s\equiv 1+{\sigma }_s=1+\frac{{\alpha }_{\mathrm{cv}}}{k_s^{3/2}}$$

$$\text{S}\equiv {\theta }_{\mathrm{norm}}^2\text{u}{\text{u}}^{⊺}\circ \text{K},$$

$$-2\log \text{mathpzc}P(\text{uptheta})\equiv \log \left|2\pi \text{S}\right|+(\text{uptheta}-{\text{uptheta}}_0{)}^{⊺}{\text{S}}^{-1}(\text{uptheta}-{\text{uptheta}}_0).$$

$$-2\log \text{mathpzc}P(\text{uptheta}|\mathbf{\Phi }{)}_{|\mathbf{\Phi }={\mathbf{\Phi }}_{\mathrm{O}}}=\log \left|2\pi \mathbf{\Gamma }\right|+(\text{uptheta}-\text{upgamma}{)}^{⊺}{\mathbf{\Gamma }}^{-1}(\text{uptheta}-\text{upgamma}),$$

$$\boxed{\text{upgamma}\equiv {\text{uptheta}}_0+\mathbf{\Gamma }(\nabla {\text{f}}_0{)}^{⊺}{\text{C}}_0^{-1}({\mathbf{\Phi }}_{\mathrm{O}}-{\text{f}}_0),}$$

$$\boxed{\mathbf{\Gamma }\equiv {\left[(\nabla {\text{f}}_0{)}^{⊺}{\text{C}}_0^{-1}\nabla {\text{f}}_0+{\text{S}}^{-1}\right]}^{-1}.}$$

$$\begin{array}{rl}-2\log \text{mathpzc}P(k_{\mathrm{corr}},{\theta }_{\mathrm{norm}}|{\text{uptheta}}_{\mathrm{fid}})\equiv & \log \left|2\pi \mathbf{\Gamma }\right|\\ & +({\text{uptheta}}_{\mathrm{fid}}-\text{upgamma}{)}^{⊺}{\mathbf{\Gamma }}^{-1}({\text{uptheta}}_{\mathrm{fid}}-\text{upgamma}),\end{array}$$

$$\text{mathpzc}P(\text{uptheta}|\text{upomega})={\text{updelta}}_{\mathrm{D}}\left(\text{uptheta}-\text{mathpzc}T(\text{upomega})\right),$$

$$\text{mathpzc}P(\text{upomega}|\mathbf{\Phi }{)}_{|\mathbf{\Phi }={\mathbf{\Phi }}_{\mathrm{O}}}={\hat{L}}_{\text{upomega}}^N(\text{upomega})\frac{\text{mathpzc}P(\text{upomega})}{Z_{\mathbf{\Phi }}},$$

$${\hat{L}}_{\text{upomega}}^N(\text{upomega})$$

Full text

I Introduction

The cosmic large-scale structure constitutes one of the major sources of information for modern cosmology. According to the current paradigm, all observable structures originate from tiny primordial fluctuations, which evolved via gravitational amplification into the presently-observed cosmic web (see e.g. Peebles, 1980; Peacock, 1999). The statistics of the initial density field are measured to be extremely close to Gaussian-distributed (see e.g. Planck Collaboration, 2018a). As a consequence, the primordial matter power spectrum is a very powerful cosmological probe: it is a sufficient statistical summary under the assumption that fluctuations are Gaussian, and – even if this assumption is violated – it remains close to capturing all of the information for all models allowed by observations. A particularly important cosmological signature, imprinted on the primordial matter power spectrum at the time of recombination, is baryon acoustic oscillations (BAOs). It constitutes a fixed comoving length scale (a “standard ruler”), which, when measured at different cosmic times, gives information on the expansion history of the Universe, including the late-time era of accelerated expansion. BAOs are thus one of the main probes to determine the equation of state of a possible dark energy component (see e.g. Eisenstein, 2005; Albrecht et al., 2006; Percival et al., 2007). A large variety of early-Universe models exhibit a deterministic relation between physical parameters of interest and the primordial matter power spectrum. The latter is therefore an interesting intermediate product for cosmological analyses, allowing parameter inference and model selection to be performed a posteriori without (or with minimal) loss of information. It can be seen as a largely agnostic and model-independent parametrisation of cosmological theories, which relies only on weak assumptions (isotropy and gaussianity).

For a long time, measuring the cosmological matter power spectrum has been one of the main goals of galaxy survey data analysis. However, inferring its shape accurately is a challenging task. Various systematic effects such as redshift uncertainties, complex survey geometries, selection effects, missing observations and foreground contamination can greatly hinder the measurement and analysis (see e.g. Ross et al., 2012; Jasche & Lavaux, 2017). This problem is particularly important for the next generation of optical surveys, such as provided by ESA’s Euclid space mission or the Large Synoptic Survey Telescope (LSST), which are expected to be dominated by systematic rather than statistical uncertainty (see e.g. Laureijs et al., 2011; LSST Science Collaboration, 2012). Even if systematic effects arising from the survey strategy were fully understood and controlled, many theoretical challenges would still be present: galaxy biasing, anisotropic clustering in redshift space, and non-linear structure growth at late times, which reduce the detectability of cosmological signatures such as BAOs (see e.g. Meiksin, White & Peacock, 1999; Eisenstein, Seo & White, 2007). Because of the limited reliability of data models, fits of the galaxy power spectrum focus on linear and mildly non-linear scales, using typically a largest wavenumber of $k_{\mathrm{max}}=0.3$ $h$/Mpc (e.g. Ross et al., 2015). However, the number of modes used in the analysis scales as $k_{\mathrm{max}}^{3}$, meaning that any improvement of data models at small scales (such as what can be achieved via numerical simulations instead of perturbation theory) gives access to much more cosmological information.

As a response to theoretical and observational challenges, many approaches to measure the power spectrum have been proposed. They can be divided into two broad categories: backward-modelling approaches (often associated with frequentist statistics, counting the frequencies of measurements in mock catalogues, and associated covariance matrices) and likelihood-based forward-modelling approaches (often associated with Bayesian statistics). Backward-modelling approaches suggest to directly account and correct for relevant effects in observational data as a pre-processing step. For example, the BAO “reconstruction” technique removes redshift-space distortions and corrects the density field from bulk motions using the inverse Zel’dovich approximation (e.g. Eisenstein et al., 2007; Padmanabhan et al., 2012; Doumler et al., 2013; Burden, Percival & Howlett, 2015; White, 2015). The end product is an estimator for the primordial matter power spectrum, which can be close to optimal if all relevant effects are modelled (Smith & Marian, 2015; Seljak et al., 2017). As backward-modelling approaches require substantial expert knowledge input, results are often largely model-dependent and with difficult propagation of uncertainties. Importantly, these approaches often rely on fiducial values for the parameters that are the target of the analysis. More recently, several thorough forward-modelling Bayesian approaches have been proposed to jointly infer the three-dimensional matter density field and its power spectrum from galaxy observations, while properly accounting for uncertainties and systematics (Jasche et al., 2010; Jasche & Wandelt, 2013b; Jasche & Lavaux, 2015; Granett et al., 2015; Jasche & Lavaux, 2017). Considerable effort is also put into reconstructing the primordial density field from present observations (Jasche & Wandelt, 2013a; Wang et al., 2013, 2014; Jasche, Leclercq & Wandelt, 2015; Lavaux & Jasche, 2016; Jasche & Lavaux, 2019; Bos, Kitaura & van de Weygaert, 2019). These methods are likelihood-based, meaning that they solve the exact inference problem by sampling from the target distribution via sophisticated Markov Chain Monte Carlo (MCMC) methods (Gibbs sampling and/or Hamiltonian Monte Carlo). In order to make the likelihood tractable, they have to involve approximations of the data model.

This paper has similar scientific aims but follows a different spirit: the presented method uses likelihood-free forward-modelling. It introduces one variant of approximate Bayesian computation (ABC) and treats the data model as a black-box simulator, i.e. performs inference without necessity to incorporate any knowledge of the data-generating processes into the analysis. This feature renders ABC ideal to infer the primordial matter power spectrum with arbitrarily complex models of galaxy surveys, including a physical treatment of structure formation and the details of observational processes, which cannot be trivially accounted for in likelihood-based statistical approaches. A popular ABC algorithm is likelihood-free rejection sampling, often coupled to Population Monte Carlo (e.g. Ishida et al., 2015; Akeret et al., 2015; Jennings & Madigan, 2017). More sophisticated approaches known as delfi (Alsing, Wandelt & Feeney, 2018) and bolfi (Leclercq, 2018) have also been recently introduced in cosmology. The ABC approach introduced in this work differs from all of the above in two aspects: (i) it allows the treatment of a much larger number of parameters (one hundred in this work), which correspond to primordial power spectrum amplitudes at different wavenumbers. To do so, (ii) it assumes the availability of a reasonably good guess of the target parameters, based on previous observations. This situation is fairly typical in cosmology, where the allowed space for parameters is already strongly constrained by previous experiments such as the Planck satellite (Planck Collaboration, 2016, 2018b). Under these assumptions, we derive an effective likelihood for the problem. We use an expansion point in parameter space and linearise the black-box around it. The use of finite differencing to compute the gradient of the black-box makes evaluations of the effective likelihood computationally feasible. When further assuming that the prior is Gaussian, we find that the effective posterior distribution for the primordial power spectrum is a Gaussian with mean and covariance matrix given by two simple “filter equations” (equations (25) and (26)), which constitute the main result of this work. Finally, we show how to infer parameters of specific cosmological models using the linearised black-box. We propose to call the algorithm “simulator expansion for likelihood-free inference” (selfi).

In order to illustrate the performance of the method, we apply it to a black-box which emulates realistic cosmological data. This black-box is built using Simbelmynë (Leclercq, Jasche & Wandelt, 2015), a hierarchical probabilistic simulator to generate synthetic galaxy survey data.111Simbelmynë is publicly available at https://bitbucket.org/florent-leclercq/simbelmyne. The data model involves a full cosmological $N$-body simulation (performed using our implementation of cola, Tassev, Zaldarriaga & Eisenstein, 2013) to evolve the three-dimensional initial density field. It includes a treatment of galaxy bias, redshift-space distortions, survey geometry, selection effects, and instrumental noise. The statistical summary chosen is the estimated power spectrum of the galaxy number count field, as is standard in large-scale structure data analysis, but can be readily extended to include more information. As a result, the inferred primordial matter power spectrum is unbiased across the entire range of Fourier modes considered, and includes in particular BAOs, which were not included in the expansion point. Our analysis demonstrates that by using a fully numerical data model in conjunction with our statistical approach, one can safely fit the galaxy power spectrum even far in the non-linear regime, up to at least $k_{\mathrm{max}}\approx 0.5$ $h$/Mpc, which provides a factor of $\sim 5$ increase in the number of modes used, with respect to state-of-the-art backward-modelling techniques. We stress that any possible refinement of the data model used in this work does not change the statistical method, and therefore does not affect the validity of the previous statement.

This paper is organised as follows. In section II, we discuss the statistical method and derive the equations for black-box simulation-based inference of the primordial matter power spectrum and cosmological parameters. In section III, we describe the data-generating model used to test our method. The results obtained by combining the two are discussed in IV. We discuss the application of our method and prospects for cosmological data analysis, and provide our conclusions in section V. Details of the statistical derivations are given in the appendices.

II Method

This section describes our method for simulation-based inference of the primordial matter power spectrum and cosmological parameters from black-box galaxy surveys. In section II.1, we design an effective likelihood for black-box models. In section II.2, we exploit previous knowledge, as could have been obtained by earlier cosmological probes, in order to linearise the black-box. We discuss the parametrisation of the primordial matter power spectrum and its prior distribution in section II.3. The equations for the effective posterior distribution are given in section II.4. In section II.5, we describe how to optimally choose the hyperparameters appearing in the prior. The inference of cosmological parameters from the linearised black-box is discussed in section II.6.

II.1 Design of an effective likelihood for black-box models

Table 1 provides an overview of the different variables appearing in this section and their interpretation in the context of galaxy survey data analysis.

II.1.1 The data model

We assume given a black-box model that provides realistic predictions for artificial observations when provided with all necessary input parameters. These consist of the target vector $\boldsymbol{\uptheta}\in\mathbb{R}^{S}$ parametrising the primordial power spectrum, and of nuisance parameters $\boldsymbol{\uppsi}\in\mathbb{R}^{T}$, independent of $\boldsymbol{\uptheta}$. Nuisance parameters account for the entire stochasticity of the data model, such as initial phases, noise realisations, sample variance, etc. In our case, nuisance parameters will be all the random numbers generated by the galaxy survey simulator for initial conditions and instrumental noise; there are typically $\mathcal{O}(10^{7})$ of those. Once realisations of $\boldsymbol{\uptheta}$ and $\boldsymbol{\uppsi}$ are specified, the output of the simulation $\textbf{d}\in\mathbb{R}^{D}$ is a deterministic numerical function $\mathpzc{S}$, i.e. $\mathpzc{P}(\textbf{d}|\boldsymbol{\uptheta},\boldsymbol{\uppsi})=\updelta_{\mathrm{D}}(\textbf{d}-\mathpzc{S}(\boldsymbol{\uptheta},\boldsymbol{\uppsi}))$, where the symbol $\mathpzc{P}$ denotes a probability distribution function (pdf) and $\updelta_{\mathrm{D}}$ a Dirac delta distribution. We refer to such realisations as mock observations. As usual in ABC approaches, the (often high-dimensional) raw prediction d can be compressed to a set of summary statistics $\boldsymbol{\Phi}\in\mathbb{R}^{P}$. We assume that this compression is a deterministic function $\mathpzc{C}$ of d, i.e. $\mathpzc{P}(\boldsymbol{\Phi}|\textbf{d})=\updelta_{\mathrm{D}}(\boldsymbol{\Phi}-\mathpzc{C}(\textbf{d}))$. It can be included in the model, so that the black-box is $\mathpzc{B}\equiv\mathpzc{C}\circ\mathpzc{S}$ and

[TABLE]

A graphical representation of the Bayesian hierarchical data model is presented in figure 1.

II.1.2 The exact Bayesian problem

Denoting by $\boldsymbol{\Phi}_{\mathrm{O}}$ the summary statistics of the observations, the inference problem considered is

[TABLE]

where the likelihood is

[TABLE]

and the normalisation constant is $Z_{\boldsymbol{\Phi}}\equiv\mathpzc{P}(\boldsymbol{\Phi})_{|\boldsymbol{\Phi}=\boldsymbol{\Phi}_{\mathrm{O}}}$. By marginalising over $\boldsymbol{\uppsi}$ and using equation (1), we have

[TABLE]

From equation (4), it is clear that the likelihood involves an intractable integral, the computation of which would require exactly hitting the observed summaries $\boldsymbol{\Phi}_{\mathrm{O}}$ with the black-box. We are therefore not able to explicitly formulate the true likelihood distribution for the considered problem.

II.1.3 The effective likelihood

To overcome this difficulty, in this section we derive an effective likelihood that allows us to perform inference by requiring black-box model evaluations only. For every $\boldsymbol{\uptheta}$, we can generate an ensemble of $N$ mock data realisations $\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}^{(i)}=\mathpzc{B}(\boldsymbol{\uptheta},\boldsymbol{\uppsi}^{(i)})$ with $i\in\{1,...,N\}$, by drawing independent and identically-distributed realisations of nuisance parameters $\boldsymbol{\uppsi}$ from the probability distribution $\mathpzc{P}(\boldsymbol{\uppsi})$ and evaluating the black-box. The approach then consists in explicitly conditioning all probabilities on $\{\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}^{(i)}\}$. Using Bayes’ theorem, we have:

[TABLE]

Using for $\boldsymbol{\Phi}$ the observed data $\boldsymbol{\Phi}_{\mathrm{O}}$, we thus have the new inference problem

[TABLE]

where we have defined the first factor on the right-hand side of equation (5) evaluated at $\boldsymbol{\Phi}_{\mathrm{O}}$ to be the effective likelihood (a computable approximation of the true likelihood):

[TABLE]

To arrive at a more explicit expression for the effective likelihood, we assume that observed data $\boldsymbol{\Phi}_{\mathrm{O}}$ and mock realisations $\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}^{(i)}$ are drawn from a common (but unknown) virtual signal $\textbf{s}\in\mathbb{R}^{P}$. The hierarchical representation of this effective data model is presented in figure 2. We assume that s carries the deterministic information on the target parameters $\boldsymbol{\uptheta}$, which is common to the data and to the mock observations; therefore, $\boldsymbol{\Phi}_{\mathrm{O}}$ and $\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}^{(i)}$ only differ by stochastic uncertainties described by the nuisance parameters, which carry no information of interest. Intuitively, s represents the “true” version of the summaries $\boldsymbol{\Phi}$, which is degraded by nuisances $\boldsymbol{\uppsi}$. Under this assumption, we introduce the probability distribution $\mathpzc{P}(\boldsymbol{\Phi}|\textbf{s})$, from which data realisations are drawn independently once the virtual signal s is given. Since we do not want to infer s explicitly, marginalisation yields:

[TABLE]

In some models, it may be possible to derive s from the model parameters $\boldsymbol{\uptheta}$. However, this is generally not true for numerical simulators, where the expected summaries have to be estimated through averaging mock realisations. In this work, in absence of prior information on s, we set $\mathpzc{P}(\textbf{s}|\boldsymbol{\uptheta})$ to a constant. The joint distribution of mock observations for a given virtual signal factorises, so that one can write

[TABLE]

In order to marginalise over s, we need to postulate a parametric form for the pdf $\mathpzc{P}(\boldsymbol{\Phi}|\textbf{s})$. In this work, we assume a Gaussian distribution centred on the virtual signal s with a covariance matrix $\boldsymbol{\Sigma}_{\boldsymbol{\uptheta}}$ quantifying the stochastic uncertainties inherent to the observations,222The investigation of different choices for $\mathpzc{P}(\boldsymbol{\Phi}|\textbf{s})$ is left to future investigations.

[TABLE]

The virtual signal can therefore be interpreted as the expectation $\mathrm{E}$ of $\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}$ once $\boldsymbol{\uptheta}$ is specified: $\textbf{s}=\mathrm{E}\left[\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}\right]$. We also have $\boldsymbol{\Sigma}_{\boldsymbol{\uptheta}}=\mathrm{E}\left[(\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}-\textbf{s})(\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}-\textbf{s})^{\intercal}\right]$.

Under these assumptions, marginalisation over s gives the effective likelihood as $\widehat{L}^{N}(\boldsymbol{\uptheta})=\exp\left[\hat{\ell}^{N}(\boldsymbol{\uptheta})\right]$ (the details of the computation are given in appendix A), with:

[TABLE]

where

[TABLE]

is the ensemble mean of mock observations and $\mathrm{E}^{N}$ stands for the empirical average over the set. The covariance matrix of $\widehat{L}^{N}(\boldsymbol{\uptheta})$ and its inverse are defined by

[TABLE]

where given a sufficiently large number of mock observations, a computable approximation of $\boldsymbol{\Sigma}_{\boldsymbol{\uptheta}}$ is estimated in the following way:

[TABLE]

The effective likelihood also requires an estimator $\boldsymbol{\hat{\Sigma}}_{\boldsymbol{\uptheta}}^{-1}$ of the inverse covariance matrix $\boldsymbol{\Sigma}_{\boldsymbol{\uptheta}}^{-1}$ (see equation (13)). As argued by Sellentin & Heavens (2016) and Jeffrey & Abdalla (2018), the proper Bayesian treatment would in fact consist in replacing the Gaussian likelihood by an alternative, corrected distribution. Nevertheless, keeping a Gaussian effective likelihood is an essential requirement of the present ABC technique; we will therefore be content with the Hartlap, Simon & Schneider (2007) correction, which consists of replacing the true inverse covariance matrix by a scaled inverse sample covariance matrix:

[TABLE]

where $P$ is the number of summary statistics. This correction debiases the expectation value $\left\langle\left(\boldsymbol{\hat{\Sigma}}_{\boldsymbol{\uptheta}}\right)^{-1}\right\rangle=\alpha^{-1}\boldsymbol{\Sigma}_{\boldsymbol{\uptheta}}^{-1}$ of the estimator, under the assumption that $\boldsymbol{\Sigma}_{\boldsymbol{\uptheta}}^{-1}$ is inverse-Wishart distributed.

The limiting approximation of $\widehat{L}^{N}(\boldsymbol{\uptheta})$ when $N\rightarrow\infty$ is $\widetilde{L}(\boldsymbol{\uptheta})=\mathpzc{P}(\boldsymbol{\Phi}|\textbf{s})_{|\boldsymbol{\Phi}=\boldsymbol{\Phi}_{\mathrm{O}}}=\exp\left[\tilde{\ell}(\boldsymbol{\uptheta})\right]$, with (as intended)

[TABLE]

There are a number of interesting similarities and differences between the effective likelihood (equation (11)) and the virtual signal pdf (equation (16)). First, both are Gaussian distributions with respect to the observed data $\boldsymbol{\Phi}_{\mathrm{O}}$. Second, the mean of the effective likelihood is given by the empirical average of simulated summaries, $\boldsymbol{\hat{\Phi}}_{\boldsymbol{\uptheta}}=\mathrm{E}^{N}\left[\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}\right]$. This is an unbiased estimator of the virtual signal $\textbf{s}=\mathrm{E}\left[\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}\right]$ for a sufficiently large number of simulations. Importantly, the replacement of the expectation $\mathrm{E}$ by an empirical average $\mathrm{E}^{N}$ was not an assumption (contrary to the synthetic likelihood, Wood, 2010; Price et al., 2018), but naturally appeared in the derivation (see appendix A). Finally, the covariance of the effective likelihood is $\frac{N+1}{N}\boldsymbol{\hat{\Sigma}}_{\boldsymbol{\uptheta}}$, a multiple of the estimated covariance $\boldsymbol{\hat{\Sigma}}_{\boldsymbol{\uptheta}}$. The numerical prefactor, similar to Bessel’s correction for the empirical sample variance, can be understood as follows. For a small number of simulations, the latent space associated to the virtual signal increases the observed scatter. For example, for $N=1$, the covariance to be used in equation (11) is $2\boldsymbol{\hat{\Sigma}}_{\boldsymbol{\uptheta}}$, which reflects the fact that observed and simulated data $\boldsymbol{\Phi}_{\mathrm{O}}$ and $\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}$ can be drawn from opposite ends of the scatter around s of $\mathpzc{P}(\boldsymbol{\Phi}|\textbf{s})$. However, when $N\rightarrow\infty$, $\frac{N+1}{N}\boldsymbol{\hat{\Sigma}}_{\boldsymbol{\uptheta}}\longrightarrow\boldsymbol{\Sigma}_{\boldsymbol{\uptheta}}$, the intrinsic covariance of $\mathpzc{P}(\boldsymbol{\Phi}|\textbf{s})$. This result is reasonable, since for large sample sizes the ensemble mean $\boldsymbol{\hat{\Phi}}_{\boldsymbol{\uptheta}}$ converges to s, the mean of $\mathpzc{P}(\boldsymbol{\Phi}|\textbf{s})$; in the same limit, we expect the covariance of the effective likelihood to approach $\boldsymbol{\Sigma}_{\boldsymbol{\uptheta}}$, the covariance of $\mathpzc{P}(\boldsymbol{\Phi}|\textbf{s})$.

Note that the scaling factor $\frac{N+1}{N}$ directly arises from the presence of the virtual signal, as shown in appendix A. Alternatively, we could have directly assumed a parametric form for the pdf $\mathpzc{P}(\boldsymbol{\Phi}|\boldsymbol{\uptheta})$ without introducing the latent variable s. For instance, the synthetic likelihood directly amounts to postulating equation (11), but without the scaling factor $\frac{N+1}{N}$ for the estimated covariance matrix. This means that the end result for $\hat{\ell}^{N}(\boldsymbol{\uptheta})$ has little sensitivity to our treatment and assumption for $\mathpzc{P}(\textbf{s}|\boldsymbol{\uptheta})$. Furthermore, since $\frac{N+1}{N}\longrightarrow 1$ when $N\rightarrow\infty$, our result with the virtual signal and the synthetic likelihood are equivalent, provided that the number of simulations is large enough.

Determining $\boldsymbol{\hat{\Phi}}_{\boldsymbol{\uptheta}}$ and $\boldsymbol{\hat{\Sigma}}_{\boldsymbol{\uptheta}}$ requires $N$ model evaluations per target parameters $\boldsymbol{\uptheta}$. In the next section, we show how the evaluation of this effective likelihood can be made more efficient when exploiting prior information on $\boldsymbol{\uptheta}$.

II.2 Linearisation of black-box models

The numerical cost of the computation of the effective likelihood described in section II.1 may be prohibitively large when the full parameter space has to be explored. However, such an extensive exploration is not always required: often, one has sufficient information to be only interested in a small region of parameter space around a specific prediction. This requirement is fulfilled for inferences of the primordial power spectrum from galaxy surveys. The target parameters $\boldsymbol{\uptheta}$ typically consist of power spectrum amplitudes in about a hundred different bands of wavevectors; thus the corresponding parameter space is very large. However, the CMB already provides exquisite measurements of the primordial cosmological power spectrum (e.g. Planck Collaboration, 2018a). Any large deviations from these previous measurements would most likely be rejected on methodological grounds (sample size, uncontrolled systematics, etc.) and not be attributed to new physics. It therefore seems reasonable to focus the inference of the primordial power spectrum from galaxy surveys within a narrow region in parameter space around a previous estimate obtained from CMB results.

Following this reasoning, we focus on searching for solutions corresponding to small deviations $\Delta\boldsymbol{\uptheta}$ around an expansion point $\boldsymbol{\uptheta}_{0}$. Consequently, the target parameters are given by $\boldsymbol{\uptheta}=\boldsymbol{\uptheta}_{0}+\Delta\boldsymbol{\uptheta}$. Assuming that $\boldsymbol{\hat{\Phi}}_{\boldsymbol{\uptheta}}$ is differentiable with respect to $\boldsymbol{\uptheta}$, we perform a first-order Taylor expansion in $\Delta\boldsymbol{\uptheta}$ around $\boldsymbol{\uptheta}_{0}$:

[TABLE]

where the defined function f is a linearised version of the averaged black-box. The first term corresponds to the mean mock observations at the expansion point $\boldsymbol{\uptheta}_{0}$, i.e. $\textbf{f}_{0}\equiv\boldsymbol{\hat{\Phi}}_{\boldsymbol{\uptheta}_{0}}$, while the second term involves $\nabla\textbf{f}_{0}$, a $P\times S$ matrix corresponding to the gradient of mean mock observations at $\boldsymbol{\uptheta}_{0}$, whose components are $\left(\nabla\textbf{f}_{0}\right)_{ps}\equiv\frac{\partial\boldsymbol{\hat{\Phi}}_{\boldsymbol{\uptheta}_{0}p}}{\partial\boldsymbol{\uptheta}_{s}}$. In this work, we estimate the gradient $\nabla\textbf{f}_{0}$ via finite differencing: given a small step size $h$, each column is approximated by

[TABLE]

where $\updelta_{\mathrm{K}}$ is the Kronecker delta. To avoid obtaining a noisy gradient, nuisance parameters $\boldsymbol{\uppsi}$ are kept at the same values in the simulations used to compute the $\textbf{f}(\boldsymbol{\uptheta}_{s})$ and $\textbf{f}_{0}$. Further, we neglect the dependence of the covariance matrix on $\boldsymbol{\uptheta}$ and we estimate it at the expansion point according to equation (14), i.e. $\boldsymbol{\hat{\Sigma}}_{\boldsymbol{\uptheta}}^{\prime}\approx\boldsymbol{\hat{\Sigma}}_{\boldsymbol{\uptheta}_{0}}^{\prime}\equiv\textbf{C}_{0}$. Importantly, fully characterising f under these assumptions only requires evaluations of the simulator, thus ensuring that the data model remains a black-box.

Using the linearised data model f as a proxy for the data model described in section II.1.3 simplifies the expression of the effective likelihood (equation (11)) to

[TABLE]

Evaluating $\textbf{f}_{0}$ and $\textbf{C}_{0}$ requires $N_{0}$ model evaluations at the expansion point $\boldsymbol{\uptheta}_{0}$. The computation of the gradient further requires $N_{s}\times S$ model evaluations. Therefore, the linearised data model is fully characterised by a fixed total of $N_{0}+N_{s}\times S$ model evaluations. $N_{0}$ and $N_{s}$ are a user choice, but should be of the order of the dimensionality of the data space $P$. Particularly at the expansion point $\boldsymbol{\uptheta}_{0}$, a minimum would be $N_{0}\geq P+3$ but it can be worth investing more simulations, in order for the estimated covariance matrix $\textbf{C}_{0}$ and its inverse to be precise. Since all nuisance parameters are kept fixed in the computation of the $\textbf{f}(\boldsymbol{\uptheta}_{s})$, $N_{s}$ can in principle be smaller than $P$ (and as small as $1$); however $N_{s}\gtrsim P$ will yield a safer evaluation of $\nabla\textbf{f}_{0}$.

It is important to note that all required model evaluations are done once and for all. Once the linearised data model is known, it is not necessary to perform additional black-box evaluations in order to perform inference from new data. This feature makes the present approach similar to supervised machine learning algorithms, which are only trained once before being applied to multiple data sets. Furthermore, all data model evaluations can be done in parallel, or even on different machines, making the approach very suitable for grid computing.

II.3 The power spectrum prior distribution

In this paper, we aim at inferring the primordial matter power spectrum $P(k)$, which is a continuous function of wavenumber $k$. We parametrise it by its amplitudes at a sufficient number $S$ of support wavenumbers $k_{s}$. As we are particularly interested in BAOs, we fix a “wiggle-less” power spectrum $P_{0}(k)$ and we work with the “wiggle function” $\theta(k)\equiv P(k)/P_{0}(k)$ as target function. Formally, with $\textbf{P}_{0}\in\mathbb{R}^{S}$ the vector of components $(\textbf{P}_{0})_{s}\equiv P_{0}(k_{s})$, the inference variable $\boldsymbol{\uptheta}$ is defined as the $S$-dimensional vector of components $P(k_{s})/(\textbf{P}_{0})_{s}$.

In order to set up the Bayesian problem, we have to formulate the prior probability distribution of $\boldsymbol{\uptheta}$. In this work, we include into our prior the following assumptions:

the power spectrum is Gaussian-distributed, 2. 2.

it is strongly constrained to live close to $P_{0}$, 3. 3.

it is a smooth function of wavenumber, 4. 4.

and the power spectrum $P_{0}$ is subject to cosmic variance.

It follows from assumptions 1 and 2 that $\mathpzc{P}(\boldsymbol{\uptheta})$ shall be a Gaussian distribution with mean $\boldsymbol{\uptheta}_{0}\equiv\boldsymbol{1}_{\mathbb{R}^{S}}$, which is also to be used as the expansion point.

We now discuss assumptions 3 and 4, in order to build the covariance matrix S of $\mathpzc{P}(\boldsymbol{\uptheta})$. Let us define the matrix K as a radial basis function, i.e. by its coefficients

[TABLE]

The hyperparameter $k_{\mathrm{corr}}$ determines the length scale on which power spectrum amplitudes of different wavenumber correlate with each other. A large value of $k_{\mathrm{corr}}$ corresponds to a strong correlation between $\boldsymbol{\uptheta}_{s}$ and $\boldsymbol{\uptheta}_{s^{\prime}}$, even if their corresponding wavenumbers are far from each other. On the other hand, for $k_{\mathrm{corr}}\ll k_{s},k_{s^{\prime}}$, K becomes the identity matrix and $\boldsymbol{\uptheta}_{s}$ and $\boldsymbol{\uptheta}_{s^{\prime}}$ do not correlate with each other at all. The previous discussion implies that the a priori smoothness of the wiggle function can be changed by tuning $k_{\mathrm{corr}}$. In a realistic scenario, the prior covariance is not scale-independent as is K ($(\textbf{K})_{ss}=1$ for all $s$). We rather want the standard deviations $(\textbf{S})^{1/2}_{ss}$ to account for the cosmic variance affecting the power spectrum $P_{0}$ (and hence the mean $\boldsymbol{\uptheta}_{0}$). In terms of power spectrum amplitudes, cosmic variance at a scale $k$ is given by $P(k)^{2}/N_{k}$, where $N_{k}\propto k^{3}$ is the number of modes of wavenumber $k$ in the considered cosmological volume. Thus, in order to account for cosmic variance in terms of $\boldsymbol{\uptheta}$, the coefficients K shall be multiplied by the coefficients of $\textbf{u}\textbf{u}^{\intercal}$, where

[TABLE]

and $\alpha_{\mathrm{cv}}$ is a hyperparameter characterising the “strength” of cosmic variance given the considered volume. Finally, the amplitude of the covariance matrix S can be captured by an overall scaling $\theta_{\mathrm{norm}}^{2}$. The final expression for the prior covariance matrix is therefore

[TABLE]

where $\circ$ is the Hadamard product. The standard deviations on the diagonal are $(\textbf{S})^{1/2}_{ss}=\theta_{\mathrm{norm}}(1+\sigma_{s})$, as intended.

The resulting prior on $\boldsymbol{\uptheta}$ is characterised by a set of three hyperparameters $\{k_{\mathrm{corr}},\alpha_{\mathrm{cv}},\theta_{\mathrm{norm}}\}$, and given as

[TABLE]

As demonstrated in section IV, a prior of this form results in a smooth posterior mean for the primordial matter power spectrum, while incorporating reasonable uncertainties around the expansion point.

II.4 The power spectrum effective posterior distribution

Using the effective likelihood with the linearised black-box data model, given in equation (19), and the prior given in equation (23), we arrive at the final expression for the effective posterior distribution (see equation (6)). It is a Gaussian distribution,

[TABLE]

with mean

[TABLE]

and covariance matrix

[TABLE]

The proof of this result (detailed in appendix B) uses the same algebra as the derivation of the Wiener filter (Wiener, 1964). Equations (25) and (26) are the main result of this work. Equation (25) provides a simple “filter equation” to infer the primordial matter power spectrum from galaxy observations via complex black-box simulations. Corresponding uncertainties are quantified by the covariance matrix given in equation (26). Table 2 summarises the variables appearing in equations (25) and (26) and their interpretation in the context of galaxy survey data analysis.

II.5 Optimisation of the prior hyperparameters

As discussed in section II.3, the chosen prior distribution involves three hyperparameters $\{k_{\mathrm{corr}},\alpha_{\mathrm{cv}},\theta_{\mathrm{norm}}\}$. By definition, $\alpha_{\mathrm{cv}}$ characterises the strength of cosmic variance in the considered cosmological volume, such that the number of modes at a given scale $k$ is $N_{k}=k^{3}/\alpha_{\mathrm{cv}}^{2}$. $\alpha_{\mathrm{cv}}$ can therefore simply be measured for a given simulator setup (i.e. box size and mesh). On the contrary, $k_{\mathrm{corr}}$ and $\theta_{\mathrm{norm}}$ are free hyperparameters and the power spectrum reconstruction ($\boldsymbol{\upgamma}$ and $\boldsymbol{\Gamma}$) generally depends on their values. In this section, we propose a procedure to find optimal values for $k_{\mathrm{corr}}$ and $\theta_{\mathrm{norm}}$.

$k_{\mathrm{corr}}$ indicates the a priori smoothness of reconstructed wiggle functions $\theta(k)$ and $\theta_{\mathrm{norm}}$ how much they can deviate from the expansion point. Together, these two parameters characterise the functional shapes of allowed target functions (much like kernels and their hyperparameters in techniques such as Gaussian process regression, see e.g. Rasmussen & Williams, 2006). In analogy with hyperparameter optimisation in machine learning, we propose to optimise $k_{\mathrm{corr}}$ and $\theta_{\mathrm{norm}}$ to reproduce the shape of a fiducial wiggle function $\boldsymbol{\uptheta}_{\mathrm{fid}}$, using as likelihood the effective posterior distribution derived in section II.4. More precisely, the likelihood for $k_{\mathrm{corr}}$ and $\theta_{\mathrm{norm}}$ is defined by

[TABLE]

where $\boldsymbol{\upgamma}$ and $\boldsymbol{\Gamma}$, defined by equations (25) and (26), are functions of $k_{\mathrm{corr}}$ and $\theta_{\mathrm{norm}}$ (through S). At this point, it is of course possible to include a hyperprior on $(k_{\mathrm{corr}},\theta_{\mathrm{norm}})$, if desired. The maximum likelihood estimator (or maximum a posteriori estimator) then provides the optimal values of $k_{\mathrm{corr}}$ and $\theta_{\mathrm{norm}}$ to be used to infer $\boldsymbol{\uptheta}$, its functional shape being assumed to be that of $\boldsymbol{\uptheta}_{\mathrm{fid}}$.

Note that evaluating the likelihood given in equation (27) is cheap once the linearised black-box f has been computed: no additional data model evaluation is required, only low-dimensional matrix operations. This allows the optimal values of $k_{\mathrm{corr}}$ and $\theta_{\mathrm{norm}}$ to be found using standard optimisers.

II.6 From the power spectrum to cosmological parameters

As argued in the introduction, the primordial matter power spectrum can be seen as a largely model-independent parametrisation of the underlying theory. The goal of this section is to go from the power spectrum to parameters of specific cosmological models. This last step in the analysis can be seen as adding a layer to the Bayesian hierarchical model (see figure 3): $\boldsymbol{\upomega}$ is a vector of cosmological parameters which generates the primordial power spectrum coefficients $\boldsymbol{\uptheta}$. This generative process is usually deterministic, i.e.

[TABLE]

where $\mathpzc{T}$ is a deterministic function of cosmological parameters, typically a Boltzmann solver or a fitting function. Given this assumption, we have $\mathpzc{P}(\boldsymbol{\Phi}|\boldsymbol{\upomega})=\mathpzc{P}(\boldsymbol{\Phi}|\mathpzc{T}(\boldsymbol{\upomega}))$, and the inference of cosmological parameters gives

[TABLE]

where the likelihood for cosmological parameters is defined as

[TABLE]

Noting $\widehat{L}_{\boldsymbol{\upomega}}^{N}(\boldsymbol{\upomega})\equiv\exp\left[\hat{\ell}_{\boldsymbol{\upomega}}^{N}(\boldsymbol{\upomega})\right]$ and using the linearised data model, we have (see equation (19))

[TABLE]

Going from the primordial power spectrum to cosmological parameters is conceptually easy, but often hard in practice, when one only has samples of the power spectrum likelihood (or posterior). In this case, the likelihood to be used for cosmological parameter inference is naively represented by a sum of Dirac delta distributions. As it usually lives in a parameter space containing hundreds to thousands of dimensions, the number of samples is always a limiting factor. For this reason, techniques that effectively broaden the obtained samples, such as Blackwell-Rao estimators (Wandelt, Larson & Lakshminarayanan, 2004) or kernel density estimates have been introduced. We note that our ABC technique does not suffer from this complication: the likelihood for cosmological parameters (equation (31)) is a Gaussian centered on the data $\boldsymbol{\Phi}_{\mathrm{O}}$ with a fixed covariance matrix $\textbf{C}_{0}$. Furthermore, getting predictions for $\boldsymbol{\Phi}$ amounts to evaluating $\textbf{f}(\mathpzc{T}(\boldsymbol{\upomega}))$, which does not require additional black-box evaluations once f is known. The Bayesian problem of inferring cosmological parameters $\boldsymbol{\upomega}$ (equation (29)) can therefore easily be solved by standard Markov Chain Monte Carlo techniques.

III Data model

In this section, we describe the generation of mock observations used to test the performance of our method. We emphasise that the statistical method presented in section II is applicable to any black-box model, which can feature arbitrarily complex processes. Therefore, the details of the simulator presented in this section are of no relevance to the performance of the statistical method. It is only used in section IV as a showcase, to highlight the performance of our method to handle a complex simulator.

III.1 Primordial power spectrum parametrisation

Throughout this paper, we work with a cubic equidistant grid with comoving side length of $1~{}\mathrm{Gpc}/h$ and $256^{3}$ voxels, spanning scales between $k_{s,\mathrm{min}}=6.28\times 10^{-3}$ $h$/Mpc and $k_{s,\mathrm{max}}=1.4$ $h$/Mpc. We use $S=100$ support wavenumbers $k_{s}$. The first eight are fixed to the values required by the Fourier grid given our setup. The remaining support wavenumbers are logarithmically spaced up to $k_{\mathrm{max}}$. Any vector $\boldsymbol{\uptheta}$ in parameter space is defined at the scales of these support wavenumbers.

Between two consecutive support wavenumbers, we interpolate power spectra $P(k)$ using a one-dimensional spline fit, using $n=5$ as the degree of the smoothing spline. We checked that this setup yields vanishing differences in the representation of cosmological power spectra, at all wavenumbers of the Fourier grid used in this work.

III.2 Galaxy surveys

The data model used in this work is a non-linear process meant to approximate the large variety of physical and observational phenomena at play in galaxy surveys. To do so, it uses the Simbelmynë cosmological code, an end-to-end generative process for galaxy survey data given a specified primordial power spectrum $P(k)$. The flat $\Lambda$CDM model is assumed, and fiducial cosmological parameters used are the Planck 2015 values (Planck Collaboration, 2016, table 4, last column), given in table 3. These are used whenever the distance-redshift relation is needed, as well as for the gravitational evolution.333In principle, our approach would require treating cosmological parameters as nuisance parameters in the inference of the primordial power spectrum, and to marginalise over them. In this work, for simplicity, we keep them fixed to the values used to predict the expansion point.

Simbelmynë first generates a realisation of the initial density contrast via the convolution approach (see e.g. Peacock & Heavens, 1985). Specifically, the code generates a white noise field w, such that in each cell $x$ the value $w_{x}$ is drawn from the zero-mean unit-variance Gaussian distribution. The white noise field is multiplied by the square root of the desired cosmological power spectrum in Fourier space to give

[TABLE]

Transformation back to configuration space yields a initial density contrast field $\boldsymbol{\updelta}^{\mathrm{i}}$. One realisation of such an initial density field is shown in the left panel of figure 4.

Generated initial density fields then act as inputs to numerical structure formation simulations. The initial grid of $256^{3}$ voxels is populated by $512^{3}$ dark matter particles placed on a regular lattice. These particles are evolved to the redshift of $z=19$ via second order Lagrangian Perturbation theory (2LPT) (see e.g. Moutarde et al., 1991; Bouchet et al., 1995; Bouchet, 1996), then with an efficient implementation of cola (COmoving Lagrangian Acceleration, Tassev, Zaldarriaga & Eisenstein, 2013) from $z=19$ to $z=0$. A particle-mesh grid of $1024^{3}$ voxels and 20 timesteps linearly-spaced in the scale factor are used for the evolution with cola. We checked that, at the scales of interest, this setup yields negligible difference in the representation of final density fields with respect to the prediction of the fully non-linear code gadget-2 (Springel, 2005). In particular, final density fields contain the additional power expected from non-linear structure formation, at 1% precision up to $k=0.5$ $h/\mathrm{Mpc}$ and 5% up to $k=1$ $h/\mathrm{Mpc}$. We place the observer at the centre of the box. The maximal distance to the observer is $866$ Mpc/$h$, which we consider sufficiently small to neglect light-cone effects. Therefore, in this paper, we only use the final snapshot of our simulations at redshift zero, and ignore the evolution of matter within the survey volume.

Using their final peculiar velocities with respect to the observer $v_{r}$, dark matter particles are placed in redshift space according to the non-linear mapping

[TABLE]

where $z_{\mathrm{cosmo}}$ is the true cosmological redshift, $z_{\mathrm{obs}}$ is the “observed” redshift and $\mathrm{c}$ is the speed of light. Note that we do not work in the plane-parallel approximation. The particles are then binned to a $256^{3}$-voxels grid with the cloud-in-cell scheme (Hockney & Eastwood, 1981) to give the final density contrast field $\boldsymbol{\updelta}^{\mathrm{f}}$. One realisation of the final redshift-space density is shown in the middle panel of figure 4.

The galaxy density $\boldsymbol{\uprho}^{\mathrm{g}}$ is predicted using a linear bias model, used in various previous studies (e.g. Verde et al., 2002; Ross et al., 2015), such that in any cell $x$,

[TABLE]

In this work, we use $b=1.2$ and $\bar{N}=0.119$ (corresponding to an observed galaxy number density $\bar{n}=2\times 10^{-3}~{}(h/\mathrm{Mpc})^{3}$, achievable for instance with the Euclid spectroscopic survey, see e.g. Majerotto et al., 2012). For the sake of simplicity, $b$ and $\bar{N}$ are fixed, but they could be straightforwardly treated as additional nuisance parameters and marginalised over.

The last step corresponds to a virtual observation of the galaxy field, accounting for observational effects expected in actual surveys. To do so, we use the three-dimensional survey response operator (or window) W, consisting of the product of the radial selection function $R(r)$ and the angular survey mask and completeness function $C(\hat{\textbf{n}})$ for any line-of-sight $\hat{\textbf{n}}$, i.e. $\textbf{W}(\hat{\textbf{n}},r)\equiv R(r)\,C(\hat{\textbf{n}})$. This operator accounts for the fact that we are only looking at certain parts of the sky and that we have different detection probabilities of galaxies depending on their distance. We obtained a simple model for W on a grid matching our simulations ($256^{3}$ voxels covering a volume of $(1~{}\mathrm{Gpc}/h)^{3}$) as follows. For the angular completeness, we mask the galactic plane by setting $C(\hat{\textbf{n}})$ to [math] for galactic latitudes $-10\degree\leq b\leq 10\degree$ and $1$ otherwise, the system of coordinates being defined such that the observer is at the origin and the plane of equation $z=0$ is the galactic plane. For the radial selection function $R(r)$, we use a Schechter luminosity function (Schechter, 1976) with previously-published parameters for the $r$-band: $\alpha=-1.05$ and $M_{\ast}=-20.44$ (Blanton et al., 2003), a limiting apparent magnitude of $m=18.5$ and absolute magnitude cuts $-25\leq M\leq-21$. This choice makes our synthetic survey complete up to a luminosity distance of $D_{\mathrm{L}}=794$ Mpc/$h$, corresponding to $r=684$ Mpc/$h$ (see figure 5).

We emulate the survey via galaxy number counts $N_{x}^{\mathrm{g}}$ in each cell $x$ of the box. In order to account for instrumental noise, our model is a non-uniform Gaussian process, meaning that in every cell $N_{x}^{\mathrm{g}}$ is drawn from a Gaussian distribution with mean $\mu^{\mathrm{g}}_{x}$ and standard deviation $\sigma_{x}$,

[TABLE]

The mean $\mu^{\mathrm{g}}_{x}\equiv W_{x}\rho^{\mathrm{g}}_{x}=W_{x}\bar{N}(1+b\,\delta^{\mathrm{f}}_{x})$ characterises the expected number of galaxies. It accounts for all physical effects via $\delta^{\mathrm{f}}_{x}$ and for the survey response operator $W_{x}$. The standard deviation is defined as $\sigma_{x}\equiv\sigma\sqrt{W_{x}\bar{N}}$, where we set the overall noise level $\sigma\equiv 10^{-1}$. A slice through one realisation of the galaxy number count field $\textbf{N}^{\mathrm{g}}$ is shown in the right panel of figure 4.

III.3 Summary statistic

As is standard practice in cosmological data analysis, we do not work at the level of the entire galaxy number count map $\textbf{N}^{\mathrm{g}}$, but use a compression $\mathpzc{C}$ of the map to well-chosen statistical summaries. In this work, we limit ourselves to an estimator of the final power spectrum of the survey. It should be remarked however, that the method can use any other statistical summary and even combinations of those.

We obtain the binned data power spectrum by taking the squared modulus of the Fourier transform of $\textbf{N}^{\mathrm{g}}$, summing the contributions over all the Fourier cells within wavenumber shell $k_{r}$, and normalizing:

[TABLE]

The overall constant factor $C=1000^{3}/256^{2}~{}(\mathrm{Mpc}/h)^{3}$ arises from our Fourier transform convention; $N_{k_{r}}$ represents the number of modes within the shell $k_{r}$; and the factor $-2$ arises from the assumption that the data power spectrum is inverse-$\Gamma$ distributed with shape parameter $N_{k_{r}}/2$ and scale parameter $C\times\left|N^{\mathrm{g}}_{k}\right|^{2}/2$ (see Jasche et al., 2010).

As the galaxy data model relies on a full cosmological simulation, the predictions for the final power spectrum remain reasonable even at scales that have experienced substantial non-linearity. In spite of the approximations made, we trust our data model at the percent level up to $k=0.5$ $h$/Mpc. For this reason, we use $P=43$ $k_{r}$-bins in the range $[0.02,0.5]$ $h$/Mpc, ensuring that each bin contains at least $100$ modes. These bins are logarithmically spaced for $k_{r}\geq 0.04$ $h$/Mpc. For convenience, we normalise the output of the black-box using the expansion point, so that $\boldsymbol{\Phi}$ is the $P$-dimensional vector of components $A\times P^{\mathrm{f}}(k_{r})/P_{0}(k_{r})$ with $A=50$. The estimation of $P^{\mathrm{f}}(k)$ is completely deterministic once the galaxy number counts are given. Therefore, we now have a complete model to generate artificial realisations of $\boldsymbol{\Phi}$ for a given primordial power spectrum $P(k)$ and specific realisations of initial phases and noise.

Due to the breakdown of models based on perturbation theory in the non-linear regime and to the difficulties in incorporating the impact of small-scale observational processes, state-of-the-art large-scale structure analyses are typically limited to $k_{\mathrm{max}}\lesssim 0.3$ $h$/Mpc (e.g. Ross et al., 2015). Pushing the analysis to $k_{\mathrm{max}}=0.5$ $h$/Mpc represents an increase by a factor of $\sim 5$ in the number of modes used (scaling as $k_{\mathrm{max}}^{3}$), which is expected to yield substantial improvements in the inference results.

III.4 Idealised data model

For testing purposes, we also define an idealised data model corresponding to a Gaussian random field. More specifically, using exactly the same setup as before, the black-box here simply consists of producing the initial density field $\boldsymbol{\updelta}_{\mathrm{i}}$ (see section III.2), scaling it to redshift zero using the linear growth factor, and measuring its normalised power spectrum $\boldsymbol{\Phi}$ (see section III.3). The use of the two different black-boxes (Gaussian random field and realistic mock survey) within our method will quantify the effect of non-linear gravity, redshift-space distortions, and survey complications on primordial power spectrum inference, in particular the detectability of BAOs.

IV Results

This section describes the results obtained by applying the statistical method proposed in section II in conjunction with the data-generating process described in section III to an artificial galaxy survey, itself generated using the same process.

We use for $P_{0}(k)$ the “wiggle-less” BBKS power spectrum (Bardeen et al., 1986) under Planck 2015 cosmology (see table 3). Unknown ground truth cosmological parameters $\boldsymbol{\upomega}_{\mathrm{gt}}$ are drawn from the (marginalised, Gaussian) Planck priors:

[TABLE]

The “wiggly” ground truth power spectrum $P_{\textrm{gt}}(k)$ is generated with the Eisenstein & Hu (1998) (EH) fitting function, using these cosmological parameters. It is used to simulate observed data $\boldsymbol{\Phi}_{\mathrm{O}}$, with unknown nuisance parameters (phase realisation and instrumental noise). For later use, the fiducial “wiggly” power spectrum $P_{\mathrm{fid}}(k)$ is also generated with the EH prescription, using Planck cosmology. The target parameters $(\boldsymbol{\uptheta})_{s}\equiv P(k_{s})/P_{0}(k_{s})$ are the values of the wiggle function at the $S=100$ support wavenumbers defined in section III.1. We note $\boldsymbol{\uptheta}_{\mathrm{gt}}$ and $\boldsymbol{\uptheta}_{\mathrm{fid}}$ the vectors of component $P_{\mathrm{gt}}(k_{s})/P_{0}(k_{s})$ and $P_{\mathrm{fid}}(k_{s})/P_{0}(k_{s})$, respectively.

IV.1 Diagnostics of the black-box

We created an ensemble of $N_{0}=150$ mock realisations at the expansion point $\boldsymbol{\uptheta}_{0}=\boldsymbol{1}_{\mathbb{R}^{S}}$ using different nuisance parameters. These are used to compute the average black-box at the expansion point, $\textbf{f}_{0}\equiv\boldsymbol{\hat{\Phi}}_{\boldsymbol{\uptheta}_{0}}$, and the covariance matrix $\textbf{C}_{0}\equiv\boldsymbol{\hat{\Sigma}}^{\prime}_{\boldsymbol{\uptheta}_{0}}$, using their definitions (equations (12), (13), and (14)). The results are shown in figure 6. There, the left panel shows individual realisations $\boldsymbol{\Phi}_{\boldsymbol{\uptheta}_{0}}$ and the observed data vector $\boldsymbol{\Phi}_{\mathrm{O}}$. The average black-box $\textbf{f}_{0}$ is also plotted, with a credible region corresponding to two standard deviations (i.e. $2\sqrt{\mathrm{diag}(\textbf{C}_{0}})$). The full estimated covariance matrix $\textbf{C}_{0}$ is shown in the right panel. As expected, the measured variance is larger on large scales due to cosmic variance, with some anti-correlations between pairs of bins. The effect of the mask, which increases power at the largest scales found in the simulation box, is also clearly visible.

Using a step size of $h=10^{-2}$ and an ensemble of $N_{s}=100$ mock realisations at each of the expansion points $\boldsymbol{\uptheta}_{s}$, we measured the gradients of the black-box $(\nabla\textbf{f}_{0})^{\intercal}_{s}$ along all directions of parameter space (see equation (18)). The nuisance parameters (phase realisation and noise) are kept at fixed values (the ones corresponding to the first $N_{s}$ realisations generated at the expansion point) for this calculation. The results are shown in figure 7, where the left panel shows $(\nabla\textbf{f}_{0})^{\intercal}_{s}$ for individual values of $s$ and the right panel shows the full rectangular matrix $\nabla\textbf{f}_{0}$. Some interesting phenomena can be observed. At large scales, (see e.g. for $k_{20}=0.0364$ $h$/Mpc) exciting one initial mode ony triggers an answer in the bins closest to this scale; the gradient therefore resembles a multiple of the identity function. This is the result expected from linear perturbation theory. However, at small scales, the non-linear simulator couples modes. This implies that the response is smaller in amplitude but distributed over a much larger ranges of scales (see e.g. for $k_{70}=0.3780$ $h$/Mpc). In the non-linear regime, the gradient is typically negative at large scales, crosses zero slightly before the excited scale $k_{s}$, then becomes positive at smaller scales.

As discussed in section II.2, the linearised black-box $\textbf{f}(\boldsymbol{\uptheta})$ is fully characterised by $\textbf{f}_{0}$, $\textbf{C}_{0}$ and $\nabla\textbf{f}_{0}$. In this work, we used a total of $N_{0}+N_{s}\times S=10,150$ simulations to get very precise estimates of $\textbf{C}_{0}$ and $\nabla\textbf{f}_{0}$, although using fewer would have been possible.

IV.2 The prior and its optimisation

As discussed in section II.3, we choose a Gaussian prior centered on the expansion point $\boldsymbol{\uptheta}_{0}$ with a covariance matrix S given by equation (22) and characterised by three hyperparameters $\{k_{\mathrm{corr}},\alpha_{\mathrm{cv}},\theta_{\mathrm{norm}}\}$.

Following the method presented in section II.5, we found optimal values for the prior hyperparameters. The strength of cosmic variance within our simulation volume shall satisfy $\alpha_{\mathrm{cv}}=\sqrt{k^{3}/N_{k}}$ at all scales $k$, where $N_{k}$ is the number of modes. In our Fourier grid (described in section III.1), we measured up to the Nyquist frequency $\alpha_{\mathrm{cv}}=8.848\times 10^{-4}$, value that we adopt.

Assuming that the target function $\boldsymbol{\uptheta}$ follows the functional shape of the fiducial wiggle function $\boldsymbol{\uptheta}_{\mathrm{fid}}$ calculated with Planck cosmology, the likelihood for $k_{\mathrm{corr}}$ and $\theta_{\mathrm{norm}}$ is given by equation (27). We further assume broad, uncorrelated Gaussian hyperpriors on $k_{\mathrm{corr}}$ and $\theta_{\mathrm{norm}}$: $k_{\mathrm{corr}}\sim\mathpzc{G}(0.020,0.015^{2})$ [$h$/Mpc] and $\theta_{\mathrm{norm}}\sim\mathpzc{G}(0.2,0.3^{2})$. The posterior surface is plotted in figure 8. We found the maximum a posteriori values using the popular optimiser L-BFGS (Byrd et al., 1995). In our run, these are $k_{\mathrm{corr}}=0.0156$ $h$/Mpc and $\theta_{\mathrm{norm}}=0.0478$. The resulting prior covariance matrix is shown in figure 9. For the idealised data model (see section III.4), the same procedure is applied and the optimal parameters are found to be $k_{\mathrm{corr}}=0.0158$ $h$/Mpc and $\theta_{\mathrm{norm}}=0.0535$.

IV.3 The effective posterior

Using the effective likelihood, characterised by the observed data $\boldsymbol{\Phi}_{\mathrm{O}}$ and the linearised black-box described in section IV.1, as well as the optimised prior discussed in section IV.2, we obtained the effective posterior on $\boldsymbol{\uptheta}$. It is a Gaussian, with mean and covariance matrix given by the “filter equations” (25) and (26).

Figure 10 shows the inferred primordial wiggle function $\boldsymbol{\upgamma}$ in comparison with the expansion point $\boldsymbol{\uptheta}_{0}$, the ground truth $\boldsymbol{\uptheta}_{\mathrm{gt}}$, and the fiducial function used to optimise prior hyperparameters, $\boldsymbol{\uptheta}_{\mathrm{fid}}$. $2\sigma$ credible regions are shown for the prior and the posterior (i.e. $2\sqrt{\mathrm{diag}(\textbf{S})}$ and $2\sqrt{\mathrm{diag}(\boldsymbol{\Gamma})}$, respectively). The top panel corresponds to the result obtained using the idealised data model (a Gaussian random field, see section III.4), and the bottom panel to the result obtained using the realistic mock survey data model. In both cases, the inference is unbiased since the ground truth always lies within the $2\sigma$ credible intervals of the reconstruction. As discussed in the introduction, this effective posterior can be seen as a largely model-independent parametrisation of the theory, containing all the available cosmological information under weak assumptions.

As can be read from the figure, the inferred vector contains the BAO wiggles, even far within the Silk damping tail (Silk, 1968). All visible oscillations are fully reconstructed in the idealised case, and in the realistic case, up to scales of $k\approx 0.3$ $h$/Mpc. In particular, 5 acoustic oscillations are unambiguously identified, which is competitive with the latest cosmic microwave background experiments and has been so far out of reach of galaxy surveys. Note that this result is obtained given a simulation volume of $(1~{}\mathrm{Gpc}/h)^{3}$ and that further improvements could be obtained with a larger volume, as will be probed by upcoming surveys. The inferred wiggle function shows higher uncertainty in regions of small and large wavenumbers. This is due to cosmic variance and noise, respectively. Cosmic variance reflects the limited number of modes that we have at the largest scales in our simulation volume. This effect limits the significance of the determination of the cosmological power spectrum at these scales. For this reason, at $k\lesssim 0.05$ $h$/Mpc the reconstruction is driven towards the prior mean, which is also the expansion point and the default answer in the absence of informative data. As expected, this effect is visible in the idealised as well as in the realistic case. On the other hand, noise (understood as the combined effect of the specific phase realisation of the data, non-linear gravity, redshift-space distortions and instrumental noise) acts on smaller scales ($0.2~{}h/\mathrm{Mpc}\lesssim k\lesssim 0.5~{}h/\mathrm{Mpc}$). Some of the primordial information is effectively destroyed at these scales – or at least is not captured by the statistical summaries $\boldsymbol{\Phi}$. The reconstruction is therefore also driven towards the prior mean, and the uncertainty is increased, because the data are less informative than in the idealised case. As expected, we recover the prior at $k\gtrsim 0.5$ $h$/Mpc (in fact a little below, due to mode coupling), since the data $\boldsymbol{\Phi}$ do not contain measurements at these scales.

It is important to note that the prior $\mathpzc{P}(\boldsymbol{\uptheta})$, used in this section to regularise the inference of the primordial power spectrum, does not appear in the inference of cosmological parameters in the next section (only the effective likelihood does). Thus, no bias is introduced in cosmological parameter inference when the amplitude of reconstructed BAO wiggles seems to undershoot the ground truth.

Since the proposed method is fully Bayesian, we do not simply obtain a point estimate, but a complete probability distribution, which provides a detailed quantification of uncertainties. Figure 11 shows the covariance matrix $\boldsymbol{\Gamma}$ of the Gaussian effective posterior; it can be compared with the prior covariance matrix S shown in figure 9.

IV.4 Cosmological parameters

The last step in the analysis is to infer parameters $\boldsymbol{\upomega}$ of specific cosmological models, given the observed data $\boldsymbol{\Phi}_{\mathrm{O}}$. In this section, we assume a flat $\Lambda$CDM model, characterised by $5$ parameters $\{h,\Omega_{\mathrm{b}},\Omega_{\mathrm{m}},n_{\mathrm{S}},\sigma_{8}\}$. For simplicity, our prescription to generate primordial power spectra given cosmological parameters is the EH fitting function. The generative process $\mathpzc{T}$ (from $\boldsymbol{\upomega}$ to $\boldsymbol{\uptheta}$) is therefore the EH fitting function divided by $P_{0}$, sampled at the support wavenumbers $k_{s}$. It could easily be generalised to extensions of the flat $\Lambda$CDM model and to include a Boltzmann solver. The linearised data model f (from $\boldsymbol{\uptheta}$ to $\boldsymbol{\Phi}$) has already been characterised in section IV.1.

The effective likelihood $\widehat{L}_{\boldsymbol{\upomega}}^{N}(\boldsymbol{\upomega})$ for cosmological parameters is given by equation (31). Consistently with the expansion point used to linearise the black-box, we complement $\widehat{L}_{\boldsymbol{\upomega}}^{N}(\boldsymbol{\upomega})$ with a Gaussian prior $\mathpzc{P}(\boldsymbol{\upomega})$ centred on Planck cosmological parameters, but with broader variance: the diagonal covariance matrix given in equation (37) is multiplied by a factor of $3$. We explored the effective posterior $\mathpzc{P}(\boldsymbol{\upomega}|\boldsymbol{\Phi})_{|\boldsymbol{\Phi}=\boldsymbol{\Phi}_{\mathrm{O}}}$ (equation (29)) via MCMC (performed using the emcee code, Foreman-Mackey et al., 2013), ensuring sufficient convergence.

Results are shown in figure 12. Prior contours are shown in blue, and posterior contours are shown in red and purple for two different realisations of the data $\boldsymbol{\Phi}_{\mathrm{O}}$. The two data realisations have been generated using different ground truth cosmological parameters $\boldsymbol{\upomega}_{\mathrm{gt}}$ (shown as dashed and dotted lines, respectively), as well as different nuisance parameters (phase realisation and noise). The plot demonstrates that cosmological parameter inference is unbiased and robust to nuisances imprinted in the data.

V Discussion and conclusion

The biggest challenge in galaxy survey analyses arises from the requirement of non-linear data modelling. In this work, we described the development of a novel simulation-based Bayesian approach, selfi, which can be used to infer the primordial matter power spectrum and cosmological parameters from galaxy surveys. The main results are the “filter equations” (25) and (26). They can be applied to get an effective posterior for any model where the mean and covariance of the data are estimated from arbitrarily complex forward models. Essentially everything is obtained from a simulator, which can be treated as a black-box, without necessity to include any knowledge of its internal mechanisms into the statistical analysis.

We derived the “filter equations” under two assumptions: the availability of a black-box able to generate artificial data, and of strong prior constraints in parameter space, obtained from a previous experiment. We built an effective likelihood for this scenario and made its evaluation efficient by linearising the black-box around an expansion point. We devised a method to optimise the hyperparameters appearing in our power spectrum prior. Finally, we derived the cheap likelihood for parameters of specific cosmological models, to be used in our framework.

The approach presented in this paper relies on likelihood-free forward-modelling via ABC. It complements statistically exact, likelihood-based forward-modelling techniques. The principal differences are as follows.

• •

First, the numerical complexity of likelihood-based MCMC approaches typically requires to approximate complex data models to allow for fast execution speeds. In this work, we rather aimed at performing approximate inference, but with full-scale black-box models. This approach allows a much more accurate modelling of cosmological data, including in particular the complicated systematics they experience.

• •

Second, for MCMC methods, the number of data model evaluations is not fixed a priori, as some proposed samples are rejected during runtime. One has to assess the convergence of the chain. In contrast, our method only requires a fixed number of realisations to characterise the effective likelihood with the linearised black-box, all of which are used to obtain the inference result. For $S$ target parameters and $N_{0}$ nuisance parameters realisations at the expansion point, $N_{s}$ nuisance parameters realisations along each direction in parameter space, one has to perform $N_{0}+N_{s}\times S$ data model evaluations. Increasing $N_{0}$ or $N_{s}$ to get a better estimate of the required covariance matrix or gradient only increases the overall computational cost linearly.

• •

Third, while MCMC have to be computed sequentially, all the simulations required in the proposed method can be run simultaneously in parallel, or even on different machines. This allows a fast application of the method and makes it particularly suitable for grid computing.

• •

Finally, the linearised black-box is trained once and for all independently of the data. This means that if one acquires new data from the same survey, no additional black-box evaluations are required to perform inference, whereas likelihood-based techniques would require a new MCMC. Furthermore, if the cosmological simulations used are stored, they can even be used to perform inference from a different survey, by just replacing the part of the black-box corresponding to survey specifications.

selfi also differs from other approaches to ABC (such as likelihood-free rejection sampling, Population Monte Carlo, delfi, or bolfi), which are limited by their inability to scale with the number of target parameters. By relying on an expansion of the simulator, selfi allows the likelihood-free inference of $S\gtrsim 100$ parameters, as is necessary for a model-independent parametrisation of theory in cosmology.

In this work, we demonstrated that a “non-wiggly” expansion point $\boldsymbol{\uptheta}_{0}$ is sufficient to recover the target wiggle function $\boldsymbol{\uptheta}$ in the domain allowed by Planck priors. However, it shall be noted that if the solution is farther from the expansion point, then the method can be iterated. In this case, the posterior mean $\boldsymbol{\upgamma}$ would be used as the new expansion point to train a new linearised black-box, used to obtain a new posterior. Using a sufficient number of iterations, we expect the effective posterior to converge to the true function, even if it strongly deviates from the first expansion point. We leave the detailed investigation of this idea to future studies.

In this paper, we showed a successful application of finite differencing to obtain the gradient of the averaged black-box $\nabla\textbf{f}_{0}$, but our equations could be used with other techniques, such as automatic differentiation. The data covariance matrix at the expansion point $\textbf{C}_{0}$ also needs to be evaluated; for this task, and for certain summary statistics, variance reduction techniques such as the use of fixed and paired simulations (Angulo & Pontzen, 2016; Villaescusa-Navarro et al., 2018) or hybrid estimators (Hall & Taylor, 2019) could further be exploited.

As a proof of concept, we applied our technique in conjunction with the artificial galaxy survey simulator Simbelmynë, emulating relevant effects at play: non-linear gravitational structure formation, redshift-space distortions, a survey mask and selection function, and instrumental noise. As a result, the inferred primordial power spectrum is unbiased with a distinct identification of BAO wiggles, even far in the Silk damping tail. We also demonstrated that unbiased inference of cosmological parameters is possible. In spite of the non-linear evolution of structures on small scales, we are able to use the power spectrum of the galaxy field as summary statistic, up to $k_{\mathrm{max}}=0.5$ $h$/Mpc. This represents an increase by a factor of $\sim 5$ in the number of modes used with respect to state-of-the-art perturbation theory and backward-modelling techniques, with perspectives for further improvements. Assuming that posteriors are Gaussian and modes are independent, this increase translates into a reduction of the size of credible contours by a factor $\sim\sqrt{5}$.

The data model used in this work remains simplified with respect to some of the complications found in real galaxy surveys. However, due to the flexible nature of the method, it is straightforward to include additional aspects in the inference process: one only has to exchange the black-box for a more sophisticated one. We developed a python code reflecting this versatility, pyselfi, which we publicly released, together with documentation and the data necessary to reproduce the results of the present paper.444Currently, the code’s homepage is hosted at http://pyselfi.florent-leclercq.eu; the sources are available on GitHub at https://github.com/florent-leclercq/pyselfi; and the documentation is on Read the Docs at https://pyselfi.readthedocs.io. The application of this method to more complex models and to real survey data is left for future research.

In conclusion, the method developed constitutes a computationally efficient and easily applicable framework to infer the primordial matter power spectrum and cosmological parameters from complex black-box mock observations. It allows the use of fully non-linear data models, as required for an optimal analysis of galaxy surveys. Other applications may include the cosmic microwave background, weak gravitational lensing, or the 21 cm signal of hydrogen. The prize for using full forward-modelling in these problems is a potentially vast gain of precision in cosmological constraints.

Appendix A Derivation of the effective likelihood

In this appendix, we derive the approximate likelihood given in equation (11), starting from equation (9) and the assumptions detailed in section II.1.3. We have

[TABLE]

Using the parametric form for $\mathpzc{P}(\boldsymbol{\Phi}|\textbf{s})$ (equation (10)) yields $\mathpzc{P}(\boldsymbol{\Phi},\{\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}^{(i)}\}|\boldsymbol{\uptheta})\propto\int\exp\left[\hat{\ell}_{1}(\boldsymbol{\Phi},\{\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}^{(i)}\},\textbf{s})\right]\,\mathrm{d}\textbf{s}$, with

[TABLE]

In order to evaluate the integral, we complete the square with respect to s in the argument of the exponential,

[TABLE]

where we have recognised $\boldsymbol{\hat{\Phi}}_{\boldsymbol{\uptheta}}=\frac{1}{N}\sum_{i=1}^{N}\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}^{(i)}$ (equation (12)) and introduced $\boldsymbol{\upeta}\equiv(\boldsymbol{\Phi}+N\boldsymbol{\hat{\Phi}}_{\boldsymbol{\uptheta}})/(N+1)$. After integration over s, the last term gives a constant factor, so that $\mathpzc{P}(\boldsymbol{\Phi},\{\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}^{(i)}\}|\boldsymbol{\uptheta})\propto\exp\left[\hat{\ell}_{2}(\boldsymbol{\Phi},\{\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}^{(i)}\})\right]$, with

[TABLE]

We now complete the square with respect to $\boldsymbol{\Phi}$ to obtain

[TABLE]

In order to obtain a computable approximation of the likelihood, the unknown covariance $\boldsymbol{\Sigma}_{\boldsymbol{\uptheta}}$ in $\hat{\ell}_{2}(\boldsymbol{\Phi},\{\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}^{(i)}\})$ has to be approximated by $\boldsymbol{\hat{\Sigma}}_{\boldsymbol{\uptheta}}$, defined by equation (14). The covariance of the effective likelihood is therefore $\boldsymbol{\hat{\Sigma}}_{\boldsymbol{\uptheta}}^{\prime}\equiv\frac{N+1}{N}\boldsymbol{\hat{\Sigma}}_{\boldsymbol{\uptheta}}$. The unknown inverse covariance $\boldsymbol{\Sigma}_{\boldsymbol{\uptheta}}^{-1}$ is also replaced by its unbiased computable approximation $\boldsymbol{\hat{\Sigma}}_{\boldsymbol{\uptheta}}^{-1}$, defined by equation (15). Finally, we use the normalisation condition $\int\mathpzc{P}(\boldsymbol{\Phi}|\{\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}^{(i)}\},\boldsymbol{\uptheta})\,\mathrm{d}\boldsymbol{\Phi}=\int\frac{\mathpzc{P}(\boldsymbol{\Phi},\{\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}^{(i)}\}|\boldsymbol{\uptheta})}{\mathpzc{P}(\{\boldsymbol{\Phi}_{\boldsymbol{\uptheta}}^{(i)}\}|\boldsymbol{\uptheta})}\,\mathrm{d}\boldsymbol{\Phi}=1$ and evaluate at $\boldsymbol{\Phi}=\boldsymbol{\Phi}_{\mathrm{O}}$, as prescribed by equation (7), to obtain $\hat{\ell}^{N}(\boldsymbol{\uptheta})$ given by equation (11). When $N\rightarrow\infty$, $\boldsymbol{\hat{\Phi}}_{\boldsymbol{\uptheta}}\longrightarrow\textbf{s}$ and $\boldsymbol{\hat{\Sigma}}_{\boldsymbol{\uptheta}}^{\prime}\longrightarrow\boldsymbol{\Sigma}_{\boldsymbol{\uptheta}}$, thus the limiting approximation is $\tilde{\ell}(\boldsymbol{\uptheta})$ given by equation (16).

Appendix B Derivation of the effective posterior

We recall the canonical form of the Gaussian distribution with mean $\textbf{x}_{0}$ and covariance matrix X, given as

[TABLE]

where $\boldsymbol{\upxi}_{0}\equiv\textbf{X}^{-1}\textbf{x}_{0}$.

Using the linearised data model (equation (17)) in the expression of the effective likelihood (equation (19)), we get

[TABLE]

where we have defined

[TABLE]

and

[TABLE]

$(\nabla\textbf{f}_{0})^{-1}$ denoting the adjoint of the Jacobian characterising the linearised black-box (its computation will not be necessary). In canonical form, the Gaussian effective likelihood is written

[TABLE]

with $\boldsymbol{\upmu}_{0}\equiv\textbf{N}_{0}^{-1}\textbf{y}_{0}=\textbf{N}_{0}^{-1}\boldsymbol{\uptheta}_{0}+(\nabla\textbf{f}_{0})^{\intercal}\textbf{C}_{0}^{-1}(\boldsymbol{\Phi}_{\mathrm{O}}-\textbf{f}_{0})$. Similarly, the prior (equation (23)) is written

[TABLE]

where $\boldsymbol{\upeta}_{0}\equiv\textbf{S}^{-1}\boldsymbol{\uptheta}_{0}$.

Adding the two expressions, we find that the effective posterior verifies

[TABLE]

This is the canonical form of a Gaussian distribution, where the covariance matrix is identified as $\boldsymbol{\Gamma}\equiv(\textbf{N}_{0}^{-1}+\textbf{S}^{-1})^{-1}$, giving equation (26), and the mean is identified as

[TABLE]

giving equation (25).

Note that the above calculation is analogous to the derivation of the Wiener filter equations: assuming a linear data model ($\textbf{d}=\textbf{s}+\textbf{n}$), a prior with mean $\bar{\textbf{s}}$ and signal covariance S, and a likelihood with mean $\bar{\textbf{d}}$ and noise covariance N, the filter covariance is $(\textbf{N}^{-1}+\textbf{S}^{-1})^{-1}$ and the filtered signal is $\bar{\textbf{s}}+(\textbf{N}^{-1}+\textbf{S}^{-1})^{-1}\textbf{N}^{-1}(\textbf{d}_{\mathrm{O}}-\bar{\textbf{d}})$.

Statement of contribution

Study concept and design (JJ, FL, WE); design of prior optimisation and cosmological parameter inference (FL); original code implementation of the filter equations (WE); code rewriting and enhancements (FL); design and implementation of the data model (FL); running of the simulations (FL); drafting of the manuscript (WE); critical revision of the manuscript (FL); proofreading (FL, WE, JJ, AH); supervision (JJ, FL); support and interpretation of results (JJ, AH). All authors read and approved the final manuscript.

Acknowledgements.

FL is grateful to Guilhem Lavaux and Andrew Jaffe for useful discussions. This work has made use of a modified version of pygtc (Bocquet & Carter, 2016). Numerical computations were done on the cx1 cluster hosted by the Research Computing Service facilities at Imperial College London (doi:10.14469/hpc/2232). This work is done within the Aquila Consortium (https://aquila-consortium.org). FL acknowledges funding from the Imperial College London Research Fellowship Scheme. This research was supported by the DFG cluster of excellence “Origin and Structure of the Universe” (www.universe-cluster.de).