Primordial power spectrum and cosmology from black-box galaxy surveys
Florent Leclercq Wolfgang Enzi Jens Jasche Alan Heavens
























florent.leclercq@polytechnique.org
http://www.florent-leclercq.eu/

Imperial Centre for Inference and Cosmology (ICIC) & Astrophysics Group, Imperial College London, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, United Kingdom


enzi@mpa-garching.mpg.de


Max-Planck Institute for Astrophysics, Karl-Schwarzschild Strasse 1, D-85748 Garching, Germany


jens.jasche@fysik.su.se


The Oskar Klein Centre, Department of Physics, Stockholm University, Albanova University Center, SE 106 91 Stockholm, Sweden
Excellence Cluster Universe, Technische Universität München, Boltzmannstrasse 2, D-85748 Garching, Germany


a.heavens@imperial.ac.uk


Imperial Centre for Inference and Cosmology (ICIC) & Astrophysics Group, Imperial College London, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, United Kingdom











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 (), builds upon approximate Bayesian computation using a novel effective likelihood, and upon the linearisation of black-box models around an expansion point. Consequently, we obtain simple “filter equations” for an effective posterior of the primordial power spectrum, and a straightforward scheme for cosmological parameter inference. We demonstrate that the workload is computationally tractable, fixed a priori, and perfectly parallel. As a proof of concept, we apply our framework to a realistic synthetic galaxy survey, with a data model accounting for physical structure formation and incomplete and noisy galaxy observations. In doing so, we show that the use of non-linear numerical models allows the galaxy power spectrum to be safely fitted up to at least k_max = 0.5 h/Mpc, 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  and our data products publicly available at http://pyselfi.florent-leclercq.euhttp://pyselfi.florent-leclercq.eu.









Primordial power spectrum and cosmology from black-box galaxy surveys

    Alan Heavens
    December 30, 2023
=======================================================================












§ 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 <cit.>. The statistics of the initial density field are measured to be extremely close to Gaussian-distributed <cit.>. 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 <cit.>. A large variety of early-Universe models exhibit a deterministic relation between physical parameters of interest and the primordial matter power spectrum. The latter is therefore an interesting intermediate product for cosmological analyses, allowing parameter inference and model selection to be performed a posteriori without (or with minimal) loss of information. It can be seen as a largely agnostic and model-independent parametrisation of cosmological theories, which relies only on weak assumptions (isotropy and gaussianity).

For a long time, measuring the cosmological matter power spectrum has been one of the main goals of galaxy survey data analysis. However, inferring its shape accurately is a challenging task. Various systematic effects such as redshift uncertainties, complex survey geometries, selection effects, missing observations and foreground contamination can greatly hinder the measurement and analysis <cit.>. 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 <cit.>. 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 <cit.>. 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 <cit.>. 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 associated covariance matrices) and likelihood-based forward-modelling approaches (often associated with Bayesian statistics). Backward-modelling approaches suggest to directly account and correct for relevant effects in observational data as a pre-processing step. For example, the BAO “reconstruction” technique removes redshift-space distortions and corrects the density field from bulk motions using the inverse Zel'dovich approximation <cit.>. The end product is an estimator for the primordial matter power spectrum, which can be close to optimal if all relevant effects are modelled <cit.>. As backward-modelling approaches require substantial expert knowledge input, results are often largely model-dependent and with difficult propagation of uncertainties. Importantly, these approaches often rely on fiducial values for the parameters that are the target of the analysis. More recently, several thorough forward-modelling Bayesian approaches have been proposed to jointly infer the three-dimensional matter density field and its power spectrum from galaxy observations, while properly accounting for uncertainties and systematics <cit.>. Considerable effort is also put into reconstructing the primordial density field from present observations <cit.>. 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 <cit.>. More sophisticated approaches known as delfi <cit.> and bolfi <cit.> have also been recently introduced in cosmology. The ABC approach introduced in this work differs from all of the above in two aspects: (i) it allows the treatment of a much larger number of parameters (one hundred in this work), which correspond to primordial power spectrum amplitudes at different wavenumbers. To do so, (ii) it assumes the availability of a reasonably good guess of the target parameters, based on previous observations. This situation is fairly typical in cosmology, where the allowed space for parameters is already strongly constrained by previous experiments such as the Planck satellite <cit.>. Under these assumptions, we derive an effective likelihood for the problem. We use an expansion point in parameter space and linearise the black-box around it. The use of finite differencing to compute the gradient of the black-box makes evaluations of the effective likelihood computationally feasible. When further assuming that the prior is Gaussian, we find that the effective posterior distribution for the primordial power spectrum is a Gaussian with mean and covariance matrix given by two simple “filter equations” (equations (<ref>) and (<ref>)), which constitute the main result of this work. Finally, we show how to infer parameters of specific cosmological models using the linearised black-box. We propose to call the algorithm “simulator expansion for likelihood-free inference” ().

In order to illustrate the performance of the method, we apply it to a black-box which emulates realistic cosmological data. This black-box is built using Simbelmynë <cit.>, a hierarchical probabilistic simulator to generate synthetic galaxy survey data.[Simbelmynë is publicly available at https://bitbucket.org/florent-leclercq/simbelmyne/https://bitbucket.org/florent-leclercq/simbelmyne.] The data model involves a full cosmological N-body simulation (performed using our implementation of cola, ) to evolve the three-dimensional initial density field. It includes a treatment of galaxy bias, redshift-space distortions, survey geometry, selection effects, and instrumental noise. The statistical summary chosen is the estimated power spectrum of the galaxy number count field, as is standard in large-scale structure data analysis, but can be readily extended to include more information. As a result, the inferred primordial matter power spectrum is unbiased across the entire range of Fourier modes considered, and includes in particular BAOs, which were not included in the expansion point. Our analysis demonstrates that by using a fully numerical data model in conjunction with our statistical approach, one can safely fit the galaxy power spectrum even far in the non-linear regime, up to at least k_max≈ 0.5 h/Mpc, 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 organised as follows. In section <ref>, 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 <ref>, we describe the data-generating model used to test our method. The results obtained by combining the two are discussed in <ref>. We discuss the application of our method and prospects for cosmological data analysis, and provide our conclusions in section <ref>. Details of the statistical derivations are given in the appendices.



§ 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 <ref>, we design an effective likelihood for black-box models. In section <ref>, we exploit previous knowledge, as could have been obtained by earlier cosmological probes, in order to linearise the black-box. We discuss the parametrisation of the primordial matter power spectrum and its prior distribution in section <ref>. The equations for the effective posterior distribution are given in section <ref>. In section <ref>, we describe how to optimally choose the hyperparameters appearing in the prior. The inference of cosmological parameters from the linearised black-box is discussed in section <ref>.



 §.§ Design of an effective likelihood for black-box models





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



  §.§.§ 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 ∈ℝ^S parametrising the primordial power spectrum, and of nuisance parameters ∈ℝ^T, independent of . Nuisance parameters account for the entire stochasticity of the data model, such as initial phases, noise realisations, sample variance, etc. In our case, nuisance parameters will be all the random numbers generated by the galaxy survey simulator for initial conditions and instrumental noise; there are typically 𝒪(10^7) of those. Once realisations of  and  are specified, the output of the simulation d∈ℝ^D is a deterministic numerical function S, i.e. (d|, ) = _D(d - S(, )), where the symbol  denotes a probability distribution function (pdf) and _D a Dirac delta distribution. We refer to such realisations as mock observations. As usual in ABC approaches, the (often high-dimensional) raw prediction d can be compressed to a set of summary statistics Φ∈ℝ^P. We assume that this compression is a deterministic function C of d, i.e. (Φ |d) = _D(Φ-C(d)). It can be included in the model, so that the black-box is B≡C∘S and

    (Φ|, ) = _D(Φ - B(, )).

A graphical representation of the Bayesian hierarchical data model is presented in figure <ref>. 



  §.§.§ The exact Bayesian problem



Denoting by Φ_O the summary statistics of the observations, the inference problem considered is

    (|Φ)_|Φ=Φ_O = L() ()/Z_Φ,

where the likelihood is

    L() ≡(Φ|)_|Φ=Φ_O

and the normalisation constant is Z_Φ≡(Φ)_|Φ=Φ_O. By marginalising over  and using equation (<ref>), we have

    L()     =    ∫(Φ|, )_|Φ=Φ_O()  d
        =    ∫_D(Φ_O - B(, )) ()  d.

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



  §.§.§ 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 , we can generate an ensemble of N mock data realisations Φ_^(i) = B(, ^(i)) with i ∈{1, ... ,N}, by drawing independent and identically-distributed realisations of nuisance parameters  from the probability distribution () and evaluating the black-box. The approach then consists in explicitly conditioning all probabilities on {Φ_^(i)}. Using Bayes' theorem, we have:

    ( | Φ , {Φ_^(i)})     =    (Φ, {Φ_^(i)} | )  ()/(Φ, {Φ_^(i)})
        =    (Φ, {Φ_^(i)} |  )/( {Φ_^(i)})()/(Φ) .

Using for Φ the observed data Φ_O, we thus have the new inference problem

    (|Φ)_|Φ=Φ_O   ≈   ( | Φ , {Φ_^(i)})_|Φ=Φ_O
        =    L^N() ()/Z_Φ

where we have defined the first factor on the right-hand side of equation (<ref>) evaluated at Φ_O to be the effective likelihood (a computable approximation of the true likelihood):

    L^N() ≡(Φ, {Φ_^(i)} | )_|Φ=Φ_O/({Φ_^(i)}).


To arrive at a more explicit expression for the effective likelihood, we assume that observed data Φ_O and mock realisations Φ_^(i) are drawn from a common  (but unknown)  virtual signal s∈ℝ^P. The hierarchical representation of this effective data model is presented in figure <ref>. We assume that s carries the deterministic information on the target parameters , which is common to the data and to the 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 (Φ|s), from which data realisations are drawn independently once the virtual signal s is given. Since we do not want to infer s explicitly, marginalisation yields:

    (Φ, {Φ_^(i)} | )     =    ∫(Φ, {Φ_^(i)}, s |  )  ds
        =    ∫(Φ|s) ({Φ_^(i)} | s) (s|)  ds.

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 realisations. In this work, in absence of prior information on s, we set (s|) to a constant. The joint distribution of mock observations for a given virtual signal factorises, so that one can write

    (Φ, {Φ_^(i)} | ) ∝∫(Φ|s) [∏_n=1^N (Φ_^(i)| s)] ds.

In order to marginalise over s, we need to postulate a parametric form for the pdf (Φ|s). In this work, we assume a Gaussian distribution centred on the virtual signal s with a covariance matrix Σ_ quantifying the stochastic uncertainties inherent to the observations,[The investigation of different choices for (Φ|s) is left to future investigations.]

    -2 log(Φ|s) = log| 2πΣ_| + (Φ-s)^⊺Σ_^-1 (Φ-s).

The virtual signal can therefore be interpreted as the expectation E of Φ_ once  is specified: s = E[ Φ_]. We also have Σ_ = E[ (Φ_-s)(Φ_-s)^⊺].

Under these assumptions, marginalisation over s gives the effective likelihood as L^N() = exp[ ℓ̂^N() ] (the details of the computation are given in appendix <ref>), with:

    -2 ℓ̂^N() = log| 2πΣ̂_' | + (Φ_O - Φ̂_)^⊺Σ̂_'^-1 (Φ_O - Φ̂_),

where

    Φ̂_≡E^N [ Φ_] = 1/N∑_i=1^NΦ_^(i)

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

    Σ̂_' ≡N+1/NΣ̂_,   Σ̂_'^-1≡(N+1/N)^-1Σ̂_^-1,

where given a sufficiently large number of mock observations, a computable approximation of Σ_ is estimated in the following way:

    Σ_≈Σ̂_   ≡   E^N[ (Φ_ - Φ̂_)(Φ_ - Φ̂_)^⊺] 
        =    1/N-1∑_i=1^N (Φ_^(i) - Φ̂_) (Φ_^(i) - Φ̂_)^⊺ .

The effective likelihood also requires an estimator Σ̂_^-1 of the inverse covariance matrix Σ_^-1 (see equation (<ref>)). As argued by <cit.> and <cit.>, 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 <cit.> correction, which consists of replacing the true inverse covariance matrix by a scaled inverse sample covariance matrix:

    Σ_^-1≈Σ̂_^-1≡α( Σ̂_)^-1,   α≡N-P-2/N-1,

where P is the number of summary statistics. This correction debiases the expectation value ⟨( Σ̂_)^-1⟩ = α^-1Σ_^-1 of the estimator, under the assumption that Σ_^-1 is inverse-Wishart distributed.

The limiting approximation of L^N() when N →∞ is L() = (Φ|s)_|Φ = Φ_O = exp[ ℓ̃() ], with (as intended)

    -2 ℓ̃() ≡log| 2πΣ_| + (Φ_O-s)^⊺Σ_^-1 (Φ_O-s).


There are a number of interesting similarities and differences between the effective likelihood (equation (<ref>)) and the virtual signal pdf (equation (<ref>)). First, both are Gaussian distributions with respect to the observed data Φ_O. Second, the mean of the effective likelihood is given by the empirical average of simulated summaries, Φ̂_ = E^N [ Φ_]. This is an unbiased estimator of the virtual signal s = E[ Φ_] 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, ), but naturally appeared in the derivation (see appendix <ref>). Finally, the covariance of the effective likelihood is N+1/NΣ̂_, a multiple of the estimated covariance Σ̂_. The numerical prefactor, similar to Bessel's correction for the empirical sample variance, can be understood as follows. For a small number of simulations, the latent space associated to the virtual signal increases the observed scatter. For example, for N=1, the covariance to be used in equation (<ref>) is 2Σ̂_, which reflects the fact that observed and simulated data Φ_O and Φ_ can be drawn from opposite ends of the scatter around s of (Φ|s). However, when N →∞, N+1/NΣ̂_⟶Σ_, the intrinsic covariance of (Φ|s). This result is reasonable, since for large sample sizes the ensemble mean Φ̂_ converges to s, the mean of (Φ|s); in the same limit, we expect the covariance of the effective likelihood to approach Σ_, the covariance of (Φ|s).

Note that the scaling factor N+1/N directly arises from the presence of the virtual signal, as shown in appendix <ref>. Alternatively, we could have directly assumed a parametric form for the pdf (Φ|) without introducing the latent variable s. For instance, the synthetic likelihood directly amounts to postulating equation (<ref>), but without the scaling factor N+1/N for the estimated covariance matrix. This means that the end result for ℓ̂^N() has little sensitivity to our treatment and assumption for (s|). Furthermore, since N+1/N⟶ 1 when N →∞, our result with the virtual signal and the synthetic likelihood are equivalent, provided that the number of simulations is large enough.

Determining Φ̂_ and Σ̂_ requires N model evaluations per target parameters . In the next section, we show how the evaluation of this effective likelihood can be made more efficient when exploiting prior information on .



 §.§ Linearisation of black-box models



The numerical cost of the computation of the effective likelihood described in section <ref> 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  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 <cit.>. 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 Δ around an expansion point _0. Consequently, the target parameters are given by = _0 + Δ. Assuming that Φ̂_ is differentiable with respect to , we perform a first-order Taylor expansion in Δ around _0:

    Φ̂_≈f_0 + ∇f_0 · (-_0) ≡f() ,

where the defined function f is a linearised version of the averaged black-box. The first term corresponds to the mean mock observations at the expansion point _0, i.e. f_0 ≡Φ̂__0, while the second term involves ∇f_0, a P ×  S matrix corresponding to the gradient of mean mock observations at _0, whose components are (∇f_0)_ps≡∂Φ̂__0 p/∂_s. In this work, we estimate the gradient ∇f_0 via finite differencing: given a small step size h, each column is approximated by

    (∇f_0)^⊺_s ≈f(_s)-f_0/h,   with  _s ≡_0 + h   (_K^ss')_s' ∈ 1, S,

where _K is the Kronecker delta. To avoid obtaining a noisy gradient, nuisance parameters  are kept at the same values in the simulations used to compute the f(_s) and f_0. Further, we neglect the dependence of the covariance matrix on  and we estimate it at the expansion point according to equation (<ref>), i.e. Σ̂_' ≈Σ̂__0' ≡C_0. Importantly, fully characterising f under these assumptions only requires evaluations of the simulator, thus ensuring that the data model remains a black-box.

Using the linearised data model f as a proxy for the data model described in section <ref> simplifies the expression of the effective likelihood (equation (<ref>)) to

    -2 ℓ̂^N() ≈log| 2πC_0 | + [Φ_O - f()]^⊺C_0^-1[Φ_O - f()].


Evaluating f_0 and C_0 requires N_0 model evaluations at the expansion point _0. The computation of the gradient further requires N_s × S model evaluations. Therefore, the linearised data model is fully characterised by a fixed total of N_0 + N_s × 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 _0, a minimum would be N_0 ≥ 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(_s), N_s can in principle be smaller than P (and as small as 1); however N_s ≳ P will yield a safer evaluation of ∇f_0.

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



 §.§ The power spectrum prior distribution



In this paper, we aim at inferring the primordial matter power spectrum P(k), which is a continuous function of wavenumber k. We parametrise it by its amplitudes at a sufficient number S of support wavenumbers k_s. As we are particularly interested in BAOs, we fix a “wiggle-less” power spectrum P_0(k) and we work with the “wiggle function” θ(k) ≡ P(k)/P_0(k) as target function. Formally, with P_0 ∈ℝ^S the vector of components (P_0)_s ≡ P_0(k_s), the inference variable  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 . In this work, we include into our prior the following assumptions:


  * the power spectrum is Gaussian-distributed,

  * it is strongly constrained to live close to P_0,

  * it is a smooth function of wavenumber,

  * and the power spectrum P_0 is subject to cosmic variance.

It follows from assumptions 1 and 2 that () shall be a Gaussian distribution with mean _0 ≡1_ℝ^S, which is also to be used as the expansion point.

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

    (K)_ss'≡exp[ -1/2(k_s - k_s'/k_corr)^2 ].

The hyperparameter k_corr determines the length scale on which power spectrum amplitudes of different wavenumber correlate with each other. A large value of k_corr corresponds to a strong correlation between _s and _s', even if their corresponding wavenumbers are far from each other. On the other hand, for k_corr≪ k_s, k_s', K becomes the identity matrix and _s and _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_corr. In a realistic scenario, the prior covariance is not scale-independent as is K ((K)_ss = 1 for all s). We rather want the standard deviations (S)^1/2_ss to account for the cosmic variance affecting the power spectrum P_0 (and hence the mean _0). In terms of power spectrum amplitudes, cosmic variance at a scale k is given by P(k)^2/N_k, where N_k ∝ 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 , the coefficients K shall be multiplied by the coefficients of uu^⊺, where

    (u)_s ≡ 1 + σ_s = 1+α_cv/k_s^3/2

and α_cv is a hyperparameter characterising the “strength” of cosmic variance given the considered volume. Finally, the amplitude of the covariance matrix S can be captured by an overall scaling θ_norm^2. The final expression for the prior covariance matrix is therefore

    S≡θ_norm^2  uu^⊺∘K,

where ∘ is the Hadamard product. The standard deviations on the diagonal are (S)^1/2_ss = θ_norm (1 + σ_s), as intended.

The resulting prior on  is characterised by a set of three hyperparameters {k_corr, α_cv, θ_norm}, and given as

    -2log() ≡log| 2πS| + (-_0)^⊺S^-1 ( - _0) .

As demonstrated in section <ref>, 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.



 §.§ The power spectrum effective posterior distribution





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

    -2log(|Φ)_|Φ=Φ_O = log| 2πΓ| + (-)^⊺Γ^-1 (-),

with mean

    ≡_0 + Γ  (∇f_0)^⊺ C_0^-1 (Φ_O-f_0),

and covariance matrix

    Γ≡[ (∇f_0)^⊺ C_0^-1∇f_0 + S^-1]^-1.

The proof of this result (detailed in appendix <ref>) uses the same algebra as the derivation of the Wiener filter <cit.>. Equations (<ref>) and (<ref>) are the main result of this work. Equation (<ref>) 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 (<ref>). Table <ref> summarises the variables appearing in equations (<ref>) and (<ref>) and their interpretation in the context of galaxy survey data analysis.



 §.§ Optimisation of the prior hyperparameters



As discussed in section <ref>, the chosen prior distribution involves three hyperparameters {k_corr, α_cv, θ_norm}. By definition, α_cv characterises the strength of cosmic variance in the considered cosmological volume, such that the number of modes at a given scale k is N_k = k^3/α_cv^2. α_cv can therefore simply be measured for a given simulator setup (i.e. box size and mesh). On the contrary, k_corr and θ_norm are free hyperparameters and the power spectrum reconstruction ( and Γ) generally depends on their values. In this section, we propose a procedure to find optimal values for k_corr and θ_norm.

k_corr indicates the a priori smoothness of reconstructed wiggle functions θ(k) and θ_norm how much they can deviate from the expansion point. Together, these two parameters characterise the functional shapes of allowed target functions (much like kernels and their hyperparameters in techniques such as Gaussian process regression, see e.g. ). In analogy with hyperparameter optimisation in machine learning, we propose to optimise k_corr and θ_norm to reproduce the shape of a fiducial wiggle function _fid, using as likelihood the effective posterior distribution derived in section <ref>. More precisely, the likelihood for k_corr and θ_norm is defined by

    -2log(k_corr,θ_norm|_fid) ≡   log| 2πΓ| 
        + (_fid-)^⊺Γ^-1 (_fid-),

where  and Γ, defined by equations (<ref>) and (<ref>), are functions of k_corr and θ_norm (through S). At this point, it is of course possible to include a hyperprior on (k_corr,θ_norm), if desired. The maximum likelihood estimator (or maximum a posteriori estimator) then provides the optimal values of k_corr and θ_norm to be used to infer , its functional shape being assumed to be that of _fid.

Note that evaluating the likelihood given in equation (<ref>) is cheap once the linearised black-box f has been computed: no additional data model evaluation is required, only low-dimensional matrix operations. This allows the optimal values of k_corr and θ_norm to be found using standard optimisers.



 §.§ From the power spectrum to cosmological parameters





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

    (|) = _D(  - T() ),

where T is a deterministic function of cosmological parameters, typically a Boltzmann solver or a fitting function. Given this assumption, we have (Φ|) = (Φ|T()), and the inference of cosmological parameters gives

    (|Φ)_|Φ=Φ_O = L_^N() ()/Z_Φ ,

where the likelihood for cosmological parameters is defined as

    L_^N()    ≡   (Φ|T())_|Φ=Φ_O
        =    L^N(T()) .

Noting L_^N() ≡exp[ ℓ̂_^N() ] and using the linearised data model, we have (see equation (<ref>))

    -2 ℓ̂_^N() ≈   log| 2πC_0 | 
        + [Φ_O - f(T())]^⊺C_0^-1[Φ_O - f(T())].


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 <cit.> or kernel density estimates have been introduced. We note that our ABC technique does not suffer from this complication: the likelihood for cosmological parameters (equation (<ref>)) is a Gaussian centered on the data Φ_O with a fixed covariance matrix C_0. Furthermore, getting predictions for Φ amounts to evaluating f(T()), which does not require additional black-box evaluations once f is known. The Bayesian problem of inferring cosmological parameters  (equation (<ref>)) can therefore easily be solved by standard Markov Chain Monte Carlo techniques.



§ DATA MODEL



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



 §.§ Primordial power spectrum parametrisation



Throughout this paper, we work with a cubic equidistant grid with comoving side length of 1 Gpc/h and 256^3 voxels, spanning scales between k_s,min = 6.28 × 10^-3 h/Mpc and k_s,max = 1.4 h/Mpc. We use S=100 support wavenumbers k_s. The first eight are fixed to the values required by the Fourier grid given our setup. The remaining support wavenumbers are logarithmically spaced up to k_max. Any vector  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.



 §.§ 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 ΛCDM model is assumed, and fiducial cosmological parameters used are the Planck 2015 values <cit.>, given in table <ref>. These are used whenever the distance-redshift relation is needed, as well as for the gravitational evolution.[In principle, our approach would require treating cosmological parameters as nuisance parameters in the inference of the primordial power spectrum, and to marginalise over them. In this work, for simplicity, we keep them fixed to the values used to predict the expansion point.]


Simbelmynë first generates a realisation of the initial density contrast via the convolution approach <cit.>. 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

    δ_k^i≡√(P(k))  w_k .

Transformation back to configuration space yields a initial density contrast field ^i. One realisation of such an initial density field is shown in the left panel of figure <ref>.



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

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

    1+z_obs = (1+z_cosmo)(1+z_pec), with  z_pec≡ -v_r/c,

where z_cosmo is the true cosmological redshift, z_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 <cit.> to give the final density contrast field ^f. One realisation of the final redshift-space density is shown in the middle panel of figure <ref>.

The galaxy density ^g is predicted using a linear bias model, used in various previous studies <cit.>, such that in any cell x,

    ρ_x^g = N̅  (1 + b  δ_x^f) .

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

The last step corresponds to a virtual observation of the galaxy field, accounting for observational effects expected in actual surveys. To do so, we use the three-dimensional survey response operator (or window) W, consisting of the product of the radial selection function R(r) and the angular survey mask and completeness function C(n̂) for any line-of-sight n̂, i.e. W(n̂,r) ≡ R(r)   C(n̂). This operator accounts for the fact that we are only looking at certain parts of the sky and that we have different detection probabilities of galaxies depending on their distance. We obtained a simple model for W on a grid matching our simulations (256^3 voxels covering a volume of (1 Gpc/h)^3) 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 <cit.> with previously-published parameters for the r-band: α=-1.05 and M_∗ = -20.44 <cit.>, 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 figure <ref>).



We emulate the survey via galaxy number counts N_x^g in each cell x of the box. In order to account for instrumental noise, our model is a non-uniform Gaussian process, meaning that in every cell N_x^g is drawn from a Gaussian distribution with mean μ^g_x and standard deviation σ_x,

    N_x^g↶G(μ^g_x|σ_x)

The mean μ^g_x ≡ W_x ρ^g_x = W_x N̅ (1+ b δ^f_x) characterises the expected number of galaxies. It accounts for all physical effects via δ^f_x and for the survey response operator W_x. The standard deviation is defined as σ_x ≡σ√(W_xN̅), where we set the overall noise level σ≡ 10^-1. A slice through one realisation of the galaxy number count field N^g is shown in the right panel of figure <ref>.



 §.§ 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^f(k_r) ≡ C ×∑_|k| ∈ k_r| N_k^g|^2/N_k_r-2 .

The overall constant factor C = 1000^3/256^2 (Mpc/h)^3 arises from our Fourier transform convention; N_k_r represents the number of modes within the shell k_r; and the factor -2 arises from the assumption that the data power spectrum is inverse-Γ distributed with shape parameter N_k_r/2 and scale parameter C ×| N^g_k |^2/2 <cit.>.

As the galaxy data model relies on a full cosmological simulation, the predictions for the final power spectrum remain reasonable even at scales that have experienced substantial non-linearity. In spite of the approximations made, we trust our data model at the percent level up to k = 0.5 h/Mpc. For this reason, we use P=43 k_r-bins in the range [0.02,0.5] h/Mpc, ensuring that each bin contains at least 100 modes. These bins are logarithmically spaced for k_r ≥ 0.04 h/Mpc. For convenience, we normalise the output of the black-box using the expansion point, so that Φ is the P-dimensional vector of components A × 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 realisations of Φ for a given primordial power spectrum P(k) and specific realisations of initial phases and noise.

Due to the breakdown of models based on perturbation theory in the non-linear regime and to the difficulties in incorporating the impact of small-scale observational processes, state-of-the-art large-scale structure analyses are typically limited to k_max≲ 0.3 h/Mpc <cit.>. Pushing the analysis to k_max = 0.5 h/Mpc represents an increase by a factor of ∼ 5 in the number of modes used (scaling as k_max^3), which is expected to yield substantial improvements in the inference results.



 §.§ Idealised data model



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



§ RESULTS



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

We use for P_0(k) the “wiggle-less” BBKS power spectrum <cit.> under Planck 2015 cosmology (see table <ref>). Unknown ground truth cosmological parameters _gt are drawn from the (marginalised, Gaussian) Planck priors:

    [   h; Ω_b; Ω_m; n_S; σ_8 ]↶G[
    [  0.6774; 0.04860;  0.3089;  0.9667;  0.8159;         ],
    diag[  0.0046^2; 0.00030^2;  0.0062^2;  0.0040^2;  0.0086^2 ]].

The “wiggly” ground truth power spectrum P_gt(k) is generated with the <cit.> (EH) fitting function, using these cosmological parameters. It is used to simulate observed data Φ_O, with unknown nuisance parameters (phase realisation 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 ()_s ≡ P(k_s)/P_0(k_s) are the values of the wiggle function at the S=100 support wavenumbers defined in section <ref>. We note _gt and _fid the vectors of component P_gt(k_s)/P_0(k_s) and P_fid(k_s)/P_0(k_s), respectively.



 §.§ Diagnostics of the black-box







We created an ensemble of N_0=150 mock realisations at the expansion point _0 = 1_ℝ^S using different nuisance parameters. These are used to compute the average black-box at the expansion point, f_0 ≡Φ̂__0, and the covariance matrix C_0 ≡Σ̂'__0, using their definitions (equations (<ref>), (<ref>), and (<ref>)). The results are shown in figure <ref>. There, the left panel shows individual realisations Φ__0 and the observed data vector Φ_O. The average black-box f_0 is also plotted, with a credible region corresponding to two standard deviations (i.e. 2 √(diag(C_0))). The full estimated covariance matrix C_0 is shown in the right panel. As expected, the measured variance is larger on large scales due to cosmic variance, with some anti-correlations between pairs of bins. The effect of the mask, which increases power at the largest scales found in the simulation box, is also clearly visible.

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

As discussed in section <ref>, the linearised black-box f() is fully characterised by f_0, C_0 and ∇f_0. In this work, we used a total of N_0 + N_s × S = 10,150 simulations to get very precise estimates of C_0 and ∇f_0, although using fewer would have been possible.



 §.§ The prior and its optimisation







As discussed in section <ref>, we choose a Gaussian prior centered on the expansion point _0 with a covariance matrix S given by equation (<ref>) and characterised by three hyperparameters {k_corr, α_cv, θ_norm}.

Following the method presented in section <ref>, we found optimal values for the prior hyperparameters. The strength of cosmic variance within our simulation volume shall satisfy α_cv = √(k^3/N_k) at all scales k, where N_k is the number of modes. In our Fourier grid (described in section <ref>), we measured up to the Nyquist frequency α_cv = 8.848 × 10^-4, value that we adopt.

Assuming that the target function  follows the functional shape of the fiducial wiggle function _fid calculated with Planck cosmology, the likelihood for k_corr and θ_norm is given by equation (<ref>). We further assume broad, uncorrelated Gaussian hyperpriors on k_corr and θ_norm: k_corr∼G(0.020, 0.015^2) [h/Mpc] and θ_norm∼G(0.2, 0.3^2). The posterior surface is plotted in figure <ref>. We found the maximum a posteriori values using the popular optimiser L-BFGS <cit.>. In our run, these are k_corr = 0.0156 h/Mpc and θ_norm = 0.0478. The resulting prior covariance matrix is shown in figure <ref>. For the idealised data model (see section <ref>), the same procedure is applied and the optimal parameters are found to be k_corr = 0.0158 h/Mpc and θ_norm = 0.0535.



 §.§ The effective posterior







Using the effective likelihood, characterised by the observed data Φ_O and the linearised black-box described in section <ref>, as well as the optimised prior discussed in section <ref>, we obtained the effective posterior on . It is a Gaussian, with mean and covariance matrix given by the “filter equations” (<ref>) and (<ref>).

Figure <ref> shows the inferred primordial wiggle function  in comparison with the expansion point _0, the ground truth _gt, and the fiducial function used to optimise prior hyperparameters, _fid. 2σ credible regions are shown for the prior and the posterior (i.e. 2 √(diag(S)) and 2 √(diag(Γ)), respectively). The top panel corresponds to the result obtained using the idealised data model (a Gaussian random field, see section <ref>), 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σ credible intervals of the reconstruction. As discussed in the introduction, this effective posterior can be seen as a largely model-independent parametrisation of the theory, containing all the available cosmological information under weak assumptions.

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

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

Since the proposed method is fully Bayesian, we do not simply obtain a point estimate, but a complete probability distribution, which provides a detailed quantification of uncertainties. Figure <ref> shows the covariance matrix Γ of the Gaussian effective posterior; it can be compared with the prior covariance matrix S shown in figure <ref>.



 §.§ Cosmological parameters





The last step in the analysis is to infer parameters  of specific cosmological models, given the observed data Φ_O. In this section, we assume a flat ΛCDM model, characterised by 5 parameters {h, Ω_b, Ω_m, n_S, σ_8 }. For simplicity, our prescription to generate primordial power spectra given cosmological parameters is the EH fitting function. The generative process T (from  to ) is therefore the EH fitting function divided by P_0, sampled at the support wavenumbers k_s. It could easily be generalised to extensions of the flat ΛCDM model and to include a Boltzmann solver. The linearised data model f (from  to Φ) has already been characterised in section <ref>.

The effective likelihood L_^N() for cosmological parameters is given by equation (<ref>). Consistently with the expansion point used to linearise the black-box, we complement L_^N() with a Gaussian prior () centred on Planck cosmological parameters, but with broader variance: the diagonal covariance matrix given in equation (<ref>) is multiplied by a factor of 3. We explored the effective posterior (|Φ)_|Φ=Φ_O (equation (<ref>)) via MCMC (performed using the emcee code, ), ensuring sufficient convergence.

Results are shown in figure <ref>. Prior contours are shown in blue, and posterior contours are shown in red and purple for two different realisations of the data Φ_O. The two data realisations have been generated using different ground truth cosmological parameters _gt (shown as dashed and dotted lines, respectively), as well as different nuisance parameters (phase realisation and noise). The plot demonstrates that cosmological parameter inference is unbiased and robust to nuisances imprinted in the data.



§ 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, , which can be used to infer the primordial matter power spectrum and cosmological parameters from galaxy surveys. The main results are the “filter equations” (<ref>) and (<ref>). They can be applied to get an effective posterior for any model where the mean and covariance of the data are estimated from arbitrarily complex forward models. Essentially everything is obtained from a simulator, which can be treated as a black-box, without necessity to include any knowledge of its internal mechanisms into the statistical analysis. 

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

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


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

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

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

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

 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,  allows the likelihood-free inference of S ≳ 100 parameters, as is necessary for a model-independent parametrisation of theory in cosmology.

In this work, we demonstrated that a “non-wiggly” expansion point _0 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 linearised black-box, used to obtain a new posterior. Using a sufficient number of iterations, we expect the effective posterior to converge to the true function, even if it strongly deviates from the first expansion point. We leave the detailed investigation of this idea to future studies.

In this paper, we showed a successful application of finite differencing to obtain the gradient of the averaged black-box ∇f_0, but our equations could be used with other techniques, such as automatic differentiation. The data covariance matrix at the expansion point C_0 also needs to be evaluated; for this task, and for certain summary statistics, variance reduction techniques such as the use of fixed and paired simulations <cit.> or hybrid estimators <cit.> could further be exploited.

As a proof of concept, we applied our technique in conjunction with the artificial galaxy survey simulator Simbelmynë, emulating relevant effects at play: non-linear gravitational structure formation, redshift-space distortions, a survey mask and selection function, and instrumental noise. As a result, the inferred primordial power spectrum is unbiased with a distinct identification of BAO wiggles, even far in the Silk damping tail. We also demonstrated that unbiased inference of cosmological parameters is possible. In spite of the non-linear evolution of structures on small scales, we are able to use the power spectrum of the galaxy field as summary statistic, up to k_max = 0.5 h/Mpc. 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 ∼√(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, , which we publicly released, together with documentation and the data necessary to reproduce the results of the present paper.[Currently, the code's homepage is hosted at http://pyselfi.florent-leclercq.eu/http://pyselfi.florent-leclercq.eu; the sources are available on GitHub at https://github.com/florent-leclercq/pyselfi/https://github.com/florent-leclercq/pyselfi; and the documentation is on Read the Docs at https://pyselfi.readthedocs.io/en/latest/https://pyselfi.readthedocs.io.] The application of this method to more complex models and to real survey data is left for future research.

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



§ DERIVATION OF THE EFFECTIVE LIKELIHOOD



In this appendix, we derive the approximate likelihood given in equation (<ref>), starting from equation (<ref>) and the assumptions detailed in section <ref>. We have

    (Φ, {Φ_^(i)} | ) ∝∫(Φ|s) [∏_n=1^N (Φ_^(i)| s)] ds.

Using the parametric form for (Φ|s) (equation (<ref>)) yields (Φ, {Φ_^(i)} | ) ∝∫exp[ ℓ̂_1(Φ, {Φ_^(i)}, s) ]  ds, with

    -2ℓ̂_1(Φ, {Φ_^(i)}, s)    ≡    (Φ-s)^⊺Σ_^-1(Φ-s) 
           + ∑_i=1^N (Φ_^(i)-s)^⊺Σ_^-1(Φ_^(i)-s) 
           + (N+1) log |2πΣ_| .

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

    -2ℓ̂_1(Φ, {Φ_^(i)}, s)     =    Φ^⊺Σ_^-1Φ + ∑_i=1^N Φ_^(i)⊺Σ_^-1Φ_^(i)
           - 2(Φ + NΦ̂_)^⊺Σ_^-1s
           + (N+1)  s^⊺Σ_^-1s
           + (N+1) log |2πΣ_|  
        =    Φ^⊺Σ_^-1Φ + ∑_i=1^N Φ_^(i)⊺Σ_^-1Φ_^(i)
           - ^⊺ (N+1)  Σ_^-1
           + (N+1)[(s - )^⊺ Σ_^-1 (s-) 
             + log |2πΣ_|],

where we have recognised Φ̂_ = 1/N∑_i=1^N Φ_^(i) (equation (<ref>)) and introduced ≡ (Φ + NΦ̂_)/(N+1). After integration over s, the last term gives a constant factor, so that (Φ, {Φ_^(i)} | ) ∝exp[ ℓ̂_2(Φ, {Φ_^(i)}) ], with

    -2ℓ̂_2(Φ, {Φ_^(i)})    ≡   Φ^⊺Σ_^-1Φ + ∑_i=1^N Φ_^(i)⊺Σ_^-1Φ_^(i)
           -(Φ+NΦ̂_)^⊺1/N+1 Σ_^-1 (Φ+NΦ̂_) .

We now complete the square with respect to Φ to obtain

    -2ℓ̂_2(Φ, {Φ_^(i)})     =     (Φ - Φ̂_)^⊺( N+1/NΣ_)^-1  (Φ - Φ̂_) 
           + constant terms.

In order to obtain a computable approximation of the likelihood, the unknown covariance Σ_ in ℓ̂_2(Φ, {Φ_^(i)}) has to be approximated by Σ̂_, defined by equation (<ref>). The covariance of the effective likelihood is therefore Σ̂_' ≡N+1/NΣ̂_. The unknown inverse covariance Σ_^-1 is also replaced by its unbiased computable approximation Σ̂_^-1, defined by equation (<ref>). Finally, we use the normalisation condition ∫(Φ | {Φ_^(i)} , )  dΦ = ∫(Φ , {Φ_^(i)} | )/({Φ_^(i)} | ) dΦ = 1 and evaluate at Φ=Φ_O, as prescribed by equation (<ref>), to obtain ℓ̂^N() given by equation (<ref>). When N →∞, Φ̂_⟶s and Σ̂_' ⟶Σ_, thus the limiting approximation is ℓ̃() given by equation (<ref>).



§ DERIVATION OF THE EFFECTIVE POSTERIOR



We recall the canonical form of the Gaussian distribution with mean x_0 and covariance matrix X, given as

    -2log(x)     =    log| 2πX| + (x-x_0)^⊺X^-1 (x-x_0) 
        =    log| 2πX| + _0^⊺X_0 -2 _0^⊺x + x^⊺X^-1x,

where _0 ≡X^-1x_0.

Using the linearised data model (equation (<ref>)) in the expression of the effective likelihood (equation (<ref>)), we get

    -2ℓ̂^N()     =    log| 2πC_0 | + [ Φ_O - f_0 - ∇f_0 · (-_0 )]^⊺·
                  C_0^-1[ Φ_O - f_0 - ∇f_0 · (-_0 )] 
        =    log| 2πC_0 | + (y_0 - )^⊺N_0^-1 (y_0 - ),

where we have defined

    N_0 ≡[ (∇f_0)^⊺C_0^-1∇f_0 ]^-1

and

    y_0 ≡_0 + (∇f_0)^-1· (Φ_O-f_0),

(∇f_0)^-1 denoting the adjoint of the Jacobian characterising the linearised black-box (its computation will not be necessary). In canonical form, the Gaussian effective likelihood is written

    -2ℓ̂^N() = log| 2πC_0 | + _0^⊺N_0 _0 - 2_0  + ^⊺N_0^-1,

with _0 ≡N_0^-1y_0 = N_0^-1_0 + (∇f_0)^⊺C_0^-1 (Φ_O-f_0).
Similarly, the prior (equation (<ref>)) is written

    -2log() = log| 2πS| + _0^⊺S_0 -2 _0^⊺ + ^⊺S^-1

where _0 ≡S^-1_0.

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

    -2log(|Φ)_|Φ=Φ_O    =     -2(_0+_0)^⊺
           + ^⊺ (N_0^-1 + S^-1) 
           + constant terms.

This is the canonical form of a Gaussian distribution, where the covariance matrix is identified as Γ≡ (N_0^-1 + S^-1)^-1, giving equation (<ref>), and the mean is identified as

    =    Γ(_0+_0) 
        =    ΓN_0^-1_0 + Γ(∇f_0)^⊺C_0^-1 (Φ_O-f_0) + ΓS^-1_0,

giving equation (<ref>).

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 s̅ and signal covariance S, and a likelihood with mean d̅ and noise covariance N, the filter covariance is (N^-1+S^-1)^-1 and the filtered signal is s̅ + (N^-1+S^-1)^-1N^-1(d_O-d̅).



§ STATEMENT OF CONTRIBUTION


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

FL is grateful to Guilhem Lavaux and Andrew Jaffe for useful discussions. This work has made use of a modified version of https://pypi.org/project/pyGTC/pygtc <cit.>. Numerical computations were done on the cx1 cluster hosted by the Research Computing Service facilities at Imperial College London (http://doi.org/10.14469/hpc/2232doi:10.14469/hpc/2232). This work is done within the Aquila Consortium (https://aquila-consortium.orghttps://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.dewww.universe-cluster.de).





§ REFERENCES



