# Primordial power spectrum and cosmology from black-box galaxy surveys

Florent Leclercq, Wolfgang Enzi, Jens Jasche, Alan Heavens

► To cite this version:

Florent Leclercq, Wolfgang Enzi, Jens Jasche, Alan Heavens. Primordial power spectrum and cosmology from black-box galaxy surveys. Monthly Notices of the Royal Astronomical Society, 2019, 490 (3), pp.4237-4253. 10.1093/mnras/stz2718. hal-05347497

HAL Id: hal-05347497

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

Submitted on 6 Nov 2025

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

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

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

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

# Primordial power spectrum and cosmology from black-box galaxy surveys

Florent Leclercq<sup>1,★</sup> Wolfgang Enzi,<sup>2</sup> Jens Jasche<sup>3,4</sup> and Alan Heavens<sup>1</sup>

<sup>1</sup>Imperial Centre for Inference and Cosmology (ICIC) & Astrophysics Group, Imperial College London, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK

<sup>2</sup>Max-Planck Institute for Astrophysics, Karl-Schwarzschild Strasse 1, D-85748 Garching, Germany

<sup>3</sup>The Oskar Klein Centre, Department of Physics, Stockholm University, Albanova University Center, SE-106 91 Stockholm, Sweden

<sup>4</sup>Excellence Cluster Universe, Technische Universität München, Boltzmannstrasse 2, D-85748 Garching, Germany

Accepted 2019 September 24. Received 2019 August 7; in original form 2019 March 1

## 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 linearization 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_{max} = 0.5 \ h \ \text{Mpc}^{-1}$, outperforming state-of-the-art backward-modelling techniques by a factor of ~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.

**Key words:** methods: statistical – cosmological parameters – large-scale structure of Universe.

## 1 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 I 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

<footnote>* E-mail: florent.leclercq@polytechnique.org</footnote>
selection to be performed a posteriori without (or with minimal)
loss of information. It can be seen as a largely agnostic and model-
independent parametrization 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
2007a). 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_{max} = 0.3\ h\ Mpc^{-1}$
(e.g. Ross et al. 2015). However, the number of modes used in
the analysis scales as $k_{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 associ-
ated 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.
2007b; 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, sev-
eral 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; Granett et al. 2015; Jasche & Lavaux
2015, 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; Bos, Kitaura & van de Weygaert
2019; Jasche & Lavaux 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. Akeret et al. 2015; Ishida 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 (100 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 XIII 2016; Planck Collaboration I 2018b). Under these assumptions, we derive an effective likelihood for the problem. We use an expansion point in parameter space and linearize 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 linearized 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 that 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.<fnref n="1" /> The data model involves a full cosmological *N*-body
simulation [performed using our implementation of COMOVING
LAGRANGIAN ACCELERATION (COLA); Tassev, Zaldarriaga & Eisen-
stein 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

<footnote n="1">¹SIMBELMYNË is publicly available at https://bitbucket.org/florent-leclercq/simbelmyne.</footnote>
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_{max} \approx 0.5 h \text{ Mpc}^{-1}$, which provides a factor of ~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 organized as follows. In Section 2, 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 3, we describe the data-generating model used to test our method. The results obtained by combining the two are discussed in Section 4. We discuss the application of our method and prospects for cosmological data analysis, and provide our conclusions in Section 5. Details of the statistical derivations are given in the appendices.

## 2 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 2.1, we design an effective likelihood for black-box models. In Section 2.2, we exploit previous knowledge, as could have been obtained by earlier cosmological probes, in order to linearize the black box. We discuss the parametrization of the primordial matter power spectrum and its prior distribution in Section 2.3. The equations for the effective posterior distribution are given in Section 2.4. In Section 2.5, we describe how to optimally choose the hyperparameters appearing in the prior. The inference of cosmological parameters from the linearized black box is discussed in Section 2.6.

### 2.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.

#### 2.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 $\theta \in \mathbb{R}^S$ parametrizing the primordial power spectrum, and of nuisance parameters $\psi \in \mathbb{R}^T$, independent of $\theta$. Nuisance parameters account for the entire stochasticity of the data model, such as initial phases, noise realizations, 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 $O(10^7)$ of those. Once realizations of $\theta$ and $\psi$ are specified, the output of the simulation $d \in \mathbb{R}^D$ is a deterministic numerical function $S$, i.e. $\mathcal{P}(d|\theta, \psi) = \delta_D(d - S(\theta, \psi))$, where the symbol $\mathcal{P}$ denotes a probability distribution function (pdf) and $\delta_D$ a Dirac delta distribution. We refer to such realizations as mock observations. As usual in ABC approaches, the (often high-dimensional) raw prediction $d$ can be compressed to a set of summary statistics $\Phi \in \mathbb{R}^P$. We assume that this compression is a deterministic function $C$ of $d$, i.e. $\mathcal{P}(\Phi|d) = \delta_D(\Phi - C(d))$. It can be included in the model, so that the black box is $\mathcal{B} \equiv C \circ S$ and

$$ \mathcal{P}(\Phi|\theta, \psi) = \delta_D(\Phi - \mathcal{B}(\theta, \psi)). \quad (1) $$

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

#### 2.1.2 The exact Bayesian problem

Denoting by $\Phi_O$ the summary statistics of the observations, the inference problem considered is

$$ \mathcal{P}(\theta|\Phi)_{|\Phi=\Phi_O} = L(\theta) \frac{\mathcal{P}(\theta)}{Z_{\Phi}}, \quad (2) $$

where the likelihood is

$$ L(\theta) \equiv \mathcal{P}(\Phi|\theta)_{|\Phi=\Phi_O} \quad (3) $$

and the normalization constant is $Z_{\Phi} \equiv \mathcal{P}(\Phi)_{|\Phi=\Phi_O}$. By marginalizing over $\psi$ and using equation (1), we have

$$ \begin{aligned} L(\theta) &= \int \mathcal{P}(\Phi|\theta, \psi)_{|\Phi=\Phi_O} \mathcal{P}(\psi) d\psi \\ &= \int \delta_D(\Phi_O - \mathcal{B}(\theta, \psi)) \mathcal{P}(\psi) d\psi. \end{aligned} \quad (4) $$

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

#### 2.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 $\theta$, we can generate an ensemble of $N$ mock data realizations $\Phi_{\theta}^{(i)} = \mathcal{B}(\theta, \psi^{(i)})$ with $i \in \{1, \dots, N\}$, by drawing independent and identically distributed realizations of nuisance parameters $\psi$ from the probability distribution $\mathcal{P}(\psi)$ and evaluating the black box. The approach then consists in explicitly conditioning all probabilities on $\{\Phi_{\theta}^{(i)}\}$. Using Bayes' theorem, we have

$$ \begin{aligned} \mathcal{P}(\theta|\Phi, \{\Phi_{\theta}^{(i)}\}) &= \frac{\mathcal{P}(\Phi, \{\Phi_{\theta}^{(i)}\}|\theta) \mathcal{P}(\theta)}{\mathcal{P}(\Phi, \{\Phi_{\theta}^{(i)}\})} \\ &= \frac{\mathcal{P}(\Phi, \{\Phi_{\theta}^{(i)}\}|\theta) \mathcal{P}(\theta)}{\mathcal{P}(\{\Phi_{\theta}^{(i)}\}) \mathcal{P}(\Phi)}. \end{aligned} \quad (5) $$

Using for $\Phi$ the observed data $\Phi_O$, we thus have the new inference problem

$$ \begin{aligned} \mathcal{P}(\theta|\Phi)_{|\Phi=\Phi_O} &\approx \mathcal{P}(\theta|\Phi, \{\Phi_{\theta}^{(i)}\})_{|\Phi=\Phi_O} \\ &= \hat{L}^N(\theta) \frac{\mathcal{P}(\theta)}{Z_{\Phi}}, \end{aligned} \quad (6) $$

where we have defined the first factor on the right-hand side of equation (5) evaluated at $\Phi_O$ to be the effective likelihood (a computable approximation of the true likelihood):

$$ \hat{L}^N(\theta) \equiv \frac{\mathcal{P}(\Phi, \{\Phi_{\theta}^{(i)}\}|\theta)_{|\Phi=\Phi_O}}{\mathcal{P}(\{\Phi_{\theta}^{(i)}\})} \quad (7) $$

To arrive at a more explicit expression for the effective likelihood, we assume that observed data $\Phi_O$ and mock realizations $\Phi_{\theta}^{(i)}$ are drawn from a common (but unknown) virtual signal $s \in \mathbb{R}^P$. The hierarchical representation of this effective data model is presented in Fig. 2. We assume that $s$ carries the deterministic information on the target parameters $\theta$, which is common to the data and to the
<table><caption>**Table 1.** The statistical variables appearing in Section 2.1 and their interpretation in the context of galaxy survey data analysis.</caption><thead><tr><th>Symbol</th><th>Meaning</th><th>Interpretation</th></tr></thead><tbody><tr><td>θ ∈ ℝ<sup>S</sup></td><td>Target parameters</td><td>Parametrization of the primordial matter power spectrum</td></tr><tr><td>ψ ∈ ℝ<sup>T</sup></td><td>Nuisance parameters</td><td>Random numbers involved in the initial phase realization, instrumental noise, etc.</td></tr><tr><td>d ∈ ℝ<sup>D</sup></td><td>Raw data</td><td>Galaxy number counts in a three-dimensional map of the survey volume</td></tr><tr><td>Φ ∈ ℝ<sup>P</sup></td><td>Summary statistics of the data (observed or simulated)</td><td>Summaries of the galaxy number count field, such as its estimated power spectrum</td></tr><tr><td>Φ<sub>O</sub> ∈ ℝ<sup>P</sup></td><td>Summary statistics of the observations</td><td>Summaries of the observed galaxy number count field</td></tr><tr><td>Φ<sub>θ</sub> ∈ ℝ<sup>P</sup></td><td>Summary statistics of simulated data</td><td>Summaries of a galaxy number count field, simulated with primordial matter power spectrum given by θ</td></tr><tr><td>s ∈ ℝ<sup>P</sup></td><td>Virtual signal</td><td>True summaries of the galaxy number count field, if they were not degraded by nuisances</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Hierarchical representation of the black-box model used in this work. The rounded green boxes represent probability distributions and the purple squares represent deterministic functions. The variables are $θ$ (the target parameters), $ψ$ (the nuisance parameters), $d$ (the full data), and $Φ$ (the summary statistics). The orange dashed rectangle represents the data-generating process, it gives the true (unknown) likelihood $L(θ)$ when $Φ = Φ_O$ (the summary statistics of the observations).</figcaption></figure>

mock observations; therefore, $Φ_O$ and $Φ_θ^{(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 $Φ$, which is degraded by nuisances $ψ$. Under this assumption, we introduce the probability distribution $\mathcal{P}(Φ|s)$, from which data realizations are drawn independently once the virtual signal $s$ is given. Since we do not want to infer $s$ explicitly, marginalization yields:

$$ \begin{aligned} \mathcal{P}(\Phi, \{\Phi_{\theta}^{(i)}\}|\theta) &= \int \mathcal{P}(\Phi, \{\Phi_{\theta}^{(i)}\}, s|\theta) ds \quad (8) \\ &= \int \mathcal{P}(\Phi|s) \mathcal{P}(\{\Phi_{\theta}^{(i)}\}|s) \mathcal{P}(s|\theta) ds. \end{aligned} $$

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

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Hierarchical representation of the approximate model used for inference of the black-box model of Fig. 1: a virtual signal $s$ has been introduced as a latent variable. The orange dashed rectangle represents the assumption made about the data-generating process, it gives the effective likelihood $\hat{L}^N(θ)$ when conditioning on a set of simulations $\{\Phi_θ^{(i)}\}$ and using $Φ = Φ_O$.</figcaption></figure>

that one can write

$$ \mathcal{P}(\Phi, \{\Phi_{\theta}^{(i)}\}|\theta) \propto \int \mathcal{P}(\Phi|s) \left[ \prod_{n=1}^{N} \mathcal{P}(\Phi_{\theta}^{(i)}|s) \right] ds. \quad (9) $$

In order to marginalize over $s$, we need to postulate a parametric form for the pdf $\mathcal{P}(Φ|s)$. In this work, we assume a Gaussian distribution centred on the virtual signal $s$ with a covariance matrix $\Sigma_θ$ quantifying the stochastic uncertainties inherent to the observations,<fnref n="2" />

$$ -2 \log \mathcal{P}(\Phi|s) = \log |2\pi \Sigma_{\theta} + (\Phi - s)^{\top} \Sigma_{\theta}^{-1} (\Phi - s)| \quad (10) $$

The virtual signal can therefore be interpreted as the expectation $E$ of $Φ_θ$ once $θ$ is specified: $s = E[Φ_θ]$. We also have $\Sigma_θ = E[(Φ_θ - s)(Φ_θ - s)^{\top}]$.

Under these assumptions, marginalization over $s$ gives the effective likelihood as $\hat{L}^N(θ) = \exp[\hat{\ell}^N(θ)]$ (the details of the computation are given in Appendix A), with

$$ -2\hat{\ell}^N(\theta) = \log |2\pi \hat{\Sigma}'_{\theta} + (\Phi_O - \hat{\Phi}_{\theta})^{\top} \hat{\Sigma}'_{\theta}^{-1} (\Phi_O - \hat{\Phi}_{\theta})|, \quad (11) $$

<footnote n="2">²The investigation of different choices for $\mathcal{P}(Φ|s)$ is left to future investigations.</footnote>
where

$$ \hat{\Phi}_{\theta} \equiv E^{N}[\Phi_{\theta}] = \frac{1}{N} \sum_{i=1}^{N} \Phi_{\theta}^{(i)} \quad (12) $$

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

$$ \hat{\Sigma}_{\theta} \equiv \frac{N+1}{N} \hat{\Sigma}_{\theta}, \quad \hat{\Sigma}_{\theta}^{-(1)} \equiv \left( \frac{N+1}{N} \right)^{-1} \hat{\Sigma}_{\theta}^{-1}, \quad (13) $$

where given a sufficiently large number of mock observations, a computable approximation of $\Sigma_{\theta}$ is estimated in the following way:

$$ \begin{aligned} \Sigma_{\theta} \approx \hat{\Sigma}_{\theta} &\equiv E^{N} \left[ (\Phi_{\theta} - \hat{\Phi}_{\theta})(\Phi_{\theta} - \hat{\Phi}_{\theta})^{\top} \right] \\ &= \frac{1}{N-1} \sum_{i=1}^{N} (\Phi_{\theta}^{(i)} - \hat{\Phi}_{\theta})(\Phi_{\theta}^{(i)} - \hat{\Phi}_{\theta})^{\top}. \end{aligned} \quad (14) $$

The effective likelihood also requires an estimator $\hat{\Sigma}_{\theta}^{-1}$ of the inverse covariance matrix $\Sigma_{\theta}^{-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:

$$ \Sigma_{\theta}^{-1} \approx \hat{\Sigma}_{\theta}^{-1} \equiv \alpha (\hat{\Sigma}_{\theta})^{-1}, \quad \alpha \equiv \frac{N-P-2}{N-1}, \quad (15) $$

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

The limiting approximation of $\hat{L}^N(\theta)$ when $N \to \infty$ is $\tilde{L}(\theta) = \mathcal{R}(\Phi|s)|_{\Phi=\Phi_O} = \exp[\tilde{\ell}(\theta)]$, with (as intended)
$$ -2\tilde{\ell}(\theta) \equiv \log|2\pi\Sigma_{\theta} + (\Phi_O - s)^{\top}\Sigma_{\theta}^{-1}(\Phi_O - s). \quad (16) $$

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 $\Phi_O$. Secondly, the mean of the effective likelihood is given by the empirical average of simulated summaries, $\hat{\Phi}_{\theta} = E^N[\Phi_{\theta}]$. This is an unbiased estimator of the virtual signal $s = E[\Phi_{\theta}]$ for a sufficiently large number of simulations. Importantly, the replacement of the expectation E by an empirical average $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} \hat{\Sigma}_{\theta}$, a multiple of the estimated covariance $\hat{\Sigma}_{\theta}$. 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 with the virtual signal increases the observed scatter. For example, for $N=1$, the covariance to be used in equation (11) is $2\hat{\Sigma}_{\theta}$, which reflects the fact that observed and simulated data $\Phi_O$ and $\hat{\Phi}_{\theta}$ can be drawn from opposite ends of the scatter around $s$ of $\mathcal{R}(\Phi|s)$. However, when $N \to \infty$, $\frac{N+1}{N} \hat{\Sigma}_{\theta} \to \Sigma_{\theta}$, the intrinsic covariance of $\mathcal{R}(\Phi|s)$. This result is reasonable, since for large sample sizes the ensemble mean $\hat{\Phi}_{\theta}$ converges to $s$, the mean of $\mathcal{R}(\Phi|s)$; in the same limit, we expect the covariance of the effective likelihood to approach $\Sigma_{\theta}$, the covariance of $\mathcal{R}(\Phi|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 $\mathcal{P}(\Phi|\theta)$ 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(\theta)$ has little sensitivity to our treatment and assumption for $\mathcal{P}(s|\theta)$. Furthermore, since $\frac{N+1}{N} \to 1$ when $N \to \infty$, our result with the virtual signal and the synthetic likelihood are equivalent, provided that the number of simulations is large enough.

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

## 2.2 Linearization of black-box models

The numerical cost of the computation of the effective likelihood described in Section 2.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 $\theta$ 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 I 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\theta$ around an expansion point $\theta_0$. Consequently, the target parameters are given by $\theta = \theta_0 + \Delta\theta$. Assuming that $\hat{\Phi}_{\theta}$ is differentiable with respect to $\theta$, we perform a first-order Taylor expansion in $\Delta\theta$ around $\theta_0$:

$$ \hat{\Phi}_{\theta} \approx f_0 + \nabla f_0 \cdot (\theta - \theta_0) \equiv f(\theta), \quad (17) $$

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

$$ (\nabla f_0)_S^\top \approx \frac{f(\theta_S) - f_0}{h}, \quad \text{with} \quad \theta_S \equiv \theta_0 + h (\delta_K^{SS})_{S' \in [1, S]}, \quad (18) $$

where $\delta_K$ is the Kronecker delta. To avoid obtaining a noisy gradient, nuisance parameters $\psi$ are kept at the same values in the simulations used to compute the $f(\theta_S)$ and $f_0$. Further, we neglect the dependence of the covariance matrix on $\theta$ and we estimate it at the expansion point according to equation (14), i.e. $\hat{\Sigma}'_{\theta} \approx \hat{\Sigma}'_{\theta_0} \equiv C_0$. Importantly, fully characterizing $f$ under these assumptions only requires evaluations of the simulator, thus ensuring that the data model remains a black box.
Using the linearized data model $f$ as a proxy for the data model described in Section 2.1.3 simplifies the expression of the effective likelihood (equation 11) to

$$-2\hat{\epsilon}^N(\theta) \approx \log|2\pi\mathbf{C}_0| + [\mathbf{\Phi}_O - f(\theta)]^\mathrm{T} \mathbf{C}_0^{-1} [\mathbf{\Phi}_O - f(\theta)]. \quad (19)$$

Evaluating $f_0$ and $C_0$ requires $N_0$ model evaluations at the expansion point $\theta_0$. The computation of the gradient further requires $N_s \times S$ model evaluations. Therefore, the linearized data model is fully characterized 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 $\theta_0$, a minimum would be $N_0 \ge P + 3$ but it can be worth investing more simulations, in order for the estimated covariance matrix $C_0$ and its inverse to be precise. Since all nuisance parameters are kept fixed in the computation of the $f(\theta_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 f_0$.

It is important to note that all required model evaluations are done once and for all. Once the linearized 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.

## 2.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 parametrize 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) = P(k)/P_0(k)$ as target function. Formally, with $P_0 \in \mathbb{R}^S$ the vector of components $(P_0)_s \equiv P_0(k_s)$, the inference variable $\theta$ is defined as the $S$-dimensional vector of components $P(k_s)/(P_0)_s$.

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

(i) the power spectrum is Gaussian-distributed,

(ii) it is strongly constrained to live close to $P_0$,

(iii) it is a smooth function of wavenumber,

(iv) and the power spectrum $P_0$ is subject to cosmic variance.

It follows from assumptions 1 and 2 that $\mathcal{P}(\theta)$ shall be a Gaussian distribution with mean $\theta_0 \equiv \mathbf{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 $\mathbf{S}$ of $\mathcal{P}(\theta)$. Let us define the matrix $\mathbf{K}$ as a radial basis function, i.e. by its coefficients

$$\mathbf{K}_{ss'} \equiv \exp \left[ -\frac{1}{2} \left( \frac{k_s - k_{s'}}{k_{\text{corr}}} \right)^2 \right]. \quad (20)$$

The hyperparameter $k_{\text{corr}}$ determines the length-scale on which power spectrum amplitudes of different wavenumber correlate with each other. A large value of $k_{\text{corr}}$ corresponds to a strong correlation between $\theta_s$ and $\theta_{s'}$, even if their corresponding wavenumbers are far from each other. On the other hand, for $k_{\text{corr}} \ll k_s, k_{s'}$, $\mathbf{K}$ becomes the identity matrix and $\theta_s$ and $\theta_{s'}$ 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_{\text{corr}}$. In a realistic scenario, the prior covariance is not scale-independent as is $\mathbf{K}$ ($\mathbf{K}_{ss} = 1$ for all $s$). We rather want the standard deviations $(\mathbf{S})_{ss}^{1/2}$ to account for the cosmic variance affecting the power spectrum $P_0$ (and hence the mean $\theta_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 $\theta$, the coefficients $\mathbf{K}$ shall be multiplied by the coefficients of $\mathbf{u}\mathbf{u}^\mathrm{T}$, where

$$\mathbf{u}_s \equiv 1 + \sigma_s = 1 + \frac{\alpha_{cv}}{k_s^{3/2}} \quad (21)$$

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

$$\mathbf{S} \equiv \theta_{\text{norm}}^2 \mathbf{u}\mathbf{u}^\mathrm{T} \circ \mathbf{K}, \quad (22)$$

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

The resulting prior on $\theta$ is characterized by a set of three hyperparameters $\{k_{\text{corr}}, \alpha_{cv}, \theta_{\text{norm}}\}$, and given as

$$-2 \log \mathcal{P}(\theta) \equiv \log|2\pi\mathbf{S}| + (\theta - \theta_0)^\mathrm{T}\mathbf{S}^{-1}(\theta - \theta_0). \quad (23)$$

As demonstrated in Section 4, 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.

## 2.4 The power spectrum effective posterior distribution

Using the effective likelihood with the linearized 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,

$$-2 \log \mathcal{P}(\theta|\mathbf{\Phi})|_{\mathbf{\Phi}=\mathbf{\Phi}_O} = \log|2\pi\mathbf{\Gamma}| + (\theta - \mathbf{\gamma})^\mathrm{T}\mathbf{\Gamma}^{-1}(\theta - \mathbf{\gamma}), \quad (24)$$

with mean

$$\mathbf{\gamma} \equiv \theta_0 + \mathbf{\Gamma} (\nabla f_0)^\mathrm{T} \mathbf{C}_0^{-1} (\mathbf{\Phi}_O - f_0), \quad (25)$$

and covariance matrix

$$\mathbf{\Gamma} \equiv [(\nabla f_0)^\mathrm{T} \mathbf{C}_0^{-1} \nabla f_0 + \mathbf{S}^{-1}]^{-1}. \quad (26)$$

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 summarizes the variables appearing in equations (25) and (26) and their interpretation in the context of galaxy survey data analysis.

## 2.5 Optimization of the prior hyperparameters

As discussed in Section 2.3, the chosen prior distribution involves three hyperparameters $\{k_{\text{corr}}, \alpha_{cv}, \theta_{\text{norm}}\}$. By definition, $\alpha_{cv}$ characterizes 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_{cv}^2$. $\alpha_{cv}$ can therefore simply be measured for
<table><caption>**Table 2.** The statistical variables appearing in SELFI and their interpretation in the context of galaxy survey data analysis.</caption><thead><tr><th>Symbol</th><th>Meaning</th><th>Interpretation</th></tr></thead><tbody><tr><td>θ<sub>0</sub> ∈ ℝ<sup>S</sup></td><td>Expansion point of the simulator in parameter space and prior mean</td><td>Fiducial primordial matter power spectrum</td></tr><tr><td>S ∈ ℝ<sup>S×S</sup></td><td>Prior covariance matrix</td><td>Prior covariance matrix of the primordial matter power spectrum</td></tr><tr><td>Φ<sub>O</sub> ∈ ℝ<sup>P</sup></td><td>Summary statistics of the observations</td><td>Summaries of the observed galaxy number count field</td></tr><tr><td>f<sub>0</sub> ∈ ℝ<sup>P</sup></td><td>Estimated average black box at the expansion point</td><td>Estimated average of summaries of simulated galaxy fields with fiducial primordial matter power spectrum</td></tr><tr><td>C<sub>0</sub> ∈ ℝ<sup>P×P</sup></td><td>Estimated covariance matrix of the black box at the expansion point</td><td>Estimated covariance of summaries of simulated galaxy fields with fiducial primordial matter power spectrum</td></tr><tr><td>∇f<sub>0</sub> ∈ ℝ<sup>P×S</sup></td><td>Estimated gradient of the average black box at the expansion point</td><td>Estimated gradient of summaries of simulated galaxy fields around the fiducial primordial matter power spectrum</td></tr><tr><td>γ ∈ ℝ<sup>S</sup></td><td>Posterior mean</td><td>Reconstructed primordial matter power spectrum, given the observed summaries Φ<sub>O</sub></td></tr><tr><td>Γ ∈ ℝ<sup>S×S</sup></td><td>Posterior covariance matrix</td><td>Uncertainties on the reconstruction of the primordial matter power spectrum, given the observed summaries Φ<sub>O</sub></td></tr></tbody></table>

a given simulator setup (i.e. box size and mesh). On the contrary, $k_{corr}$ and $\theta_{norm}$ are free hyperparameters and the power spectrum reconstruction ($\gamma$ and $\Gamma$) generally depends on their values. In this section, we propose a procedure to find optimal values for $k_{corr}$ and $\theta_{norm}$.

$k_{corr}$ indicates the a priori smoothness of reconstructed wiggle functions $\theta(k)$ and $\theta_{norm}$ how much they can deviate from the expansion point. Together, these two parameters characterize 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 optimization in machine learning, we propose to optimize $k_{corr}$ and $\theta_{norm}$ to reproduce the shape of a fiducial wiggle function $\theta_{fid}$, using as likelihood the effective posterior distribution derived in Section 2.4. More precisely, the likelihood for $k_{corr}$ and $\theta_{norm}$ is defined by

$$ -2 \log \mathcal{P}(k_{\text{corr}}, \theta_{\text{norm}} | \theta_{\text{fid}}) = \log |2\pi \Gamma| + (\theta_{\text{fid}} - \gamma)^{\top} \Gamma^{-1} (\theta_{\text{fid}} - \gamma), \quad (27) $$

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

Note that evaluating the likelihood given in equation (27) is cheap once the linearized 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_{corr}$ and $\theta_{norm}$ to be found using standard optimizers.

## 2.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 parametrization 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 Fig. 3): $\omega$ is a vector of cosmological parameters which generates the primordial power spectrum coefficients $\theta$. This generative process is usually deterministic, i.e.

$$ \mathcal{P}(\theta|\omega) = \delta_D(\theta - \mathcal{T}(\omega)), \quad (28) $$

where $\mathcal{T}$ is a deterministic function of cosmological parameters, typically a Boltzmann solver or a fitting function. Given this

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Hierarchical representation of the Bayesian model for inference of cosmological parameters $\omega$: a layer has been added above the primordial power spectrum $\theta$. The purple square represents the deterministic process generating $\theta$ given $\omega$.</figcaption></figure>

assumption, we have $\mathcal{P}(\Phi|\omega) = \mathcal{P}(\Phi|\mathcal{T}(\omega))$, and the inference of cosmological parameters gives

$$ \mathcal{P}(\omega|\Phi)|_{\Phi=\Phi_O} = \hat{L}_{\omega}^{N}(\omega) \frac{\mathcal{P}(\omega)}{Z_{\Phi}}, \quad (29) $$

where the likelihood for cosmological parameters is defined as

$$ \begin{aligned} \hat{L}_{\omega}^{N}(\omega) &\equiv \mathcal{P}(\Phi|\mathcal{T}(\omega))|_{\Phi=\Phi_O} \\ &= \hat{L}^{N}(\mathcal{T}(\omega)). \end{aligned} \quad (30) $$

Noting $\hat{L}_{\omega}^{N}(\omega) \equiv \exp[\hat{\ell}_{\omega}^{N}(\omega)]$ and using the linearized data model, we have (see equation 19)

$$ -2\hat{\ell}_{\omega}^{N}(\omega) \approx \log|2\pi C_0| + [\Phi_O - f(\mathcal{T}(\omega))]^{\top} C_0^{-1} [\Phi_O - f(\mathcal{T}(\omega))]. \quad (31) $$

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
<table><caption>**Table 3.** The cosmological parameters used in this work.</caption><thead><tr><th><i>h</i></th><th>Ω<sub>b</sub></th><th>Ω<sub>m</sub></th><th><i>n</i><sub>s</sub></th><th>σ<sub>8</sub></th></tr></thead><tbody><tr><td>0.6674</td><td>0.0486</td><td>0.3089</td><td>0.9667</td><td>0.8159</td></tr></tbody></table>

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 centred on the data $\Phi_O$ with a fixed covariance matrix $C_O$. Furthermore, getting predictions for $\Phi$ amounts to evaluating $f(\mathcal{I}(\omega))$, which does not require additional black-box evaluations once $f$ is known. The Bayesian problem of inferring cosmological parameters $\omega$ (equation 29) can therefore easily be solved by standard MCMC techniques.

## 3 DATA MODEL

In this section, we describe the generation of mock observations used to test the performance of our method. We emphasize that the statistical method presented in Section 2 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 4 as a showcase, to highlight the performance of our method to handle a complex simulator.

### 3.1 Primordial power spectrum parametrization

Throughout this paper, we work with a cubic equidistant grid with comoving side length of $1 \text{ Gpc h}^{-1}$ and $256^3$ voxels, spanning scales between $k_{s,min} = 6.28 \times 10^{-3} \text{ h Mpc}^{-1}$ and $k_{s,max} = 1.4 \text{ h Mpc}^{-1}$. 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_{max}$. Any vector $\theta$ 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.

### 3.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 XIII 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.<fnref n="3" />

SIMBELMYNË first generates a realization 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

$$ \delta_k^i \equiv \sqrt{P(k)} w_k. \qquad (32) $$

Transformation back to configuration space yields an initial density contrast field $\delta^i$. One realization of such an initial density field is shown in the left-hand panel of Fig. 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 (see e.g. Moutarde et al. 1991; Bouchet et al. 1995; Bouchet 1996), then with an efficient implementation of COLA (Tassev et al. 2013) from $z=19$ to $z=0$. A particle-mesh grid of $1024^3$ voxels and 20 time-steps 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 per cent precision up to $k = 0.5 \text{ h Mpc}^{-1}$ and 5 per cent up to $k = 1 \text{ h Mpc}^{-1}$. We place the observer at the centre of the box. The maximal distance to the observer is $866 \text{ Mpc h}^{-1}$, 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

$$ 1 + z_{\text{obs}} = (1 + z_{\text{cosmo}})(1 + z_{\text{pec}}), \quad \text{with } z_{\text{pec}} \equiv -\frac{v_r}{c}, \qquad (33) $$

where $z_{\text{cosmo}}$ is the true cosmological redshift, $z_{\text{obs}}$ is the 'observed' redshift and $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 $\delta^i$. One realization of the final redshift-space density is shown in the middle panel of Fig. 4.

The galaxy density $\rho^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$,

$$ \rho_x^g = \bar{N} (1 + b \delta_x^f). \qquad (34) $$

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} (\text{h Mpc}^{-1})^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 marginalized 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) $\mathbf{W}$, consisting of the product of the radial selection function $R(r)$ and the angular survey mask and completeness function $C(\hat{n})$ for any line-of-sight $\hat{n}$, i.e. $\mathbf{W}(\hat{n}, r) \equiv R(r) C(\hat{n})$. This operator accounts for the fact that we are only looking at certain

<footnote n="3">³In principle, our approach would require treating cosmological parameters as nuisance parameters in the inference of the primordial power spectrum, and to marginalize over them. In this work, for simplicity, we keep them fixed to the values used to predict the expansion point.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Slices through one realization of the initial density field (left) scaled to redshift zero with the linear growth factor, the corresponding evolved density field in redshift space (middle), and the field of observed galaxy number counts (right), according to the data model used in this work. The fields are defined on a grid of 256³ cells covering a total volume of (1 Gpc h⁻¹)³.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Radial selection function R(r) (solid blue line) used in our synthetic survey as a function of the comoving distance to the observer r. The model used is a Schechter luminosity function with parameters α = -1.05 and M* = -20.44, a limiting apparent magnitude of m = 18.5 and absolute magnitude cuts -25 ≤ M ≤ -21. The luminosity distance D_L (dashed red line, right vertical axis) is also shown. The survey is complete up to D_L = 794 Mpc h⁻¹, corresponding to r = 684 Mpc h⁻¹.</figcaption></figure>

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³ voxels covering a volume of (1 Gpc h⁻¹)³] as follows. For the angular completeness, we mask the galactic plane by setting C(ˆn) to 0 for galactic latitudes -10° ≤ b ≤ 10° 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: α = -1.05 and M* = -20.44 (Blanton et al. 2003), a limiting apparent magnitude of m = 18.5 and absolute magnitude cuts -25 ≤ M ≤ -21. This choice makes our synthetic survey complete up to a luminosity distance of D_L = 794 Mpc h⁻¹, corresponding to r = 684 Mpc h⁻¹ (see Fig. 5).

We emulate the survey via galaxy number counts N_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_g is drawn from a Gaussian distribution with mean μ_g and standard deviation σ_x,

$$N_g \sim G(\mu_g | \sigma_x). \qquad (35)$$

The mean μ_g = W_x ρ_x = W_x N̄(1 + b δ_x^f) characterizes the expected number of galaxies. It accounts for all physical effects via δ_x^f and for the survey response operator W_x. The standard deviation is defined as σ_x = σ√W_x N̄, where we set the overall noise level σ = 10⁻¹. A slice through one realization of the galaxy number count field N_g is shown in the right-hand panel of Fig. 4.

## 3.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 N_g, but use a compression 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 N_g, summing the contributions over all the Fourier cells within wavenumber shell k_r, and normalizing:

$$P^r(k_r) \equiv C \times \sum_{|k| \in k_r} \frac{|N_k^g|^2}{N_{k_r} - 2}. \qquad (36)$$

The overall constant factor C = 1000³/256² (Mpc h⁻¹)³ 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-Γ distributed with shape parameter N_{k_r}/2 and scale parameter C × |N_k^g|^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
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Diagnostics of the black box at the expansion point $\theta_0$. The left-hand panel shows individual mock observations $\Phi_{\theta_0}$ (grey lines), the observed data $\Phi_O$ (solid black line), and the average black-box $f_0$ (dashed purple line). The shaded region corresponds to two standard deviations. The right-hand panel shows the covariance matrix of the summaries at the expansion point, $C_0$. The dashed red lines correspond to the positions of the bins at which summaries are measured in data space.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Gradient of the black box, measured via finite differencing. The left-hand panel shows individual columns $(\nabla f_0)_s^T$ for $s = 20, 50$, and $70$, corresponding to support wavenumbers $k_s = 0.0364, 0.1484$, and $0.3780 \, h \, \text{Mpc}^{-1}$, respectively. The right-hand panel shows the gradient matrix $\nabla f_0$. The dashed grey line corresponds to the identity function; the dashed red lines correspond to the positions of the bins at which summaries are measured in data space; and the dotted green lines correspond to the support wavenumbers in parameter space.</figcaption></figure>

in the range $[0.02, 0.5] \, h \, \text{Mpc}^{-1}$, ensuring that each bin contains at least 100 modes. These bins are logarithmically spaced for $k_r \ge 0.04 \, h \, \text{Mpc}^{-1}$. For convenience, we normalize the output of the black box using the expansion point, so that $\Phi$ is the $P$-dimensional vector of components $A \times P^f(k_r)/P_0(k_r)$ with $A = 50$. The estimation of $P^f(k)$ is completely deterministic once the galaxy number counts are given. Therefore, we now have a complete model to generate artificial realizations of $\Phi$ for a given primordial power spectrum $P(k)$ and specific realizations 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_{\text{max}} \lesssim 0.3 \, h \, \text{Mpc}^{-1}$.

(e.g. Ross et al. 2015). Pushing the analysis to $k_{\text{max}} = 0.5 \, h \, \text{Mpc}^{-1}$ represents an increase by a factor of $\sim 5$ in the number of modes used (scaling as $k_{\text{max}}^3$), which is expected to yield substantial improvements in the inference results.

## 3.4 Idealized data model

For testing purposes, we also define an idealized 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 $\delta_i$ (see Section 3.2), scaling it to redshift zero using the linear growth factor, and measuring its normalized power spectrum $\Phi$ (see Section 3.3). The use of the two
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Posterior surface for the prior hyperparameters $k_{corr}$ and $\theta_{norm}$. The $1\sigma$, $2\sigma$, and $3\sigma$ credible contours are shown as solid blue lines. The dashed red lines mark the maximum a posteriori values.</figcaption></figure>

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.

## 4 RESULTS

This section describes the results obtained by applying the statistical method proposed in Section 2 in conjunction with the data-generating process described in Section 3 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 $\omega_{gt}$ are drawn from the (marginalized, Gaussian) Planck priors:

$$ \begin{pmatrix} h \\ \Omega_b \\ \Omega_m \\ n_s \\ \sigma_8 \end{pmatrix} \rightleftharpoons \mathcal{G} \left[ \begin{pmatrix} 0.6774 \\ 0.04860 \\ 0.3089 \\ 0.9667 \\ 0.8159 \end{pmatrix}, \text{diag} \left[ \begin{pmatrix} 0.0046^2 \\ 0.00030^2 \\ 0.0062^2 \\ 0.0040^2 \\ 0.0086^2 \end{pmatrix} \right] \right]. \quad (37) $$

The ‘wiggly’ ground truth power spectrum $P_{gt}(k)$ is generated with the Eisenstein & Hu (1998, EH) fitting function, using these cosmological parameters. It is used to simulate observed data $\Phi_O$, with unknown nuisance parameters (phase realization and instrumental noise). For later use, the fiducial ‘wiggly’ power spectrum $P_{fid}(k)$ is also generated with the EH prescription, using Planck cosmology. The target parameters $(\theta)_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 3.1. We note $\theta_{gt}$ and $\theta_{fid}$ the vectors of component $P_{gt}(k_s)/P_0(k_s)$ and $P_{fid}(k_s)/P_0(k_s)$, respectively.

### 4.1 Diagnostics of the black box

We created an ensemble of $N_0 = 150$ mock realizations at the expansion point $\theta_0 = 1_{RS}$ using different nuisance parameters. These are used to compute the average black box at the expansion point, $f_0 \equiv \hat{\Phi}_{\theta_0}$, and the covariance matrix $C_0 \equiv \hat{\Sigma}_{\theta_0}$, using their definitions (equations 12, 13, and 14). The results are shown in

Fig. 6. There, the left-hand panel shows individual realizations $\Phi_{\theta_0}$ and the observed data vector $\Phi_O$. The average black-box $f_0$ is also plotted, with a credible region corresponding to two standard deviations (i.e. $2\sqrt{\text{diag}(C_0)}$). The full estimated covariance matrix $C_0$ is shown in the right-hand panel. As expected, the measured variance is larger on large scales due to cosmic variance, with some anticorrelations 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 realizations at each of the expansion points $\theta_s$, we measured the gradients of the black-box $(\nabla f_0)_s^\mathrm{T}$ along all directions of parameter space (see equation 18). The nuisance parameters (phase realization and noise) are kept at fixed values (the ones corresponding to the first $N_s$ realizations generated at the expansion point) for this calculation. The results are shown in Fig. 7, where the left-hand panel shows $(\nabla f_0)_s^\mathrm{T}$ for individual values of $s$ and the right-hand panel shows the full rectangular matrix $\nabla f_0$. Some interesting phenomena can be observed. At large scales, (see e.g. for $k_{20} = 0.0364 \, h \, \text{Mpc}^{-1}$) exciting one initial mode only 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 \, \text{Mpc}^{-1}$). 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 2.2, the linearized black-box $f(\theta)$ is fully characterized by $f_0$, $C_0$ and $\nabla 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 $C_0$ and $\nabla f_0$, although using fewer would have been possible.

### 4.2 The prior and its optimization

As discussed in Section 2.3, we choose a Gaussian prior centred on the expansion point $\theta_0$ with a covariance matrix $S$ given by equation (22) and characterized by three hyperparameters $\{k_{corr}, \alpha_{cv}, \theta_{norm}\}$.

Following the method presented in Section 2.5, we found optimal values for the prior hyperparameters. The strength of cosmic variance within our simulation volume shall satisfy $\alpha_{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 3.1), we measured up to the Nyquist frequency $\alpha_{cv} = 8.848 \times 10^{-4}$, value that we adopt.

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

### 4.3 The effective posterior

Using the effective likelihood, characterized by the observed data $\Phi_O$ and the linearized black box described in Section 4.1, as well
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Prior covariance matrix $\mathbf{S}$ for the wiggle function $\theta$ (left-hand panel). The corresponding prior covariance matrix for power spectrum amplitudes, $\mathrm{diag}(\mathbf{P}_0) \cdot \mathbf{S} \cdot \mathrm{diag}(\mathbf{P}_0)$, is shown in the right-hand panel. The dotted green lines correspond to the support wavenumbers in parameter space.</figcaption></figure>

as the optimized prior discussed in Section 4.2, we obtained the effective posterior on $\theta$. It is a Gaussian, with mean and covariance matrix given by the ‘filter equations’ (25) and (26).

Fig. 10 shows the inferred primordial wiggle function $\gamma$ in comparison with the expansion point $\theta_0$, the ground truth $\theta_{gt}$, and the fiducial function used to optimize prior hyperparameters, $\theta_{hd}$. 2$\sigma$ credible regions are shown for the prior and the posterior (i.e. $2\sqrt{\mathrm{diag}(\mathbf{S})}$ and $2\sqrt{\mathrm{diag}(\mathbf{\Gamma})}$, respectively). The top panel corresponds to the result obtained using the idealized data model (a Gaussian random field, see Section 3.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 parametrization 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 idealized case, and in the realistic case, up to scales of $k \approx 0.3 \, h \, \mathrm{Mpc}^{-1}$. In particular, five 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^{-1})^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 \, \mathrm{Mpc}^{-1}$ 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 idealized as well as in the realistic case. On the other hand, noise (understood as the combined effect of the specific phase realization of the data, non-linear gravity, redshift-space distortions and instrumental noise) acts on smaller scales ($0.2 \, h \, \mathrm{Mpc}^{-1} \lesssim k$

$\lesssim 0.5 \, h \, \mathrm{Mpc}^{-1}$). Some of the primordial information is effectively destroyed at these scales – or at least is not captured by the statistical summaries $\mathbf{\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 idealized case. As expected, we recover the prior at $k \gtrsim 0.5 \, h \, \mathrm{Mpc}^{-1}$ (in fact a little below, due to mode coupling), since the data $\mathbf{\Phi}$ do not contain measurements at these scales.

It is important to note that the prior $\mathcal{P}(\theta)$, used in this section to regularize 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. Fig. 11 shows the covariance matrix $\mathbf{\Gamma}$ of the Gaussian effective posterior; it can be compared with the prior covariance matrix $\mathbf{S}$ shown in Fig. 9.

## 4.4 Cosmological parameters

The last step in the analysis is to infer parameters $\omega$ of specific cosmological models, given the observed data $\mathbf{\Phi}_O$. In this section, we assume a flat $\Lambda$CDM model, characterized by five parameters $\{h, \Omega_b, \Omega_m, n_s, \sigma_8\}$. For simplicity, our prescription to generate primordial power spectra given cosmological parameters is the EH fitting function. The generative process $\mathcal{T}$ (from $\omega$ to $\theta$) is therefore the EH fitting function divided by $P_0$, sampled at the support wavenumbers $k_s$. It could easily be generalized to extensions of the flat $\Lambda$CDM model and to include a Boltzmann solver. The linearized data model $\mathbf{f}$ (from $\theta$ to $\mathbf{\Phi}$) has already been characterized in Section 4.1.

The effective likelihood $\hat{\mathcal{L}}_{\omega}^{N}(\omega)$ for cosmological parameters is given by equation (31). Consistently with the expansion point used to linearize the black box, we complement $\hat{\mathcal{L}}_{\omega}^{N}(\omega)$ with a Gaussian prior $\mathcal{P}(\omega)$ centred on Planck cosmological parameters,
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Inference of the wiggle function $\theta(k) \equiv P(k)/P_0(k)$ as a function of wavenumber $k$, using as data model a Gaussian random field (top panel) or a realistic mock survey (bottom panel). The prior mean $\theta_0$ and the effective posterior mean $\gamma$ are represented as solid yellow and green lines, respectively, with their $2\sigma$ credible intervals (for the prior, $2\sigma = 2\theta_{\text{norm}}(1 + \alpha_{\text{cv}}/k^{3/2})$). For comparison, the ground truth $\theta_{\text{gt}}$ and the fiducial ‘wiggly’ function $\theta_{\text{fid}}$ used to optimize the prior hyperparameters are plotted as solid blue and dashed orange line, respectively. The dashed red lines correspond to the positions of the bins at which summaries are measured in data space; and the dotted green lines correspond to the support wavenumbers in parameter space. In the realistic case, in spite of survey complications which limit the information captured, the signature of BAOs is well reconstructed up to $k \approx 0.3 \, \text{h Mpc}^{-1}$, with five inferred acoustic peaks, result which could be improved using more volume. In the absence of informative data, the power spectrum reconstruction is driven towards the prior mean, but this effect does not affect cosmological parameter inference (see Section 4.3 for details).</figcaption></figure>

but with broader variance: the diagonal covariance matrix given in equation (37) is multiplied by a factor of 3. We explored the effective posterior $\mathcal{P}(\omega|\Phi)|_{\Phi=\Phi_0}$ (equation 29) via MCMC (performed using the EMCEE code; Foreman-Mackey et al. 2013), ensuring sufficient convergence.

Results are shown in Fig. 12. Prior contours are shown in blue, and posterior contours are shown in red and purple for two different realizations of the data $\Phi_0$. The two data realizations have been generated using different ground truth cosmological parameters $\omega_{\text{gt}}$ (shown as dashed and dotted lines, respectively), as well as different nuisance parameters (phase realization and noise). The plot demonstrates that cosmological parameter inference is unbiased and robust to nuisances imprinted in the data.

# 5 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
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Posterior covariance matrix **Γ** for the wiggle function θ (left-hand panel). The corresponding posterior covariance matrix for power spectrum amplitudes, diag(*P*₀) · **Γ** · diag(*P*₀), is shown in the right-hand panel. The dotted green lines correspond to the support wavenumbers in parameter space. The correlation structure between different inferred parameters exhibits a complex behaviour (see also Fig. 9 for comparison).</figcaption></figure>

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 linearizing the black
box around an expansion point. We devised a method to optimize 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 dif-
ferences are as follows.

(i) 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.

(ii) Secondly, 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 re-
alizations to characterize the effective likelihood with thelinearized
black box, all of which are used to obtain the inference result. For
*S* target parameters and *N*₀ nuisance parameters realizations at the
expansion point, *N*ₛ nuisance parameters realizations along each
direction in parameter space, one has to perform *N*₀ + *N*ₛ × *S* data
model evaluations. Increasing *N*₀ or *N*ₛ to get a better estimate of
the required covariance matrix or gradient only increases the overall
computational cost linearly.

(iii) Thirdly, 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.

(iv) Finally, the linearized 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* ≳ 100 pa-
rameters, as is necessary for a model-independent parametrization
of theory in cosmology.

In this work, we demonstrated that a ‘non-wiggly’ expansion
point θ₀ is sufficient to recover the target wiggle function θ 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 γ would
be used as the new expansion point to train a new linearized 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 differ-
encing to obtain the gradient of the averaged black-box ∇f₀, but our
equations could be used with other techniques, such as automatic
differentiation. The data covariance matrix at the expansion point
C₀ 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,
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Cosmological parameter inference using a linearized black-box model of galaxy surveys. The prior is shown in blue, and the effective posteriors for two different data realizations are shown in red and purple. The two different data realizations have different input cosmological parameters (shown as dashed and dotted lines), different phases and noise realizations. For all distributions, the $1\sigma$ and $2\sigma$ contours are shown.</figcaption></figure>

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_{max} = 0.5 \ h\ Mpc^{-1}$. This represents an increase by a factor of ~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 ~$\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.<fnref n="4" /> 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

<footnote n="4">⁴Currently, 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.</footnote>
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.

## STATEMENT OF CONTRIBUTION

Study concept and design (JJ, FL, WE); design of prior optimization 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).

## REFERENCES

Akeret J., Refregier A., Amara A., Seehars S., Hasner C., 2015, J. Cosmol. Astropart. Phys., 8, 043

Albrecht A. et al., 2006, preprint (arXiv:astro-ph/0609591)

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

Angulo R. E., Pontzen A., 2016, MNRAS, 462, L1

Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986, ApJ, 304, 15

Blanton M. R. et al., 2003, ApJ, 592, 819

Bocquet S., Carter F. W., 2016, J. Open Source Softw., 1, 46

Bos E. G. P., Kitaura F.-S., van de Weygaert R., 2019, MNRAS, 488, 2573

Bouchet F. R., 1996, in Bonometto S., Primack J. R., Provenzale A., eds, Dark Matter in the Universe. IOS Press, Oxford, GB, p. 565

Bouchet F. R., Colombi S., Hivon E., Juszkiewicz R., 1995, A&A, 296, 575

Burden A., Percival W. J., Howlett C., 2015, MNRAS, 453, 456

Byrd R. H., Lu P., Nocedal J., Zhu C., 1995, SIAM J. Sci. Comput., 16, 1190

Doumler T., Hoffman Y., Courtois H., Gottlöber S., 2013, MNRAS, 430, 888

Eisenstein D. J., 2005, New Astron Rev., 49, 360

Eisenstein D. J., Hu W., 1998, ApJ, 496, 605

Eisenstein D. J., Seo H.-J., White M., 2007a, ApJ, 664, 660

Eisenstein D. J., Seo H.-J., Sirko E., Spergel D. N., 2007b, ApJ, 664, 675

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

Granett B. R. et al., 2015, A&A, 583, A61

Hall A., Taylor A., 2019, MNRAS, 483, 189

Hartlap J., Simon P., Schneider P., 2007, A&A, 464, 399

Hockney R. W., Eastwood J. W., 1981, Computer Simulation Using Particles. McGraw-Hill, New York

Ishida E. E. O. et al., 2015, Astron. Comput., 13, 1

Jasche J., Lavaux G., 2015, MNRAS, 447, 1204

Jasche J., Lavaux G., 2017, A&A, 606, A37

Jasche J., Lavaux G., 2019, A&A, 625, A64

Jasche J., Wandelt B. D., 2013a, MNRAS, 432, 894

Jasche J., Wandelt B. D., 2013b, ApJ, 779, 15

Jasche J., Kitaura F. S., Wandelt B. D., Enßlin T. A., 2010, MNRAS, 406, 60

Jasche J., Leclercq F., Wandelt B. D., 2015, J. Cosmol. Astropart. Phys., 1, 036

Jeffrey N., Abdalla F. B., 2018, preprint (arXiv:1809.08246)

Jennings E., Madigan M., 2017, Astron. Comput., 19, 16

Laureijs R. et al., 2011, preprint (arXiv:1110.3193)

Lavaux G., Jasche J., 2016, MNRAS, 455, 3169

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

Leclercq F., Jasche J., Wandelt B., 2015, J. Cosmol. Astropart. Phys., 6, 15

LSST Science Collaboration, 2012, preprint (arXiv:1211.0310)

Majerotto E. et al., 2012, MNRAS, 424, 1392

Meiksin A., White M., Peacock J. A., 1999, MNRAS, 304, 851

Moutarde F., Alimi J.-M., Bouchet F. R., Pellat R., Ramani A., 1991, ApJ, 382, 377

Padmanabhan N., Xu X., Eisenstein D. J., Scalzo R., Cuesta A. J., Mehta K. T., Kazin E., 2012, MNRAS, 427, 2132

Peacock J. A., 1999, Cosmological Physics. Cambridge Univ. Press, Cambridge

Peacock J. A., Heavens A. F., 1985, MNRAS, 217, 805

Peebles P. J. E., 1980, The Large-Scale Structure of the Universe. Princeton Univ. Press, Princeton, NJ

Percival W. J., Cole S., Eisenstein D. J., Nichol R. C., Peacock J. A., Pope A. C., Szalay A. S., 2007, MNRAS, 381, 1053

Planck Collaboration XIII, 2016, A&A, 594, A13

Planck Collaboration I, 2018a, preprint (arXiv:1807.06205)

Planck Collaboration I, 2018b, preprint (arXiv:1807.06209)

Price L. F., Drovandi C. C., Lee A., Nott D. J., 2018, J. Comput. Graph. Stat., 27, 1

Rasmussen C. E., Williams C. K. I., 2006, Gaussian Processes for Machine Learning. Adaptative Computation and Machine Learning SSeries. University Press Group Limited

Ross A. J. et al., 2012, MNRAS, 424, 564

Ross A. J., Samushia L., Howlett C., Percival W. J., Burden A., Manera M., 2015, MNRAS, 449, 835

Schechter P., 1976, ApJ, 203, 297

Seljak U., Aslanyan G., Feng Y., Modi C., 2017, J. Cosmol. Astropart. Phys., 12, 009

Sellentin E., Heavens A. F., 2016, MNRAS, 456, L132

Silk J., 1968, ApJ, 151, 459

Smith R. E., Marian L., 2015, MNRAS, 454, 1266

Springel V., 2005, MNRAS, 364, 1105

Tassev S., Zaldarriaga M., Eisenstein D. J., 2013, J. Cosmol. Astropart. Phys., 6, 36

Verde L. et al., 2002, MNRAS, 335, 432

Villaescusa-Navarro F. et al., 2018, ApJ, 867, 137

Wandelt B. D., Larson D. L., Lakshminarayanan A., 2004, Phys. Rev. D, 70, 083511

Wang H., Mo H. J., Yang X., van den Bosch F. C., 2013, ApJ, 772, 63

Wang H., Mo H. J., Yang X., Jing Y. P., Lin W. P., 2014, ApJ, 794, 94

White M., 2015, MNRAS, 450, 3822

Wiener N., 1964, Extrapolation, Interpolation, and Smoothing of Stationary Time Series: With Engineering Applications. Principles of Electrical Engineering Series. M.I.T. Press

Wood S. N., 2010, Nature, 466, 1102
## 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 2.1.3. We have

$$ \mathcal{P}(\mathbf{\Phi}, \{\mathbf{\Phi}_{\theta}^{(i)}\}|\theta) \propto \int \mathcal{P}(\mathbf{\Phi}|s) \left[ \prod_{n=1}^{N} \mathcal{P}(\mathbf{\Phi}_{\theta}^{(i)}|s) \right] ds. \quad (A1) $$

Using the parametric form for $\mathcal{P}(\mathbf{\Phi}|s)$ (equation 10) yields $\mathcal{P}(\mathbf{\Phi}, \{\mathbf{\Phi}_{\theta}^{(i)}\}|\theta) \propto \int \exp[\hat{\ell}_1(\mathbf{\Phi}, \{\mathbf{\Phi}_{\theta}^{(i)}\}, s)] ds$, with

$$ -2\hat{\ell}_1(\mathbf{\Phi}, \{\mathbf{\Phi}_{\theta}^{(i)}\}, s) = (\mathbf{\Phi} - s)^{\mathrm{T}} \mathbf{\Sigma}_{\theta}^{-1} (\mathbf{\Phi} - s) + \sum_{i=1}^{N} (\mathbf{\Phi}_{\theta}^{(i)} - s)^{\mathrm{T}} \mathbf{\Sigma}_{\theta}^{-1} (\mathbf{\Phi}_{\theta}^{(i)} - s) + (N+1) \log |2\pi \mathbf{\Sigma}_{\theta}|. \quad (A2) $$

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

$$ \begin{aligned} -2\hat{\ell}_1(\mathbf{\Phi}, \{\mathbf{\Phi}_{\theta}^{(i)}\}, s) &= \mathbf{\Phi}^{\mathrm{T}} \mathbf{\Sigma}_{\theta}^{-1} \mathbf{\Phi} + \sum_{i=1}^{N} \mathbf{\Phi}_{\theta}^{(i)\mathrm{T}} \mathbf{\Sigma}_{\theta}^{-1} \mathbf{\Phi}_{\theta}^{(i)} \\ &\quad - 2(\mathbf{\Phi} + N\hat{\mathbf{\Phi}}_{\theta})^{\mathrm{T}} \mathbf{\Sigma}_{\theta}^{-1} s \\ &\quad + (N+1) s^{\mathrm{T}} \mathbf{\Sigma}_{\theta}^{-1} s \\ &\quad + (N+1) \log |2\pi \mathbf{\Sigma}_{\theta}| \\ &= \mathbf{\Phi}^{\mathrm{T}} \mathbf{\Sigma}_{\theta}^{-1} \mathbf{\Phi} + \sum_{i=1}^{N} \mathbf{\Phi}_{\theta}^{(i)\mathrm{T}} \mathbf{\Sigma}_{\theta}^{-1} \mathbf{\Phi}_{\theta}^{(i)} \\ &\quad - \eta^{\mathrm{T}}(N+1) \mathbf{\Sigma}_{\theta}^{-1} \eta \\ &\quad + (N+1) [(s-\eta)^{\mathrm{T}} \mathbf{\Sigma}_{\theta}^{-1} (s-\eta) \\ &\quad + \log |2\pi \mathbf{\Sigma}_{\theta}|], \end{aligned} \quad (A3) $$

where we have recognized $\hat{\mathbf{\Phi}}_{\theta} = \frac{1}{N} \sum_{i=1}^{N} \mathbf{\Phi}_{\theta}^{(i)}$ (equation 12) and introduced $\eta = (\mathbf{\Phi} + N\hat{\mathbf{\Phi}}_{\theta})/(N+1)$. After integration over $s$, the last term gives a constant factor, so that $\mathcal{P}(\mathbf{\Phi}, \{\mathbf{\Phi}_{\theta}^{(i)}\}|\theta) \propto \exp[\hat{\ell}_2(\mathbf{\Phi}, \{\mathbf{\Phi}_{\theta}^{(i)}\})]$, with

$$ -2\hat{\ell}_2(\mathbf{\Phi}, \{\mathbf{\Phi}_{\theta}^{(i)}\}) = \mathbf{\Phi}^{\mathrm{T}} \mathbf{\Sigma}_{\theta}^{-1} \mathbf{\Phi} + \sum_{i=1}^{N} \mathbf{\Phi}_{\theta}^{(i)\mathrm{T}} \mathbf{\Sigma}_{\theta}^{-1} \mathbf{\Phi}_{\theta}^{(i)} - (\mathbf{\Phi} + N\hat{\mathbf{\Phi}}_{\theta})^{\mathrm{T}} \frac{1}{N+1} \mathbf{\Sigma}_{\theta}^{-1} (\mathbf{\Phi} + N\hat{\mathbf{\Phi}}_{\theta}). \quad (A4) $$

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

$$ -2\hat{\ell}_2(\mathbf{\Phi}, \{\mathbf{\Phi}_{\theta}^{(i)}\}) = (\mathbf{\Phi} - \hat{\mathbf{\Phi}}_{\theta})^{\mathrm{T}} \left( \frac{N+1}{N} \mathbf{\Sigma}_{\theta} \right)^{-1} (\mathbf{\Phi} - \hat{\mathbf{\Phi}}_{\theta}) + \text{constant terms.} \quad (A5) $$

In order to obtain a computable approximation of the likelihood, the unknown covariance $\mathbf{\Sigma}_{\theta}$ in $\hat{\ell}_2(\mathbf{\Phi}, \{\mathbf{\Phi}_{\theta}^{(i)}\})$ has to be approximated by $\hat{\mathbf{\Sigma}}_{\theta}$, defined by equation (14). The covariance of the effective likelihood is therefore $\hat{\mathbf{\Sigma}}'_{\theta} = \frac{N+1}{N} \hat{\mathbf{\Sigma}}_{\theta}$. The unknown inverse covariance $\mathbf{\Sigma}_{\theta}^{-1}$ is also replaced by its unbiased computable approximation $\hat{\mathbf{\Sigma}}_{\theta}^{-1}$, defined by equation (15). Finally, we use the normalization

condition $\int \mathcal{P}(\mathbf{\Phi}|\{\mathbf{\Phi}_{\theta}^{(i)}\}, \theta) d\mathbf{\Phi} = \int \frac{\mathcal{P}(\mathbf{\Phi}, \{\mathbf{\Phi}_{\theta}^{(i)}\}|\theta)}{\mathcal{P}(\mathbf{\Phi}|\{\mathbf{\Phi}_{\theta}^{(i)}\}|\theta)} d\mathbf{\Phi} = 1$ and evaluate at $\mathbf{\Phi} = \mathbf{\Phi}_O$, as prescribed by equation (7), to obtain $\hat{\ell}^N(\theta)$ given by equation (11). When $N \to \infty$, $\hat{\mathbf{\Phi}}_{\theta} \to s$ and $\hat{\mathbf{\Sigma}}'_{\theta} \to \mathbf{\Sigma}_{\theta}$, thus the limiting approximation is $\hat{\ell}(\theta)$ given by equation (16).

## APPENDIX B: DERIVATION OF THE EFFECTIVE POSTERIOR

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

$$ \begin{aligned} -2 \log \mathcal{P}(\mathbf{x}) &= \log |2\pi \mathbf{C}_0| + (\mathbf{x} - \mathbf{x}_0)^{\mathrm{T}} \mathbf{X}^{-1} (\mathbf{x} - \mathbf{x}_0) \\ &= \log |2\pi \mathbf{X}| + \xi_0^{\mathrm{T}} \mathbf{X} \xi_0 - 2 \xi_0^{\mathrm{T}} \mathbf{x} + \mathbf{x}^{\mathrm{T}} \mathbf{X}^{-1} \mathbf{x}, \end{aligned} \quad (B1) $$

where $\xi_0 = \mathbf{X}^{-1}\mathbf{x}_0$.

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

$$ \begin{aligned} -2\hat{\ell}^N(\theta) &= \log |2\pi \mathbf{C}_0| + [\mathbf{\Phi}_O - f_0 - \nabla f_0 \cdot (\theta - \theta_0)]^{\mathrm{T}} \mathbf{C}_0^{-1} [\mathbf{\Phi}_O - f_0 - \nabla f_0 \cdot (\theta - \theta_0)] \\ &= \log |2\pi \mathbf{C}_0| + (y_0 - \theta)^{\mathrm{T}} \mathbf{N}_0^{-1} (y_0 - \theta), \end{aligned} \quad (B2) $$

where we have defined

$$ \mathbf{N}_0 = [(\nabla f_0)^{\mathrm{T}} \mathbf{C}_0^{-1} \nabla f_0]^{-1} \quad (B3) $$

and

$$ y_0 = \theta_0 + (\nabla f_0)^{-1} \cdot (\mathbf{\Phi}_O - f_0), \quad (B4) $$

($\nabla f_0)^{-1}$ denoting the adjoint of the Jacobian characterizing the linearized black box (its computation will not be necessary). In canonical form, the Gaussian effective likelihood is written

$$ -2\hat{\ell}^N(\theta) = \log |2\pi \mathbf{C}_0| + \mu_0^{\mathrm{T}} \mathbf{N}_0 \mu_0 - 2\mu_0 \theta + \theta^{\mathrm{T}} \mathbf{N}_0^{-1} \theta, \quad (B5) $$

with $\mu_0 = \mathbf{N}_0^{-1} y_0 = \mathbf{N}_0^{-1} \theta_0 + (\nabla f_0)^{\mathrm{T}} \mathbf{C}_0^{-1} (\mathbf{\Phi}_O - f_0)$. Similarly, the prior (equation 23) is written

$$ -2 \log \mathcal{P}(\theta) = \log |2\pi \mathbf{S}| + \eta_0^{\mathrm{T}} \mathbf{S} \eta_0 - 2 \eta_0^{\mathrm{T}} \theta + \theta^{\mathrm{T}} \mathbf{S}^{-1} \theta, \quad (B6) $$

where $\eta_0 = \mathbf{S}^{-1}\theta_0$.

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

$$ \begin{aligned} -2 \log \mathcal{P}(\theta | \mathbf{\Phi})_{|\mathbf{\Phi}=\mathbf{\Phi}_O} &= -2(\mu_0 + \eta_0)^{\mathrm{T}} \theta \\ &\quad + \theta^{\mathrm{T}} (\mathbf{N}_0^{-1} + \mathbf{S}^{-1}) \theta \\ &\quad + \text{constant terms.} \end{aligned} \quad (B7) $$

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

$$ \begin{aligned} \gamma &= \mathbf{\Gamma}(\mu_0 + \eta_0) \\ &= \mathbf{\Gamma} \mathbf{N}_0^{-1} \theta_0 + \mathbf{\Gamma}(\nabla f_0)^{\mathrm{T}} \mathbf{C}_0^{-1} (\mathbf{\Phi}_O - f_0) + \mathbf{\Gamma} \mathbf{S}^{-1} \theta_0, \end{aligned} \quad (B8) $$

giving equation (25).

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