Towards optimal extraction of cosmological information from nonlinear data
Uros Seljak, Grigor Aslanyan, Yu Feng, Chirag Modi
TL;DR
This paper introduces a Bayesian method to optimally extract cosmological information from nonlinear data by reconstructing initial density modes, improving the estimation of the initial power spectrum and baryonic acoustic oscillations.
Contribution
It develops a Bayesian framework that expands the likelihood around the maximum posterior of linear modes, enabling near-optimal reconstruction of the initial power spectrum from nonlinear data.
Findings
Unbiased estimator of the initial power spectrum achieved.
Near-optimal reconstruction of baryonic acoustic oscillations demonstrated.
Method applicable using simplified 2-LPT or N-body simulations.
Abstract
One of the main unsolved problems of cosmology is how to maximize the extraction of information from nonlinear data. If the data are nonlinear the usual approach is to employ a sequence of statistics (N-point statistics, counting statistics of clusters, density peaks or voids etc.), along with the corresponding covariance matrices. However, this approach is computationally prohibitive and has not been shown to be exhaustive in terms of information content. Here we instead develop a Bayesian approach, expanding the likelihood around the maximum posterior of linear modes, which we solve for using optimization methods. By integrating out the modes using perturbative expansion of the likelihood we construct an initial power spectrum estimator, which for a fixed forward model contains all the cosmological information if the initial modes are gaussian distributed. We develop a method to…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Towards optimal extraction of cosmological information from nonlinear data
Uroš Seljak
Grigor Aslanyan
Yu Feng
Chirag Modi
Abstract
One of the main unsolved problems of cosmology is how to maximize the extraction of information from nonlinear data. If the data are nonlinear the usual approach is to employ a sequence of statistics (N-point statistics, counting statistics of clusters, density peaks or voids etc.), along with the corresponding covariance matrices. However, this approach is computationally prohibitive and has not been shown to be exhaustive in terms of information content. Here we instead develop a hierarchical Bayesian approach, expanding the likelihood around the maximum posterior of linear modes, which we solve for using optimization methods. By integrating out the modes using perturbative expansion of the likelihood we construct an initial power spectrum estimator, which for a fixed forward model contains all the cosmological information if the initial modes are gaussian distributed. We develop a method to construct the window and covariance matrix such that the estimator is explicitly unbiased and nearly optimal. We then generalize the method to include the forward model parameters, including cosmological and nuisance parameters, and primordial non-gaussianity. We apply the method in the simplified context of nonlinear structure formation, using either simplified 2-LPT dynamics or N-body simulations as the nonlinear mapping between linear and nonlinear density, and 2-LPT dynamics in the optimization steps used to reconstruct the initial density modes. We demonstrate that the method gives an unbiased estimator of the initial power spectrum, providing among other a near optimal reconstruction of linear baryonic acoustic oscillations.
1 Introduction
The issue of optimal map reconstruction and optimal power spectrum reconstruction from a set of noisy and sparsely sampled data has received a lot of attention in the field of large scale structure (LSS) and cosmic microwave background (CMB) anisotropies. Most of the work assumes the data are a linear transformation of the initial modes, as is the case with CMB [1, 2, 3]. In this case the minimum variance map solution is the well-known Wiener filter map [4]. Computing Wiener filter can be expensive, since it requires an inverse noise weighting of the data, where in the cosmological context noise consists both of the actual measurement noise and the actual signal (the so called sampling or cosmic variance noise). Noise covariance matrix is typically sparse (often diagonal) in the configuration space, and given sparsely sampled data its Fourier space representation is not diagonal. Signal covariance matrix is diagonal in Fourier space, and its configuration space representation is not diagonal. The sum of the two is thus not sparse in any basis, and computing the inverse covariance matrix using brute-force methods is an process, which becomes prohibitively expensive for large . Alternative methods must therefore be used to solve for the Wiener filter map [5, 6].
The reconstructed map typically has millions of data points and contains too much information to be useful on its own. What we want instead is the optimal power spectrum given the map. This is a typical hierarchical Bayesian setting, where there are many latent variables that need to be marginalized over, and only their priors remain. Here the latent variables are the modes and their prior is the power spectrum, which is a useful summary statistic, since for a gaussian field of initial modes it contains all the information present in the data, and the data compression is lossless. In the linear regime it can be computed using the optimal quadratic estimator, which is quadratic in the data. This requires the data to be first inverse covariance matrix weighted [1, 2, 3], hence it is also computationally expensive. The two problems, optimal map making and optimal power spectrum, are connected: the optimal quadratic estimator can be built out of the Wiener field map reconstruction [7].
An alternative method to optimal quadratic estimator is to use sampling methods to determine the power spectrum probability distribution. For example, a Gibbs sampling approach consists of a two step sampling procedure [8]. In the first one a sample map consistent with the data is created by adding a generalized noise realization to the Wiener filter map. In the second, sampling of the power spectrum is created consistent with the given map realization. Creating a map sample also requires an expensive inversion of the covariance matrix, and as a consequence these sampling approaches are typically slower than the optimal quadratic estimator, and converge particularly slowly in the high noise regime.
These methods, while optimal in the linear regime, are often replaced with faster and less optimal methods due to their computational cost. Traditional power spectrum estimations of CMB and LSS use a linear transformation between the data and the model (Fourier transform, FT), and assume that the covariance matrix of the data is diagonal (which allows a fast evaluation of the corresponding curvature matrix of the modes). In LSS this is the so called FKP method [9], while in CMB it is called the pseudo- method [10, 11]. This solution is not optimal, because the correct weighting of the data is to multiply it with the inverse of the covariance matrix, an evaluation that scales as and is too costly to be performed on large data sets. FKP or pseudo- weighting is particularly poor if the noise is varying significantly across the survey, in which case it is optimized for a single value of the power spectrum (the one used in FKP weighting). It also fails if there is a complicated geometry of the survey. Moreover, the method does not handle optimally if there are modes that are contaminated and create large scale correlations that need to be marginalized over. Nevertheless, these methods have become standard, because for linear problems one can show that under certain assumptions (typically valid on scales small compared to the survey size) the amount of information is nearly the same [12].
So far the discussion above has been about a linear mapping between the modes and the data. While the primary CMB is linear to a very good approximation, in LSS this is only valid on the largest scales. Most of the LSS information is contained on smaller, nonlinear, scales: in a 3-d LSS survey the number of modes scales as , where is the typical wavevector and the typical scale. The question of optimal linear map given the nonlinear data has been addressed in recent work, and Hamiltonian Monte Carlo (HMC) sampling approach has been used to create the map [13], as well as create samples consistent with the data for a given power spectrum [14, 15]. Hamiltonian sampling has an advantage over regular (e.g. Metropolis-Hastings) sampling that it can propagate far from the current sample using dynamics rather than random walk, and still the acceptance rate can be very high, both of which are necessary conditions for a rapid convergence in large dimensions due to the curse of dimensionality. However, in the current implementations the samples are still highly correlated, with correlation length of order 100-200 being reported [15]. To create each sample one needs many Hamiltonian dynamics evaluation steps (of order 10), so the total cost of evaluation of a single independent realization can easily exceed 1000 calls, each being a full forward model. It also requires knowledge of a gradient of the data model with respect to initial modes, a feature shared with the methods developed in this paper.
In contrast to the optimal linear power spectrum analyses, there has been very little work in terms of extracting the optimal summary statistics given a set of noisy and incomplete data in the nonlinear regime. For sampling methods this may be too difficult to solve, if HMC requires 1000 calls or more to create a single independent realization. In the high noise regime sampling converges very slowly onto the correct cosmological model, because most of the mode power comes from the assumed power spectrum at a previous step. As a consequence this approach may not be feasible. Most of the focus so far has instead been to extract information from the two-point function of the nonlinear modes. In this case the same quadratic estimator methods discussed above can be used [16], but they are not optimal since there is information in higher order statistics as well. Because of this the two point function analysis is sometimes supplemented with additional statistics that are most likely to be complementary. Among these are the higher order correlations, starting with the three point function/bispectrum, various void statistics, topology statistics, counts of objects like clusters or other density peaks, and reconstruction methods, primarily focused on baryonic acoustic oscillations (BAO) (see e.g. [17] for a review of different probes of LSS). These approaches share the property that they use statistical information in addition to the two point statistics of the final field, but beyond that they differ enormously in terms of their motivation and scope. Combining the different statistics creates a significant problem of modeling their joint error distribution, since their joint covariance matrix cannot be computed ab-initio, but must instead be obtained from the simulations, which are noisy and expensive. Moreover, there is no guarantee that even after combining several of these statistics one will exhaust the information content in the data: new statistics are continuously being proposed (and argued to be superior by their authors), suggesting that this question is unlikely to be settled by this ad-hoc approach. Ideally, what is needed is an ab-initio approach that is built with a guarantee to give a nearly optimal answer. This is an approach we attempt to develop in this paper.
Since the initial density modes are assumed to be gaussian, their power spectrum is a lossless summary statistic, which should contain all the information present in the data (with an exception of cosmological parameters that may be required to map from the initial modes into the final data). The lossless nature may however break down in the situations of shell crossings, where there can be more than one initial solution that maps to the same final data. In this paper we derive an initial linear power spectrum estimation, given incomplete and noisy data and given a nonlinear model between the initial modes and the data. We will show that as an intermediate step we will also need to derive the optimal initial density reconstruction. For the linear case the optimal power spectrum from the optimal density reconstruction has been derived in [7], and here we generalize the expressions to the nonlinear model. We will show that it is possible to cast the solution into an optimization problem, and we develop methods to solve it efficiently. We then generalize the method to include forward model parameters as well. As a proof of principle we apply the derived expressions to extract the initial power spectrum from a final density field obtained either in an N-body simulation of dark matter or in a 2-LPT simulation models, and using a 2-LPT approximation (e.g. [18]) as a forward model.
2 Statistical approach and heuristic derivation
Following the notation of [7] let us assume that we measure some nonlinear observations d(\hbox{\boldmath{r}}_{i}) at pixelized (2-d or 3-d) spatial positions \hbox{\boldmath{r}}_{i}. The data could be nonlinear dark matter density (although this is typically not directly observable), some projected dark matter density such as lensing shear or convergence, galaxy density or luminosity, Sunyaev-Zeldovich intensity etc. We arrange these into a vector \hbox{\boldmath{d}}=\{d(\hbox{\boldmath{r}}_{i})\}(i=1,...,N). Each measurement consists of a signal and a noise contribution, \hbox{\boldmath{d}}=\hbox{\boldmath{f(s,\lambda)}}+\hbox{\boldmath{d}}_{n}, where noise is assumed to be uncorrelated with the signal. Here is a nonlinear mapping from the initial linear density modes to the final model prediction of the observation, and \hbox{\boldmath{s}}=\{s_{j}\}(j=1,...,M) are the underlying initial density mode coefficients that we wish to estimate. In the statistical language are the observed variables, are the latent variables and are the forward model parameters. While our primary motivation for these is linear matter over-density , the formalism we develop is also useful for other applications such as CMB anisotropies, hence we keep the notation more general. We will assume these coefficients to be in Fourier space, where their gaussian prior can be written in a diagonal form. The modes are complex and obey s^{*}(\hbox{\boldmath{k}})=s(-\hbox{\boldmath{k}}), where is the wavevector, but in our labeling of modes we will treat real and imaginary component as two independent modes. If the mapping is linear we can write \hbox{\boldmath{f}}(\hbox{\boldmath{s}})=\hbox{\boldmath{f^{\prime}}}\hbox{\boldmath{s}}.
We will assume the nonlinear mapping to be computable given the initial density field and given some forward model parameters . Typically this mapping will be a full N-body simulation, with some additional processing to produce a realistic model prediction for the specific observations. The forward model parameters can be matter density, massive neutrinos or other parameters that affect the growth of structure, various astrophysical modeling parameters (related to how galaxies and baryons populate dark matter halos), and observational nuisance parameters (such as shear bias in weak lensing etc.). The mapping may also include smoothing (e.g. beam smoothing) or pixelization of the data. The true underlying field has an infinite number of Fourier modes, but only a finite number of these can be estimated. Typically we will embed a given LSS survey into a periodic box larger than the survey, with zero padding in regions without the data, but in this paper we will simplify this to periodic box analysis.
The modes are latent variables with a prior of their own. This prior is parametrized as a multivariate gaussian, with covariance matrix \hbox{\boldmath{S}}=\langle\hbox{\boldmath{ss}}^{{\dagger}}\rangle, which is assumed to be diagonal in Fourier space. We will assume the power spectrum depends on parameters , which will typically be bandpowers, but this will later be generalized to any parameters that change the power spectrum. These bandpowers can also have a prior (a hyperprior in the language of hierarchical Bayesian models), but in this paper we will assume their prior is flat, and we will not even impose positivity, since these are summary statistics expected to be used in a later analysis of cosmological parameters.
The noise vector \hbox{\boldmath{d}}_{n}=\{d_{n,i}\}(i=1,...,N) is parametrized with the noise covariance matrix \hbox{\boldmath{N}}=\langle\hbox{\boldmath{d}}_{n}\hbox{\boldmath{d}}_{n}^{{\dagger}}\rangle.The i-th diagonal element corresponds to the noise variance at the spatial position . This noise matrix is assumed to be known, uncorrelated with the signal and diagonal (or sparse) in real space, so that any operations involving noise matrix scale as rather than or steeper. For simplicity we will also assume noise is gaussian distributed, but this can be generalized to other probability distributions. It is reasonably appropriate for weak lensing, where by central limit theorem averaging over shape noise produces approximately gaussian noise. Outside the survey mask (which can include holes inside the survey) we will assume N(\hbox{\boldmath{r}}_{i})=\infty and assign d(\hbox{\boldmath{r}}_{i})=0 (although assigning any other value would be just as good).
The goal of parameter inference is to derive the posterior distribution of parameters and given their prior and the data . In some cases, such as linear model, one can write analytic expression for the likelihood L(\hbox{\boldmath{d}}|\hbox{\boldmath{\Theta}},\hbox{\boldmath{\lambda}}), but its explicit evaluation requires inversion and trace or determinant of a very large non-sparse matrix, both of which are operations, which becomes too expensive when the size of the data becomes large. In the nonlinear case even writing down the posterior is a major challenge. We will argue that it is easier to solve the problem in a typical hierarchical Bayesian model approach, by working with the latent variables and solving for these first, then marginalizing over them to obtain the likelihood and the posterior of parameters (which are the same if the priors are flat, as we will assume here). This approach perhaps seems counter-intuitive, since the dimension of the latent variables is comparable to the data . However, the major simplification in the linear case is that there are no large matrix inversions required for the solution to be obtained. For the nonlinear case we will argue that finding a good solution at or close to the global minimum is achievable even in a large number of dimensions, and the problem cannot even be solved in the absence of latent variables. By writing a full probabilistic model in terms of conditional probabilities between individual variables, in this case latent variables conditional on , and data conditional on and , we are able to solve for the posterior of and given the data (combining with the (hyper)priors on and ).
One can write the joint probability of signal and noise \hbox{\boldmath{d}}_{n} as a product of individual probabilities, under the assumption that they are uncorrelated,
[TABLE]
The Bayes theorem can be applied to , and to obtain their posterior
[TABLE]
The first term is the probability distribution of the noise for some set of observations , \hbox{\boldmath{d}}_{n}=\hbox{\boldmath{d}}-\hbox{\boldmath{f}}(\hbox{\boldmath{s,\lambda}}), while the second is the prior for given (or, equivalently, ). The last term P(\hbox{\boldmath{\Theta,\lambda}}) is the prior on the initial power spectrum parameters and on forward model parameters . We will assume this prior is flat, so that we have no prior information on them. This is because we want the result of the analysis to summarize the information from the given data, without inclusion of external data. This is easy to modify at a later stage if combining different data sets. Similarly, we will not be concerned with the normalization of the posterior, so we simply use the proportionality symbol in equation 2.2. We will fix the normalization at the end using gaussian approximation for the posterior.
Maximizing the posterior in terms of at a fixed and gives \hat{\hbox{\boldmath{s}}}, the so called maximal a posterior (MAP), which is also the Wiener filter map in the linear case. We will use optimization methods with analytic gradient evaluation to perform this step. Gradient based optimization methods can be very efficient at finding the maximum, but typically do not provide a reliable curvature matrix (Hessian or its inverse). We will develop a method where the curvature matrix of around \hat{\hbox{\boldmath{s}}} is not needed in the construction of the final solution, so that any optimization method can be used.
Ultimately the optimal map is not what we really care about. What we want is the posterior of model parameters and independent of latent variables . One way to do this would be to use sampling of both , and , but with dimensions this can be extremely expensive. Another approach would be to maximize the posterior with respect to all the parameters at once. This approach is strictly not valid when the latent variables depend on the parameters we wish to determine: maximization is not the same as marginalization. In addition, this approach has technical difficulties related to the fact that for the parameters we care about we also need a very accurate curvature matrix and its log determinant, which is not provided by the optimizer. In this paper we instead perform an analytic marginalization over the latent variables.
In our approach, we work with large dimensions of the modes, . We will assume that the same maximum is always found. To be more precise, we will argue that many posterior maxima can exist and do not change the nature of the solution as long as they are all similar to each other in a well defined sense discussed in section 4). This assumption can be explicitly proven in the linear regime where the loss function is convex [7]. In the nonlinear case it depends on the nature of the problem, but cannot be proven generally. For example, in the strongly nonlinear regime, inside virialized regions, satellites on orbits can yield several identical phase space configurations. While this degeneracy will be broken by the gaussian prior on the modes, it still suggests local maxima in posterior (even if not of equal height). We will proceed not ignoring this issue, but assuming it can be calibrated out under the assumption that the maxima are always quantitatively similar.
We will thus analytically integrate over
[TABLE]
performing this marginalization around the MAP solution \hat{\hbox{\boldmath{s}}}. We will perform a perturbative expansion of the posterior around MAP, which we can integrate over. This produces two terms, one is simply the value of the loss function in the exponent of equation 2.2 at the MAP position, and the second is the volume of posterior, roughly given by the determinant of the curvature matrix at that position (and further increased by the higher order moments). However, we will never explicitly perform this integral, and instead develop a method to determine the resulting volume element using simulations, so in this sense our approach does not rely on perturbative expansion.
Once we have P(\hbox{\boldmath{\Theta,\lambda}}|\hbox{\boldmath{d}}) we can perform its maximization with respect to to find their peak posterior solution, i.e. the maximum likelihood (since we assume flat prior). By expanding the log posterior to second order we also obtain the curvature matrix under the Laplace approximation. We will use Newton’s method to find the maximum likelihood solution and then construct an explicitly unbiased estimator out of it. In this paper we focus primarily on bandpowers for , and we argue that quadratic expansion of their log likelihood is likely to be sufficient, but we also develop the method for other parameters, such as forward model parameters.
In the past work on linear problem [5], the evaluation of the curvature matrix has proven to be the most expensive part of the problem. In this paper we develop a novel method to evaluate an approximation to its ensemble average, the Fisher matrix. We use data simulation(s) and investigate their response to the change in parameters, taking advantage of the fact that in Newton’s method the response matrix is also the curvature matrix. This enables us to evaluate this term much faster than otherwise possible. In the linear case this approach remains Bayesian, since the curvature matrix equals the Fisher matrix. In the nonlinear case this is no longer the case, but we will argue that the difference are likely to be small due to the large number of modes.
It is worth emphasizing that if the initial equations 2.1 and 2.2 are exact in a probabilistic sense, and if all the steps we described above are performed exactly, then the final solution is optimal, i.e. we obtain true posterior distribution of the parameters given the data. In this paper we will discuss in detail which approximations we are invoking to solve these equations and what their impact may be on the optimality of the final result.
2.1 Heuristic derivation
The simplest power spectrum estimate one can make out of MAP is to square the modes and average them within the bandpower bin. This does not lead to an unbiased estimator, but intuitively this must be the path to the correct procedure, since it is built from the minimum variance estimator of the initial modes, which in some sense is the best we can do. This procedure has indeed been formally proven in the linear case [7]. In the next sections we present a nonlinear version of this statement. Before proceeding we give a heuristic derivation of our procedure following the arguments in [7].
In general an optimal two point function (e.g. power spectrum) analysis for the linear case requires first inverse generalized noise weighting of the data, after which one adds up all pair products of these inverse noise processed data, weighted by the expected response of each pair product to the specific power spectrum parameter (e.g. bandpower). Inverse weighting by the generalized noise means weighting by the inverse of the full covariance matrix, which consists of noise variance and signal (sampling) variance, and is in general not diagonal in any basis. In the nonlinear case the generalized noise matrix cannot even be defined a priori since it depends on the solution. However, finding MAP of initial modes \hat{\hbox{\boldmath{s}}} achieves this inverse generalized noise weighting. This may still not be sufficient for an optimal analysis, but going beyond it would require doing a full likelihood analysis in dimensional space, which we would like to avoid if possible. So our strategy will be to use this solution to construct the linear power spectrum.
While minimum variance mode reconstruction \hat{\hbox{\boldmath{s}}} has been inverse generalized noise weighted, its square cannot give an unbiased estimator of the power spectrum, and needs to be corrected. First of all, \hat{\hbox{\boldmath{s}}} contains noise contribution, and upon squaring it this will give a term that needs to be removed. Second, minimum variance modes \hat{\hbox{\boldmath{s}}} only agree with truth in the absence of noise and for a complete coverage, but are otherwise reduced in amplitude and mixed with each other. So one needs a mixing matrix to correct for these effects.
Heuristically, we can write an estimator for power at the single mode as
[TABLE]
where is the mixing matrix between the modes and is the noise bias term. This term can be determined by requiring that the equation above is unbiased at some fiducial power spectrum . We have divided by the sampling variance of the power spectrum for a single mode (to be evaluated at some fiducial power spectrum). In the limit where there is no noise the latter is the variance of that mode power (note that we are imagining one independent mode per each , while in practice these are complex modes with two independent mode realizations, but where the mode at is a complex conjugate of the mode at -\hbox{\boldmath{k}}). So in this limit this gives the inverse variance weighting of the power spectrum, where the variance is given by the signal (sampling variance). Noise variance is accounted for by the reduced value of relative to no noise case, an outcome of the minimum variance reconstruction of MAP \hat{\hbox{\boldmath{s}}}. As shown below, the above expression correctly weights by the inverse generalized noise of the power.
We will interpret as the matrix that determines the mixing between the different power spectrum estimators (where in the expression above each mode leads to an independent power spectrum estimator). We will determine both and using simulations around some fiducial model, enforcing that equation 2.4 gives an on average unbiased estimation for any small power spectrum variations around some fiducial model. We will show that we can interpret as the Fisher matrix, i.e. the ensemble averaged inverse covariance matrix for the estimators, as a consequence of Newton’s method. Since equation 2.4 leads to an estimator for each mode it is customary to average them into bandpowers, by adding up estimators within a certain range of wavevectors.
The plan of the next sections is as follows. In section 3 we construct the minimum variance solution for . We will assume and are fixed and we will not carry their dependence. In section 4 we marginalize over and construct maximum likelihood solution for , assuming is fixed. In section 5 we also allow variation of , constructing joint maximum likelihood solution of both and , and providing a general solution. We then extend this to optimal higher order estimators such as bispectrum. In section 6 we apply the method to a simple but nontrivial nonlinear example, which is followed by conclusions in section 7.
3 Maximum a posteriori (MAP) initial mode estimator
In this section we work out the minimum variance modes at fixed around some fiducial model, so we will not include their dependence in the expressions below. For simplicity we will use for the fiducial model. The conditional posterior probability for given the data and prior P_{s}(\hbox{\boldmath{s}}) is P(\hbox{\boldmath{s}}|\hbox{\boldmath{d}})\propto P_{s}(\hbox{\boldmath{s}})P_{n}(\hbox{\boldmath{d-f(s)}}), which from equation 2.1 is
[TABLE]
The MAP solution of the initial modes can be obtained by maximizing the posterior of at a fixed and , which is equivalent to minimizing the loss function ,
[TABLE]
with respect to all coefficients .
The problem of finding the minimum variance field can thus be recast into solving a minimization problem of a nonlinear function , which is an optimization problem. A generic class of optimization models is based on second order Newton’s method, which expands the nonlinear function to a quadratic order around a point \hbox{\boldmath{s}}_{m},
[TABLE]
where \Delta\hbox{\boldmath{s}}=\hbox{\boldmath{s}}-\hbox{\boldmath{s}}_{m}. We can introduce the gradient matrix around as
[TABLE]
Response function depends on the position , but we will not keep that index explicitly in our expressions. For linear models the matrix is a constant independent of position. The gradient of the cost function is
[TABLE]
We want to find a solution where \hbox{\boldmath{g}}=0, and we will denote the solution as \hat{\hbox{\boldmath{s}}},
[TABLE]
Note that is the derivative of the forward model at every position with respect to every initial mode. This can be very expensive to compute for complicated nonlinear models since it requires back-propagation, and is typically not part of the standard N-body simulation codes.
The curvature matrix equals one half of Hessian matrix and can be written as
[TABLE]
The last term above contains a second derivative of the , \hbox{\boldmath{f^{\prime\prime}}}=\partial^{2}{f}/\partial\hbox{\boldmath{s}}\partial\hbox{\boldmath{s}}, and is often neglected when performing optimization, because it fluctuates around zero uncorrelated with the model and hence tends to average to zero. It is also often small (or exactly zero for linear models), and its inclusion may make the optimization unstable [19]. We will also not keep it in our expressions below, and the procedure becomes the so called Gauss-Newton method. Keeping or dropping this term has no impact on the final solution of the optimization.
Minimizing equals setting its gradient with respect to \Delta\hbox{\boldmath{s}} to zero,
[TABLE]
which from equation 3.3 gives
[TABLE]
If we had a perfect inverse curvature matrix and the system was linear this would find the solution where \hbox{\boldmath{g}}=0. In nonlinear Newton’s method one takes this solution as a starting point, performing a line search in the direction of to find a suitable new position \hbox{\boldmath{s}}_{m+1}, which reduces the cost function along the line,
[TABLE]
We start with . If this reduces the cost function we accept it as the next iteration, otherwise we use bisection method to determine a new , repeating the procedure (if no solution is found we switch to regular steepest descent method). This line search does not require a new evaluation of the derivatives, but it does require a full forward model evaluation of for each , and hence can be expensive. The procedure is repeated until one converges to the minimum .
As stated above in linear models one can perform this step only once, starting from \hbox{\boldmath{s}}_{0}=0, to give
[TABLE]
This is the Wiener filter solution, which minimizes the variance [4]. Note that even though the solution is analytic, it requires inversion of a large matrix , which is an process and becomes prohibitively expensive for large . Optimization methods discussed here can be the best choice even for linear models.
At the position of the minimum one can write the conditional probability of given the data as
[TABLE]
where \hat{\hbox{\boldmath{s}}} is the value of found at the minimum. The covariance matrix of residuals is given by the inverse of curvature matrix
[TABLE]
where the gradient is computed at the MAP \hat{\hbox{\boldmath{s}}} and we dropped \hbox{\boldmath{f^{\prime\prime}}}\hbox{\boldmath{N}}^{-1}(\hbox{\boldmath{d}}-\hbox{\boldmath{f}}) term, which fluctuates around zero.
Function optimization (minimization/maximization) is a large field of research, and many different methods have been developed for a wide range of applications [20]. Here we will work in a large number of dimensions (with both and expected to be in millions or more), and a brute force Hessian inversion in equation 3.9 is not feasible. Two of the popular classes of optimization solvers are conjugate gradient methods (often coupled to a preconditioner), which do not require a construction of curvature matrix (or its inverse) at all, and quasi-Newton’s methods, which approximate the inverse Hessian, of which limited memory BFGS (L-BFGS) is the most popular implementation [19]. We have used Gauss-Newton’s approach in the derivation above to highlight the analogy with linear algebra derivation of [7], but ultimately the choice of optimization method will be determined by its effectiveness. In our tests we have found L-BFGS to be faster than conjugate gradient, but this could depend on the application and implementation. The derivation above does not include normalization, but in practical applications we have found to be useful to rescale (normalize) and work with the modes \hbox{\boldmath{u}}=\hbox{\boldmath{s/S}}^{1/2}.
4 Minimum variance initial power spectrum estimator
While knowing the MAP solution \hat{\hbox{\boldmath{s}}} is useful to create a map of the universe, this is still too much information for most applications. What we want to know is the statistical distribution of the linear modes, and since these are assumed to be gaussian a suitable summary statistic is their power spectrum. Basically, we want to know the parameters , which are priors for and marginalize over . In this section we address the dependence on parameters , while still keeping parameters fixed at their fiducial values (and we will not carry their dependence).
Formally we want to maximize the posterior probability of the data as a function of some power spectrum parameters that determine the initial mode power spectrum , independent of the underlying field \hat{\hbox{\boldmath{s}}}. We can introduce a derivative matrix \hbox{\boldmath{\Pi}}_{l} around some fiducial power spectrum \hbox{\boldmath{S}}^{\rm fid}, defined as
[TABLE]
which is to be evaluated at the fiducial model. In terms of this we can write
[TABLE]
For linear dependence of on and bandpowers we can use
[TABLE]
i.e. . In this case the derivative matrix takes us from , which is the power spectrum value representative over a bin, to that is the power spectrum. For concreteness we will initially bin the power spectrum into bins, summing over all the modes . These are the modes that contribute to bandpower of the power spectrum parametrized with . For a narrow bin the derivative matrix \hbox{\boldmath{\Pi}}_{l} consists of ones along the diagonal corresponding to the modes and zeros otherwise. For broader bandpower bins the derivative matrix can allow for any expected variation of the power within the bin. For example, if we only have a single bin which covers the entire range of modes then the corresponding is a power spectrum amplitude chosen at some wavemode (or some integral over the modes, such as normalization), while the derivative matrix has the shape of a fiducial power spectrum, i.e. \hbox{\boldmath{\Pi}}_{1}=\hbox{\boldmath{S}}_{\rm fid}/\Theta_{1}.
4.1 Analytic marginalization with perturbative approach
In the case of power spectrum estimation one wants to marginalize the full posterior over the parameters \hat{\hbox{\boldmath{s}}}. Since a brute force approach with sampling is too expensive, we choose to approximate it as a multi-variate gaussian as in equation 3.12. An implicit assumption in the procedure is that there is one global minimum, to which optimization converges regardless of the starting point. This is a major assumption and may not be valid, either because the global minimum is not reachable (because optimization is stuck in a local minimum), because the correct solution does not correspond to the global minimum, or because the posterior around the global minimum is not well described as a multi-variate gaussian. Here we will assume this approach is justified but there are no guarantees and indeed on very small scales the methodology may need to be modified.
Under this assumption our strategy is to marginalize over the modes, by writing it as a multi-variate gaussian around the maximum, and then complete the square. Note that the MAP \hat{\hbox{\boldmath{s}}} and its covariance is implicitly a function of the data, the forward model parameters and the power spectrum bandpower parameters we wish to estimate. In [21] the resulting expression was called a grand likelihood L(\hbox{\boldmath{d}}|\hbox{\boldmath{\Theta}}). The full problem can be solved if we linearize the response around the MAP initial mode estimator \hat{\hbox{\boldmath{s}}} and perform analytic integration under gaussian approximation. We begin with equation 3.1, expanding around \hbox{\boldmath{s}}=\hat{\hbox{\boldmath{s}}}+\delta\hbox{\boldmath{s}},
[TABLE]
If we collect all the terms that involve \delta\hbox{\boldmath{s}} they will vanish by the requirement that the gradient is zero, equation 3.6, by definition of \hat{\hbox{\boldmath{s}}}: the posterior surface is quadratic at the minimum of the loss function. We are left with
[TABLE]
where is the curvature matrix of equation 3.7 and more generally we have defined
[TABLE]
We see that \hbox{\boldmath{D}}=2\hbox{\boldmath{D}}_{2}. In terms of derivatives of we have
[TABLE]
and
[TABLE]
We have consistently dropped the terms with \hbox{\boldmath{d}}-\hbox{\boldmath{f}}(\hat{\hbox{\boldmath{s}}}), which fluctuate around zero. All the terms in the last line of equation 4.5 are scalars.
Next we assume a perturbative expansion in terms of higher order derivatives and/or low noise (we will define the expansion parameter below), bringing the higher order terms \hbox{\boldmath{D}}_{n} down from the exponential,
[TABLE]
We can now obtain the perturbative expansion of the marginalized likelihood function by integrating out \delta\hbox{\boldmath{s}},
[TABLE]
We have defined
[TABLE]
The lowest order term \det\hbox{\boldmath{D}} is given by the gaussian integral. The next order term is a 1-loop contribution consisting of and inside \hbox{\boldmath{D}}_{4}, and we used Wick’s theorem to perform the gaussian integrals. We can also define
[TABLE]
where denotes higher order terms in \tilde{\hbox{\boldmath{D}}} and we used equation 3.13 for .
We see that to marginalize over a set of parameters one has first to find their maximum posterior \hat{\hbox{\boldmath{s}}}, as derived in previous section via optimization, after which the integration over the parameters can be performed analytically assuming a perturbative expansion. This procedure implicitly assumes a single maximum posterior and a simple posterior around \hat{\hbox{\boldmath{s}}}. When this is not satisfied, for example in the presence of multiple peaked posterior, one must either perform the appropriate statistical average over these, if they are widely separated, or show that all the local peaks are close to each other and thus essentially giving the same solution. For now we will proceed assuming a single peaked posterior. Note that integration over produced an additional factor of \tilde{\hbox{\boldmath{D}}}^{-1/2} relative to doing the maximum likelihood over all parameters simulteneously.
Let us look at the structure of the resulting expression in equation 4.10. The main term that depends on the data is in the exponential, \exp[-\chi^{2}(\hat{\hbox{\boldmath{s}}})/2]. This term consists of the prior term and of the goodness of data fit, both evaluated at the MAP position \hat{\hbox{\boldmath{s}}}. The other term is \det\hbox{\boldmath{SN}}\tilde{\hbox{\boldmath{D}}}. This term represents the volume of the priors and posteriors. We wish to explore its sensitivity to the parameters that determine , so we wish to take derivative of this term wrt .
From equation 4.12 we obtain
[TABLE]
In the low noise limit it gives \det(\tilde{\hbox{\boldmath{C}}})\approx\det(\hbox{\boldmath{f^{\prime}}}\hbox{\boldmath{S}}\hbox{\boldmath{f^{\prime}}}^{{\dagger}}), which is highly sensitive to the parameters inside . In the high noise limit we have \det(\hbox{\boldmath{SND}})\approx\det\hbox{\boldmath{N}}. In this limit we do not get any information about the parameters . This is reflected in the values of \hat{\hbox{\boldmath{s}}}, which are at the their full value in the low noise limit and become smaller and smaller as the noise increases.
Next we look at the data dependence of the curvature matrix \tilde{\hbox{\boldmath{D}}}. If the model is linear \tilde{\hbox{\boldmath{D}}} does not depend on the data, because is a constant and for . Thus this term can be evaluated using a method that is data independent, such as a simulation. For nonlinear models we have to account for the data dependence of \tilde{\hbox{\boldmath{D}}}. Let us first look at the leading term , given in equation 3.13. For nonlinear models the first derivative \hbox{\boldmath{f^{\prime}}}(\hat{\hbox{\boldmath{s}}}) is data dependent. For example, suppose we Taylor expand the forward model in terms of ,
[TABLE]
we see that
[TABLE]
and hence the covariance depends on \hat{\hbox{\boldmath{s}}} through \hbox{\boldmath{R}}_{2}\hat{\hbox{\boldmath{s}}}+.... Depending on the realization value of \hat{\hbox{\boldmath{s}}} the associated error can be larger or smaller than the first order term \hbox{\boldmath{R}}_{1}. We can see this from equation 3.6. The solution to this equation balances between the prior term \hbox{\boldmath{S}}^{-1}\hat{\hbox{\boldmath{s}}}, which on its own would give \hat{\hbox{\boldmath{s}}}=0, and the data term \hbox{\boldmath{f^{\prime{\dagger}}(\hat{s}){N}^{-1}[d-f(\hat{s})}}], which on its own would solve for \hbox{\boldmath{f}}(\hat{\hbox{\boldmath{s}}})=\hbox{\boldmath{d}}. The balance between the two is determined by , which contains data realization dependence inside \hbox{\boldmath{f^{\prime}}}(\hat{\hbox{\boldmath{s}}}). Thus the solution for \hat{\hbox{\boldmath{s}}} automatically accounts for generalized noise weighting of the modes including realization dependence, at least at the lowest order.
Let us look next at the 1-loop term {\rm tr}(\hbox{\boldmath{D}}^{-1}\hbox{\boldmath{D}}_{4}\hbox{\boldmath{D}}^{-1}). This term represents an increase in volume of posterior of \delta\hbox{\boldmath{s}} around \hat{\hbox{\boldmath{s}}} relative to the gaussian approximation. Using equations 3.13 and 4.8 we find this term scales as, schematically,
[TABLE]
Thus the 1-loop term vanishes in the limit of either small noise or small nonlinearity. Both of these conditions are likely to be satisfied on large scales, where noise relative to power is small and modes are nearly linear. As we push the marginalization to smaller scales we get larger and larger corrections from this term. This term is also data dependent. This means that one may have a situation where Taylor expansion suggests the posterior is wide (narrow), because the variance of the mode \hbox{\boldmath{D}}^{1/2} is large (small), but the full posterior may be narrower (wider) once we include higher order terms, such as curtosis computed here. At this order our procedure will cease to be optimal in terms of optimal mode weighting, but we will still make it unbiased, as described below. It is unclear however how suboptimal due to higher order corrections our method is. We will show below that our method is essentially giving equal weight to all reconstructed modes, and this may well be valid even beyond the formal applicability regime derived here.
Optimization methods in large number of dimensions do not provide \det\hbox{\boldmath{D}}, although approximations to it exist in certain methods (e.g. BFGS). It seems even less promising to be able to compute analytically \det\tilde{\hbox{\boldmath{D}}}. While this term is data realization dependent, it is computed independently of the \exp[-\chi^{2}(\hat{\hbox{\boldmath{s}}})] term, and is essentially determined by the posterior volume of \delta\hbox{\boldmath{s}} around \hat{\hbox{\boldmath{s}}}. In a large number of dimensions the average of this term will not be strongly data realization dependent: for any given realization we will have about the same level of mode to mode fluctuations that will on average have the same posterior volume and give the same value of \det\tilde{\hbox{\boldmath{D}}} regardless of where \hat{\hbox{\boldmath{s}}} is. This suggests we may evaluate this term average using Monte Carlo methods, with a realistic simulation of the data: even though the MAP values \hat{\hbox{\boldmath{s}}} in a simulation will have no relation to the corresponding MAP of the data, the posterior volume of \delta\hbox{\boldmath{s}} around it will be nearly the same as for the data, as long as the simulation is close enough to the data. We may give up some optimality by not properly weighting the modes by their full posterior volume, but our procedure is exact in the linear case and in the low noise case, and as argued above equal weighting of the modes may be the valid approach even beyond these two formal limits. In any case, we cannot evaluate the posterior volume properly using the methods developed here, which are perturbative anyways and so not valid in general. This will be the basis for our efficient method of evaluating the optimal estimator and its covariance matrix. Note that even though we use random simulations in this paper to compute the volume term, we may also construct a simulation that is close to the actual data realization, further reducing the dependence of the volume term on the data realization.
4.2 Maximizing the bandpower likelihood
Using the Bayes theorem and assuming a flat prior on we have interpreted the posterior P(\hbox{\boldmath{\Theta}}|\hbox{\boldmath{d}}) as proportional to the likelihood L(\hbox{\boldmath{d}}|\hbox{\boldmath{\Theta}}). We can justify using flat prior as the choice which gives unbiased estimators of cosmological parameters, when the summary statistics are used to determine them. This interpretation is only valid when the noise bias and Fisher matrix are evaluated using a fiducial power spectrum not affected by the measurements. When the fiducial model is too far from the final solution one can repeat the procedure with a fiducial model that is closer to the actual solution. Note that we do not impose positivity of the bandpowers with the prior: this is justified if the bandpowers are viewed as summary statistics that are subsequently used to determine cosmological parameters (where physical priors can be imposed).
The next step is to maximize this likelihood with respect to the bandpowers . Given a set of measurements we wish to find the most probable set of bandpowers from equation 4.10, where , and \hat{\hbox{\boldmath{s}}} implicitly depend on the bandpowers . To find the most probable set of parameters one needs to find the maximum of the likelihood function L(\hbox{\boldmath{\Theta}}).
We will find the maximum posterior using Newton’s method. We wish to expand the log likelihood in terms of to a quadratic order around some fiducial values,
[TABLE]
where the terms are evaluated at the fiducial model \hbox{\boldmath{\Theta}}_{\rm fid}. We will seek a solution where the first derivative vanishes using Newton’s method. The last term of equation 4.17 defines the curvature matrix as the second derivatives of log likelihood with respect to the parameters.
To get the first derivative we need to take a derivative of log likelihood in equation 4.10 with respect to . The dependence on is in \hat{\hbox{\boldmath{s}}}, and \hbox{\boldmath{f}}(\hat{\hbox{\boldmath{s}}}). We have
[TABLE]
where the first relation follows from equation 3.4 and the last relation follows from equation 3.6. We also have
[TABLE]
defining
[TABLE]
All the matrix operations involve diagonal matrices, so using equation 4.3 we get
[TABLE]
where the sum is over all modes within the bandpower and is the fiducial power spectrum amplitude at mode , which for narrow power spectrum bins is simply . We see that the terms with d\hat{\hbox{\boldmath{s}}}/d\Theta_{l} cancel out, which gives us, using equation 4.10
[TABLE]
where we defined
[TABLE]
Note that has no dependence on . This term is non-zero because squaring the MAP field (as done in first term of rhs in equation 4.22) creates a bias. In line with the linear analysis we will continue to call this term noise bias, but unlike the linear case this term may depend on both the noise and the signal. The first term in equation 4.22 states that for narrow bins the derivative is given by adding up all the squares of the modes. We can see that in this term the curvature matrix (or its generalization \tilde{\hbox{\boldmath{D}}}) does not enter, and its evaluation is not required. Since we cannot compute the curvature matrix we also cannot take its derivative with respect to parameters, which is needed to compute the noise bias in equation 4.23. In this paper we instead compute the noise bias by evaluating its ensemble average, using a method described further below. This will mean that we have given up on the information contained in this term. We expect this information to be subdominant, at least in the linear and low noise limits.
The maximum likelihood solution for \hat{\hbox{\boldmath{\Theta}}} is given by
[TABLE]
To evaluate we would need to evaluate equation 4.23. This gives
[TABLE]
Using equation 4.12 and dropping higher order terms we obtain
[TABLE]
We wish to solve equation 4.24 to obtain ML estimate of , which appears inside . This is a complicated nonlinear equation in terms of and hence one cannot write the solution in a closed form.
Instead of solving this equation directly we will adopt Newton’s method, where we use the quadratic expansion of log-likelihood in equation 4.17, evaluate all the elements at \hbox{\boldmath{\Theta}}_{\rm fid} where the log-likelihood derivative is not zero, and then use Newton’s method to find the values \hat{\hbox{\boldmath{\Theta}}} where the first derivative is zero. To do this we must therefore evaluate the curvature matrix. Instead of the actual curvature matrix for parameters we will often work with the Fisher matrix, defined as the ensemble average of the curvature around the maximum,
[TABLE]
Brackets denote ensemble averaging. It may be possible to compute the actual curvature matrix for the data (discussed further below and in appendix B), but in this paper we will not distinguish between the two and simply call it Fisher matrix, even though this is not exactly the same as the curvature matrix in the Bayesian analysis. For modes that are gaussian distributed one can interpret the inverse of the Fisher matrix \hbox{\boldmath{F}}^{-1} as an estimate of the covariance matrix of the parameters \hat{\hbox{\boldmath{\Theta}}},
[TABLE]
Note that for linear modes the likelihood function can be approximated as gaussian around the maximum provided that sufficient number of independent modes contribute to each , by central limit theorem. Appendix A discusses an approximation that goes beyond this limit.
We are finally in position to use Newton’s method to write the maximum likelihood estimator for \Delta\hbox{\boldmath{\Theta}}. The solution to the quadratic log likelihood at the peak can be obtained by setting the derivative of equation 4.17 with respect to to zero, which upon inserting equation 4.22 gives
[TABLE]
where is the noise bias contribution. To solve this equation one needs both and . In principle we could get the latter by evaluating equation 4.27. This would give an expression in terms of a matrix multiplications and trace,
[TABLE]
which is formally of order , if we had the matrices. But we do not actually have the matrix \tilde{\hbox{\boldmath{C}}} or \tilde{\hbox{\boldmath{D}}}, so we cannot evaluate this even if we wanted to. We will thus evaluate the Fisher matrix using a different approach using ensemble averaging of simulations. However, in the special case of no mode coupling the Fisher matrix is diagonal and we can combine equations 4.24, 4.25 and 4.30 to obtain
[TABLE]
In this case we get the actual curvature matrix from the data itself. We will test the accuracy of this approximation in section 6 for the specific case of periodic boundary conditions, which is the only case where the assumption of uncorrelated modes may be valid.
From equation 4.29 we see that the raw power spectrum of the reconstructed field requires noise bias subtraction. We also see that after noise bias subtraction we obtain the bandpower estimates convolved with the Fisher matrix \hbox{\boldmath{F}}\Delta\hat{\hbox{\boldmath{\Theta}}}. Fisher matrix thus describes both the covariance matrix and the bandpower mixing. It is the latter interpretation that will be the basis for our fast Fisher matrix evaluation method. As promised, equation 4.21 agrees with our heuristic derivation in equation 2.4. The procedure is implied to be iterative: if the chosen fiducial model is not sufficiently close to the final answer one should choose a new fiducial model and repeat the procedure. The fiducial model should not be simply chosen to be the observed one given by \hbox{\boldmath{\Theta}}_{\rm fid}+\Delta\hat{\hbox{\boldmath{\Theta}}}, since that contains noise and sampling variance fluctuations. Instead, it should be a smooth model that is sufficiently close to it. Equation 4.29 suggests that all the information from the data is inside E(\hbox{\boldmath{S}}_{\rm fid},\hat{\hbox{\boldmath{s}}}). In reality, and may also be data dependent (although this is explicitly not the case for the linear model). This is discussed further below.
The form given in equation 4.29 is not the only possibility, and may not be the most practical. More generally, one can write
[TABLE]
where is a matrix that can be chosen depending on the form in which we wish to provide the estimator \Delta\hbox{\boldmath{\Theta}}^{,}, which can be viewed as a window convolved version of \Delta\hbox{\boldmath{\Theta}}. There are three natural choices [5] that we discuss next.
First choice is \hbox{\boldmath{M}}=\hbox{\boldmath{F}}^{-1}, which corresponds to the fully deconvolved estimator (\Delta\hbox{\boldmath{\Theta}}^{,}=\Delta\hbox{\boldmath{\Theta}} from equation 4.29), with anti-correlated errorbars (which can explode if the sampling of modes is too fine compared to the width of the window).
A second choice is
[TABLE]
which corresponds to the minimum variance errors on , but which are all convolved with a window given by . We will use this choice in this paper because it does not require a matrix inversion. Division by ensures the estimator is centered around the true value. However, this can cause problems if the matrix is singular, which can happen if there is no information about a given bandpower in the data. In this case it is better not to divide by this factor and simply use .
A third choice is to decorrelate the estimators (within the disconnected covariance matrix approximation) [22],
[TABLE]
With any of these definitions we can write the expectation value and variance of the estimator as
[TABLE]
where the latter contains only the connected part of covariance matrix and we defined the window matrix \hbox{\boldmath{W}}=\hbox{\boldmath{MF}} and the covariance matrix \hbox{\boldmath{C}}=\hbox{\boldmath{MFM}}. All of these choices are equivalent in terms of their information content: there is no gain or loss of information in choosing one over the other. All the data compression has been done in equation 4.32. The first choice corresponds to identity window matrix, the second is minimum variance and does not require any matrix inversion (equation 4.35) and is thus the easiest to compute, while the third choice provides a diagonal covariance matrix of the estimator, so bandpowers are easy to combine.
The estimator is unbiased as long as the Fisher matrix is computed using a fiducial power spectrum, not affected by the data. If the chosen fiducial power spectrum is very different from the one preferred by the data the procedure is still unbiased, but the estimator may not be optimal and the Fisher matrix may not be accurately interpreted as a covariance matrix. In these situations one may repeat the analysis with a chosen fiducial power spectrum closer to the one preferred by the data. In the case of high precision cosmology where the deviations between the true and fiducial power spectrum are small, the corresponding error induced on the covariance matrix also becomes small.
The estimator presented here (quadratic in MAP) approximates the posterior for bandpowers as a multi-variate gaussian. This allows us to normalize the posterior, and assuming bandpowers we have
[TABLE]
In the linear case this is valid in the limit of large number of modes by central limit theorem. Typically one combines different bandpowers to determine a small number of cosmological paramaters, which makes central limit theorem even more valid. However, in special cases, such as clustering on very large scales with few measured modes, the small number of modes requires one to use a more accurate posterior probability distribution. We discuss this situation in appendix A. The connected part of the covariance matrix is not given by the Fisher matrix construction presented above. However, since we are estimating the linear power spectrum, which is assumed to be gaussian, it is likely that there is no significant connected part, at least for low noise. To address this properly one needs to compare the disconnected covariance matrix to the one from a large number of mocks where the whole procedure has been simulated, and we will not investigate this in any detail in this paper.
4.3 An efficient evaluation of noise bias and Fisher matrix
What is left is to evaluate the noise bias and Fisher matrix. This is often the most expensive part of the calculation, even for the linear model [5]. Here we are solving the nonlinear problem in which the minimization cannot be solved with linear algebra anyways, so we have to devise a new scheme to evaluate these terms.
Before proceeding it is useful to clarify the nature of the whole procedure. We began by estimating the modes that minimize the variance of the data given the noise and the prior. In doing so we have chosen a fiducial power spectrum \hbox{\boldmath{S}}_{\rm fid} for the prior. We then square these modes and add them up within the bandpower . We want to make this estimator unbiased, and with the correct window function describing correlations between the bandpowers. Our primary task is to have the estimator that is unbiased for any choice of the fiducial power spectrum \hbox{\boldmath{\Theta}}_{\rm fid}. In practice, we expect the estimator will be almost as good as minimum variance even if the fiducial power spectrum chosen in the minimization procedure is not exactly the true power spectrum. Fisher matrix plays two distinct roles in the parameter estimation, as a window function (equation 4.35) and as a covariance matrix (equation 4.36). Here we will use its interpretation as the window function to derive it, quantifying the sensitivity of one bandpower to another. This only involves two point functions, and does not rely on the assumption of gaussianity. But Newton’s method guarantees that the resuting answer is also the covariance matrix.
Since a direct evaluation is not computationally feasible, we will be evaluating Fisher matrix and noise bias terms using simulations. The first question is whether these terms depend on the actual realization. For example, the noise bias is given by the derivative of the log determinant of the curvature matrix with respect to the bandpowers (equation 4.23) and this can depend on the actual realization. This is however not the case in the linear limit [7]. It is also not very important in the low noise limit, where most of the information is in the reconstructed modes and the noise bias is small. More generally, we have argued in the beginning of this section that the volume element \det\tilde{\hbox{\boldmath{D}}} of \delta\hbox{\boldmath{s}} around \hat{\hbox{\boldmath{s}}} can be determined from a simulation in the limit of large number of modes even if the values of \hat{\hbox{\boldmath{s}}} differ between a simulation and the data. We will thus proceed by assuming that there is no realization dependence of Fisher matrix and noise bias, and we will determine them using simulations that are unrelated to the actual realization given by the data. Should this assumption turn out to be violated one could explore realization dependence of these terms, for example by using realizations that are conditioned on the data (such as Gibbs sampling), but we do not pursue this further in this paper. We emphasize that our estimator will still be unbiased at the fiducial model and small variations around it, but it may not be minimum variance.
For the simulations we thus first generate a gaussian random realization of the signal in Fourier space \hbox{\boldmath{s}}_{s} at the fiducial model, such that
[TABLE]
We have explored both a gaussian realization, which contains sampling variance fluctuations, and the fixed norm of the mode, where only the phase of the mode is random, while the mode square is exactly given by the fiducial power. The latter does not suffer from sampling variance and converges faster if no sampling variance correction is used, although we will also give prescriptions how to handle sampling variance for the former case. Suppression of sampling variance is crucial and guarantees faster convergence, so the procedure here should not be viewed simply as a Monte Carlo simulation approach.
Next, we map it into the observational space using a forward model,
[TABLE]
We also generate a noise sample. Similar to the signal, we have explored both a gaussian realization with unit variance, and the real stochastic estimator, where for each real space point \hbox{\boldmath{r}}_{i} we generate a random number consisting of 1’s and -1’s. The latter has a fixed amount of power and is expected to converge faster. In both cases we multiply with \sigma(\hbox{\boldmath{r}}_{i}), the square root of the noise variance at that point N(\hbox{\boldmath{r}}_{i})=\sigma^{2}(\hbox{\boldmath{r}}_{i}), to create a random noise data vector \hbox{\boldmath{d}}_{n} with the correct variance. We use this noise realizations \hbox{\boldmath{d}}_{n} added to the signal data to create
[TABLE]
Now that we have a data simulation we pass this through the same sequence of optimization steps as the real data. This leads us to the simulated realizations of the mode estimator \hat{\hbox{\boldmath{s}}}_{s+n} for the fiducial model. The cost for performing a single realization is comparable to the cost of analyzing the data.
We want to ensure that the final estimates are unbiased. From equation 4.29 the noise bias is
[TABLE]
This equation is simply stating that for a fiducial simulation the estimate of bandpower relative to fiducial power has to vanish. If needed, this needs to be averaged over several realizations.
The window function interpretation of the Fisher matrix is a response of a bandpower to another bandpower at , as in equation 4.29. A way to evaluate it is to inject power into one bandpower and measure the estimator response in all the bandpowers (in statistics this is called a sensitivity analysis). Specifically, we take the fiducial power spectrum realization from equation 4.40. We create another gaussian realization as a small perturbation to it, by injecting a small amount of power into the modes at a single bandpower . For each of the modes within the bandpower we have a choice to generate either a random realization of the phase, or one given by . We add it to fiducial realization \hbox{\boldmath{s}}_{s}
[TABLE]
We want the injected power to be \Delta\Theta_{l^{\prime}}\hbox{\boldmath{\Pi}}_{l^{\prime}} (where should be small so that linearization of the Fisher matrix is valid, but not too small to be susceptible to roundoff errors), so we choose \Delta\hbox{\boldmath{s}}_{l^{\prime}} to satisfy
[TABLE]
We have found that the injected power with modes in phase with \hbox{\boldmath{s}}_{s} works better in terms of the noise in the Fisher matrix.
We then run a forward model to generate \hbox{\boldmath{d}}_{l^{\prime},s}=f(\hbox{\boldmath{s}}_{l^{\prime},s}), to which we add exactly the same noise realization as the one in equation 4.40, \hbox{\boldmath{d}}_{l^{\prime},s+n}=\hbox{\boldmath{d}}_{l^{\prime},s}+\hbox{\boldmath{d}}_{n}. We pass it through the same sequence of optimization steps as the fiducial model simulation. Since we already have the optimization solution for the fiducial model simulation at \hbox{\boldmath{d}}_{s}+\hbox{\boldmath{d}}_{n}, and the injected modes \Delta\hbox{\boldmath{s}}_{l^{\prime}} are a small perturbation in a single bandpower, we can start the optimization at . This requires very few steps to converge, giving \hat{\hbox{\boldmath{s}}}_{l^{\prime},s+n}. We compute for each bandpower
[TABLE]
We can apply equation 4.29 to and , but still performing the optimization around the fiducial model,
[TABLE]
We are assuming that the Fisher matrix changes are at higher order, so that is the same at the fiducial model and small perturbation away from it. In general, the noise bias depends on the signal and we cannot assume b_{l}(\hbox{\boldmath{\Theta}}_{\rm fid},\hat{\hbox{\boldmath{s}}}_{s+n})=b_{l}(\hbox{\boldmath{\Theta}}_{\rm fid},\hat{\hbox{\boldmath{s}}}_{l^{\prime},s+n}). However, in the linear regime this condition is satisfied, and even when it is not we still obtain an estimator that is unbiased, but whose covariance may not be given by the resulting .
Taking the difference between the two terms in equation 4.45 we obtain
[TABLE]
where the last equality follows if the additional injected modes are uncorrelated with fiducial modes. We have defined a new matrix \tilde{\hbox{\boldmath{F}}} as the sum of the true Fisher matrix and the matrix of derivatives of noise bias, which will only approximate the true Fisher matrix (although the relation is exact in linear case where the noise bias derivative term vanishes). In the remainder of this paper we will ignore this distinction between and \tilde{\hbox{\boldmath{F}}}, although we will verify the quality of resulting covariance matrix in our tests. The Fisher matrix is then obtained from equation 4.46, averaged over several realizations if needed. It is important to recognize that our final expressions do not explicitly evaluate the curvature matrix , hence any optimization method can be used to find MAP. An alternative method using the auto and cross power of reconstructed modes is presented in appendix B.
The Fisher matrix is rapidly changing with the bandpower , mostly because typically both the number of modes and the fiducial power are rapidly changing with . However, we can also define normalized version
[TABLE]
which can be viewed as an average squared transfer function as in equation B.2. This has the advantage of being a slowly changing function of . As a consequence, one does not need to evaluate it at every , especially for narrow bandpower bins. This allows one to determine the Fisher matrix by injecting the power into sparsely separated bandpowers, and then interpolate between averaged squared transfer function to get the full Fisher matrix.
Sometimes we can adopt a simple approximation of equation 4.31 (see appendix B and [21]), which follows from the covariance interpretation of the Fisher matrix, together with the gaussian variance approximation. We will show below this works remarkably well in simple periodic box test cases.
5 Beyond the bandpowers
5.1 Modeling the forward model parameters
In the previous section we have assumed that the mapping between the initial modes, whose statistics is described with the initial power spectrum , and the data is unique, in which case all the information is in the optimal power spectrum estimators of the bandpowers of . In reality this is almost never the case. First of all, the nonlinear mapping \hbox{\boldmath{f}}(\hbox{\boldmath{s}}) depends on other parameters that determine the growth of structure, such as matter density, dark energy equation of state etc. Observations are also not directly in terms of matter over-density in comoving coordinates: for example, weak lensing observations can be expressed as an integral over the over-density times the mean matter density, and are typically observed as a function of angle, which can be converted to physical scale assuming cosmological model parameters such as dark matter, dark energy density etc. In addition, there are uncertainties in the physical modeling, such as baryonic effects, which depend on unknown astrophysical parameters and can be parametrized in terms of these. There are also data modeling uncertainties, for example, in the case of weak lensing both multiplicative and additive shear calibration biases and intrinsic alignments [23].
When it comes to galaxies the situation is even more complicated, since there are many ways to assign galaxies to dark matter halos, and one needs to allow for all models that are consistent with our present understanding of galaxy formation physics, yet currently cannot be determined from an ab initio modeling. For example, the simplest way is to parametrize galaxy clustering in terms of a free galaxy bias, which cannot be assumed to be known, but needs to be marginalized over. In addition, there are forward model parameters that control redshift space distortions (growth of structure) and mapping from comoving distances to angular and radial distance.
If the nuisance parameters affect only a few modes one can simply include these during the optimization process. For example, one may know the spatial templates that describe the systematics such as star contamination or dust extinction effects. The procedure to handle this situation is described in appendix C. In the absence of a prior the procedure is equivalent to setting the noise in these modes to a very large value and optimization sends the modes to zero. This eliminates cosmological information contained in those modes. For a few modes only this may not lead to a major loss of information. However, most of the nuisance parameters described above, such as galaxy bias, affect all the modes, and marginalizing over them during the optimization procedure would send most of the modes to zero. Instead, these parameters must be handled at the level of likelihood function analysis, similar to the power spectum bandpowers.
Our general strategy how to handle these parameters will thus be the same as the bandpowers: we write the likelihood function in terms of these parameters and the bandpowers, expand it to second order, and use Newton’s method to write down the solution at a quadratic order. Since these parameters enter in the nonlinear mapping f(\hbox{\boldmath{s}}) one must perform the corresponding partial derivatives. A very useful property of bandpowers is that they are linear in terms of covariance matrix , and as a result a single step Newton’s method usually suffices, and their posterior is well described as a multivariate gaussian (or inverse-Wishart if the bandpower consists of a few modes only). This can no longer be assumed for the parameters that affect the forward model, and one may need to iterate several times before reaching the maximum of the likelihood. The posterior may also be more complicated than a gaussian (gaussian posterior is the reason that bandpowers are preferred over the cosmological parameters as the expansion basis for , even though the number of parameters for the latter may be a lot lower). Here we will expand to quadratic order and leave the question of more accurate description of resulting posterior to future work.
To proceed we will linearize forward model response around MAP solution \hat{\hbox{\boldmath{s}}}. The likelihood is given by equation 4.10. We want to expand it in terms of the forward model parameter . The terms containing derivative of \hat{\hbox{\boldmath{s}}} cancel out as in equation 4.22. Let us look at the remaining data dependent term of the likelihood that depends on , which gives
[TABLE]
where
[TABLE]
and the derivative of with respect to is to be evaluated at fiducial \hbox{\boldmath{\lambda}}_{\rm fid} and fixed . This term is straight-forward to evaluate with the forward model: one computes forward model around the fiducial value (for example at ) and take the difference to the fiducial model to obtain
[TABLE]
Note that this must be evaluated at the actual realization consistent with the data. In the presence of noise this could be problematic: one possibility is to average over realizations consistent with the data, or use analytic expressions. An example is given below in the context of primordial non-gaussianity. Assuming the noise matrix is diagonal the sum is over all the data pixels ,
[TABLE]
This expression looks deceptively simple, but note that enumerator and denominator both vanish in the zero noise limit, so evaluating this expression can be numerically challenging in the low noise regime. We will convert it into mode space below. We also defined
[TABLE]
since matrix \tilde{\hbox{\boldmath{D}}} is the only one that depends on the forward model parameters.
From equation 5.2 it would appear that the forward model parameters are limited by noise only, and not sampling variance. However, as already mentioned above the noise term actually cancels out of this expression in the low noise limit. One can rewrite equation 5.3 using equation 3.6 into
[TABLE]
It is simplest to analyze this in the first order Taylor expansion case where \hbox{\boldmath{f}}(\hat{\hbox{\boldmath{s}}})=\hbox{\boldmath{f^{\prime}}}\hat{\hbox{\boldmath{s}}}, in which case we obtain
[TABLE]
Comparing to expression for in equation 4.20 we see the two are the same if
[TABLE]
This is what one would expect if, for example, the forward model parameter simply rescales the amplitude, in which case it is degenerate with the overall linear amplitude parameter \hbox{\boldmath{S}}\propto\Theta_{1}^{2}, for which \hbox{\boldmath{\Pi}}_{l}=2\hbox{\boldmath{S}}/\Theta_{1}. Another way to obtain the same result is by rewriting the data likelihood from equation 4.10 in the linearized case, using equations 3.11 and 3.7 (dropping \hbox{\boldmath{f}}" for linearized models) into
[TABLE]
where is defined in equation 4.12. It is clear from this expression that the dependence of the likelihood with respect to the forward model parameters in and with respect to the power spectrum parameters in appears in the same term \hbox{\boldmath{f^{\prime}}}\hbox{\boldmath{S}}\hbox{\boldmath{f^{\prime}}}^{\dagger} contained inside \tilde{\hbox{\boldmath{C}}}.
We can rewrite equation 5.7 and make it applicable to the nonlinear model,
[TABLE]
Here again the derivative has to be done on the actual reconstructed modes \hat{\hbox{\boldmath{s}}}. The interpretation of this derivative is as follows: we forward model the solution \hat{\hbox{\boldmath{s}}} with a small change in forward model , \hbox{\boldmath{f}}(\hat{\hbox{\boldmath{s}}},\lambda_{l,{\rm fid}}+d\lambda_{l}), and then perform the optimization steps around the fiducial model to get the new solution \hat{\hbox{\boldmath{s}}}_{\lambda_{l,{\rm fid}}+d\lambda_{l}}. Taking the difference with respect to \hat{\hbox{\boldmath{s}}} and dividing by gives us \partial\hat{\hbox{\boldmath{s}}}/\partial\lambda_{l}. This can be numerically noisy, and analytic expressions for the derivative may give better performance. In the linear regime these are easy to derive, and perturbative methods can be used to derive them beyond the linear theory. We will discuss an application to primordial non-gaussianity below.
The remaining analysis follows the steps outlined in the previous section. We can combine the union of all parameters into a vector \hbox{\boldmath{p}}=(\hbox{\boldmath{\Theta}},\hbox{\boldmath{\lambda}}). Defining the fiducial values of all the parameters as \hbox{\boldmath{p}}_{\rm fid} we can work with parameters relative to fiducial values \Delta\hbox{\boldmath{p}}=\hbox{\boldmath{p}}-\hbox{\boldmath{p}}_{\rm fid}. The estimator for \Delta\hbox{\boldmath{p}} parameters is then
[TABLE]
Noise bias is given by
[TABLE]
i.e. it is the estimator evaluated at the fiducial model simulation, guaranteeing the estimator is unbiased. The Fisher matrix can be computed as in the previous section by changing a single model parameter in the forward model and observing the response to all the parameters , both the bandpowers and the model parameters.
5.2 The general solution
In general, a given parameter can affect both the primordial power spectrum and the forward model (e.g. matter density, massive neutrinos etc.). Similarly, one can work directly with the cosmological parameters instead of the bandpowers. We can unify the methods of previous sections by treating all these cases with the same formalism.
Using in equation 4.20 (where in general the sum is now over all modes ), with \hbox{\boldmath{\Pi}}_{l} from equation 4.1 (replacing with ), and defined in equation 5.2 (replacing with ), the total response is defined as their sum . We can now summarize the general analysis as following:
- Use the fiducial model to create a realization of the model \hbox{\boldmath{s}}_{s} (or multiple realizations if needed), using methods described in previous section. Forward model \hbox{\boldmath{s}}_{s} to the data domain to obtain \hbox{\boldmath{f}}(\hbox{\boldmath{s}}), and add noise to create a simulated realization of the data \hbox{\boldmath{d}}_{s+n}. Perform optimization to obtain reconstructed MAP \hat{\hbox{\boldmath{s}}}_{s+n}. Forward model again and compute derivatives with respect to forward model parameters (equation 5.3). Define the noise bias as
[TABLE]
- Perturb parameters by , one at a time. The values of are to be small enough so that they give a good estimate of the derivative, and large enough to reduce the roundoff errors. If changes the initial power spectrum compute the derivative analytically to generate a new power spectrum \hbox{\boldmath{S}}^{\rm fid}+(\partial\hbox{\boldmath{S}}^{\rm fid}/\partial p_{l})\Delta p_{l}. Generate modes with mode amplitude
[TABLE]
If changes the forward model include that in the forward model to generate the simulated data. Add the same noise as for the fiducial simulation. Perform the optimization to obtain MAP \hat{\hbox{\boldmath{s}}}_{l,s+n}, starting from \hat{\hbox{\boldmath{s}}}_{s+n} for a fast convergence. Define the Fisher matrix elements as
[TABLE]
- Perform optimization on the data using the fiducial model to obtain MAP \hat{\hbox{\boldmath{s}}}. Forward model using the fiducial model and compute the derivatives with respect to all the parameters in the forward model. The parameter estimators relative to fiducial model \Delta\hat{\hbox{\boldmath{p}}} are
[TABLE]
If the parameters are significantly different from zero (i.e. away from the fiducial model) repeat the procedure around the new fiducial model.
The covariance matrix for \Delta\hat{\hbox{\boldmath{p}}}, both the cosmology and the model parameters, is given by the inverse Fisher matrix. The solution can be multiplied by another choice of matrix . The estimator will be unbiased, but it may not be the minimum variance if some of the assumptions are violated. Additional evaluations are also necessary if one wishes to explore the posterior surface of \Delta\hat{\hbox{\boldmath{p}}} beyond the multi-variate gaussian approximation given by the Fisher matrix.
5.3 Primordial non-gaussianity
The procedure described above can also be applied to primordial non-gaussianity, which can be induced in the early universe, for example during the inflationary phase. In a typical scenario this creates a quadratic coupling between the modes, although higher order couplings are also possible. For example, the popular local non-gaussianity model can be written as [24]
[TABLE]
Here is the primordial gaussian potential, which can be evolved to late times using linear theory evolution, and related to linear density field through the Poisson equation and transfer function describing growth of structure through radiation and matter epochs. We can Fourier transform the equation above and convert to an observable density field (here meant to be linear density field, although we can also use other observables, such as the CMB temperature fluctuations) to obtain
[TABLE]
where T_{\Phi}(\hbox{\boldmath{k}}) is the transfer function between the initial potential \Phi(\hbox{\boldmath{k}}) and linear density modes s(\hbox{\boldmath{k}}), \delta_{\hbox{\boldmath{k}}} is the Kronecker delta being nonzero only for \hbox{\boldmath{k}}=0 mode, and is the variance \sigma^{2}=\sum_{\hbox{\boldmath{k}}^{\prime}}|(T_{\Phi}s)(\hbox{\boldmath{k}}^{\prime})|^{2}. We have written the expressions as a sum over dimensionless modes, with the appropriate volume factors inserted to convert from the usual modes with dimension of square root of volume.
One can view equation 5.18 as part of a forward model, where one first modifies the initial conditions using a deterministic expression of equation 5.17, with one free parameter , and then uses the forward model (such as an N-body simulation) to evolve the system forward into the data domain. One can therefore use equation 5.10) as the estimator, with the derivative in equation 5.3 obtained by a numerical or analytic evaluation of the forward model at the MAP for two different values of . From equation 5.18 we have
[TABLE]
which can be inserted into equation 5.10 (since we used dimensionless version of modes , the power spectrum should also be made dimensionless by inserting appropriate factor of volume). This gives us the optimal estimator for , and agrees in the limit of no noise and with the estimator given in [25] (expanding the estimator around ), which was shown to be optimal in this limit. Note that the reconstructed modes \hat{\hbox{\boldmath{s}}} are used in our estimator in equation 5.19. This means the noise dominated modes will be suppressed, which is equivalent to inverse covariance matrix weighting of the data, guaranteeing the analysis is optimal in the presence of both the noise and the sampling variance. As argued in equation 5.3, one may be able to evaluate the derivative of equation 5.19 numerically, by running an simulation using \hat{\hbox{\boldmath{s}}}. In the low noise limit this should work well, while its properties in the high noise regime need to be explored further.
This method is an efficient implementation of N-point function evaluation, which when done brute force involves a summation of all possible products over N modes, where in the case of example in equation 5.17. In our method we need a single summation over all the modes (and a simulation to get the derivative). This is still a 3-point function calculation, since the derivative with respect to in equation 5.10 is quadratic in the field , and we multiply it with the field itself.
One must also evaluate the Fisher matrix to obtain the errors and the correlation of estimator with all other parameters, as well as noise bias, following the steps in previous subsection. Specifically, we run a simulation with , and derive the response of the estimator in equation 5.10 as the difference between this estimator and the one for the fiducial model (which should be close to 0), which gives us the appropriate normalization . By the arguments given in section 4 this also provides the inverse covariance matrix. Explicit evaluation of the likelihood function to quadratic order in and evaluating the ensemble average of its second derivative confirms this result [25]. It is straight-forward to evaluate the other parameter response to the difference between and fiducial simulation. Among these parameters are the bandpowers, as well as the mean density parameter, which has to be determined from the data. Our procedure thus gives the minimum variance estimator, together with the estimate of the corresponding Fisher matrix.
The forward model parameters may not suffer from the sampling variance errors in the same way as the bandpower parameters. This depends on the full Fisher matrix and the mixing between the forward model parameters and the bandpower parameters. The easiest way to see the sampling variance cancellation aspect is to have two different tracers, for example one that responds to on large scales and one that does not, as in the case of local at the linear level, where these can be viewed as galaxies with bias and dark matter, respectively. In this case the dark matter will uniquely determine initial large scale modes , leaving the galaxies to determine limited only by the noise of the former tracer. This is the basis of the sampling variance cancellation technique as applied to local [26]. Equations above should optimally extract all the information present in the data, automatically including all the information. At higher orders the sampling variance cancellation is implemented in the mode couplings: for example, a given large scale mode can have a variable coupling strength as a function of small scale modes in the presence of primordial non-gaussianity, allowing to identify the primordial non-gaussianity independent of the amplitude of the long wavelength mode. This is limited by the variance of the small scale modes. The method is thus not really sampling variance free, but limited by the sampling variance of small scales modes, which however may be small because of the high number of modes.
6 Application to dark matter simulations
We have tested the bandpower method described in section 4 on both 2-LPT simulations and on N-body simulations. For the latter we use a cubic box of size with a grid, using the publicly available code FastPM [27] for simulating the nonlinear structure formation and the public code Cosmo++ [28] for optimization. We have used both the full N-body simulation and a simplified second order Lagrangian perturbation theory (2-LPT) for our truth, but below we only show results for the 2-LPT, since the results were almost indistinguishable. We use standard Planck cosmology for all cosmological parameters and the outputs are evaluated at .
Generating a single realization of a full simulation using FastPM is fast. However, for the reconstruction one also needs a gradient with respect to all initial modes (equation 3.4). This can be done using backward propagation [29], but is expensive, both because FastPM is slower than 2-LPT, and because gradient calculation is expensive and memory intensive. Moreover, as we show below, the optimization requires hundreds of calls to compute the gradient, making the full gradient considerably more expensive to compute.
For the purpose of reconstruction in this paper we thus use a simpler method - we simply use 2-LPT multiplied with an appropriate transfer function. Specifically, we first compute a transfer function by simulating the density fields using PM and 2-LPT and taking the square root of the ratio of the power spectra. Then the gradients are performed using 2-LPT multiplied by the transfer function. The 2-LPT gradients are given in appendix D. If the cross-correlation coefficient between full N-body and 2-LPT is close to unity then there is no loss of information in this process. In practice, at the cross-correlation coefficient is 0.98 [30], and this approximation is extremely accurate over the range of we are interested in here.
In this work we fix the transfer function, under the assumption that it does not vary with the input power spectrum. A full gradient implementation would remove this transfer function step. The simulated fields are smoothed with a gaussian kernel of size , which roughly corresponds to the pixel size. We add white noise to the truth simulation in real space. While smoothing was required to achieve convergence in our tests, with more accurate gradients one should be able to remove this step, and this is something we plan to pursue in the future. Smoothing can often be justified as a way to implement finite resolution of the experiment, although here it should be viewed simply as a tool to improve the convergence.
Note that since we use a full N-body simulation for the truth, but we only use 2-LPT times the transfer function for the optimization, we cannot claim that the procedure extracts all the information optimally. For this reason we also performed the same analysis using 2-LPT as the truth: in this case no transfer function is needed, and within the context of 2-LPT dynamics this procedure is guaranteed to give the optimal power spectrum reconstruction. However, we observe almost no difference between the two cases, a consequence of smoothing, which eliminates signal relative to noise at high , where the cross-correlation coefficient drops below unity. We only present results from the 2-LPT simulation as the truth.
The forward model depends not only on the initial modes, but also on the assumed cosmology. Specifically, in the context of CDM one needs to assume matter density , which determines the growth rate as a function of redshift. In this paper we simplify the analysis by ignoring this dependence, and we assume is known. A more general analysis would explore the likelihood surface as a function of as well, and we plan to explore this in the future.
We show a slice of our truth simulation, with and without noise, in Fig. 1. The nonlinear power spectrum of the truth simulation is shown in Fig. 2. As a consequence of smoothing the data becomes noise dominated for .
Let us first discuss the reconstruction of the initial density field. Our optimizer starts at the initial field of zero, and converges in iterations, and at each iteration we evaluate one or two line search steps, with a total of about 300 function evaluations. A slice through the original and the reconstructed initial density fields are shown in Fig. 3. For comparison, we also show the reconstruction after and iterations. The corresponding nonlinear density fields are shown in Fig. 4. As we can see, the reconstruction works very well on large scales. However, the reconstruction is missing the small scale structure, which is expected, since the data is noise dominated on small scales and the prior, the first term in Eq. 3.2, drives the small scale modes to zero. The prior thus acts as a regularizer: in the absence of any real information in the MAP it sets the modes to zero. Lack of small scale modes is much more evident in the initial linear density field of Fig. 3. However, it can also be noticed in the nonlinear field in Fig. 4. So far we have not encountered any problems with local minima: the procedure always converges to the same solution irrespective of the initial starting point. It remains to be seen whether that will remain the case when we have to deal with strongly nonlinear regime.
Note that even after iterations the structures on large scales are already reconstructed fairly well. For this reason we would expect baryonic acoustic oscillations (BAO) to be reconstructed well after a few iterations only, similar to the BAO reconstruction methods [31, 32, 33]. After iterations the reconstructed field is indistinguishable by eye from the final iteration.
To show the convergence rate of the optimizer we have plotted as a function of number of iterations in Fig. 5. The convergence, here defined as a change of less than 1, is achieved after 236 steps. After 150 steps the is about 1000 above the minimum. While this may sound a lot, the number of modes is , so per degree of freedom this is a small change, and explains why visually the maps do not significantly change from 100 iterations to the final reconstructed map. We discuss this issue further below.
To compare the reconstructed density field with the original truth field quantitatively, we plot the cross correlation and the transfer function between these two fields in Fig. 6. The transfer function is defined as the square root of the ratio of the power spectra, while the cross correlation is the ratio of the cross-power to the square root of the product of the individual power spectra. As expected, both the transfer function and the cross-correlation are very close to unity for large scale modes, and gradually drop down to [math] for . This is expected, since the data becomes noise dominated beyond . The effect of noise is more pronounced for initial density field. This is because high nonlinear density field is still significantly determined by low linear density, as a consequence of nonlinear mode coupling and power transfer.
Let us now discuss the power spectrum reconstruction. We use equal width bins in . For the reconstruction we use a fiducial power spectrum which is a no-wiggle version of the power spectrum for the truth simulation. This means that the BAO wiggles have been removed. One of the goals of the initial power spectrum reconstruction is to reconstruct the BAO wiggles without putting them into the fiducial model, and to show that the reconstruction completely removes nonlinear smoothing it is best if the fiducial power spectrum has no wiggles at all.
We use equations 4.46 and 4.41 to compute the Fisher matrix and noise bias. Specifically, we first create a fiducial power spectrum mock simulation as in equation 4.40, and perform exactly the same optimization steps as for the ”truth” simulation. We next inject small amount of power into a single bandpower and measure the response as given by all the other bandpowers. Here we start the optimization at the values given by the fiducial model simulation, and as a consequence we only need about 50 calls to converge. To make the procedure even faster we could do this for every N-th bandpower (e.g. N=5), since the mean transfer function squared of equation 4.47 is a smooth function of bandpower and can be interpolated across them. Once we have computed the Fisher matrix we also compute the noise bias using equation 4.41.
We plot the normalized Fisher matrix (equation 4.47) in figure 7. As we can see, the window functions are peaked at the diagonal, and no significant mixing of modes occurs. This is expected when noise is low, since in that case we are reconstructing true modes. Note that we use periodic box, so there is no window function associated with the survey geometry itself: the nonlinear modes are orthogonal, so any window function deviation from a delta function would be due to the nonlinear mode coupling to noise, creating off-diagonal elements. However, we observe no such mode coupling and no off-diagonal window function even for high noise. This result is intriguing and suggests that nonlinear evolution coupled to noise does not create much off-diagonal coupling. Our test case had a very low level of noise, and it remains to be seen if the same result is confirmed with higher levels of noise. Survey mask and other effects will create a window, which will be localized to the width of the bin of order , where is the size of the survey. It may thus be possible to separate geometry effects from the nonlinear noise coupling effects, and develop rapid methods for the estimation of the covariance matrix. Clearly this needs more detailed investigation, as it could greatly simplify the calculation of the window function and covariance matrix.
Finally, we plot the original and the reconstructed power spectra in Fig. 8, all relative to no wiggle fiducial power spectrum. We show the original truth power spectrum (red solid curve), the power spectrum of the reconstructed MAP initial density field (grey dashed curve), and our reconstructed power spectrum with error bars (black points). We can see that our reconstructed “best fit” MAP density field has a power spectrum that matches very well the “truth” power spectrum on large scales, but starts to decrease at small scales. As we discussed above, this is due to noise: in the presence of noise the optimal reconstruction sends the modes to zero. However, our reconstructed power spectrum (black solid curve) is able to determine the amount of lost power, and make the estimator unbiased. Indeed, it matches very well the model power spectrum on small scales, which is an artifact of the fact that we used the no wiggle version of the same power spectrum for bias term as for the truth. In reality we would not have this luxury so the results will be biased at high , however not more than the computed errors, which properly account for this effect: we see that the errors get very large for , since the transfer function of initial modes becomes very close to 0 and the MAP reconstruction is very close to zero.
As we can see clearly from Fig. 8, our reconstruction does reproduce the linear BAO wiggles, which match very well the wiggles in the truth power spectrum. For the nonlinear power spectrum the BAO wiggles are damped, as shown in the plot with the blue dashed curve. We can see from the plot that we are reconstructing the BAO wiggles of linear power spectrum, and not of the nonlinear power spectrum, achieving very good reconstruction. The last two wiggle with are not reconstructed, because there the noise and smoothing already completely destroy the signal, and since we did not have BAO wiggles in our fiducial model they are not reproduced. The size of the errors reflects this, and the last BAO wiggle that is detectable is around .
As discussed in section 5 it is not guaranteed that our procedure gives a reliable Fisher matrix that can be used as the covariance matrix, although this is guaranteed in the linear regime and in the low noise regime. We have computed goodness of fit value of the estimators relative to the truth \Delta\hbox{\boldmath{\Theta}}^{\rm truth} (which is not 0 because our assumed fiducial model is not the truth model)
[TABLE]
using 40 bandpowers up to . The reduced of 0.94 suggests the covariance matrix is reliable.
We have also compared the diagonal terms of the Fisher matrix to the simple approximation of equation 4.31, finding good agreement between the two methods (figure 9). This, and the nearly diagonal window (figure 7) suggests that one may be able to estimate a reliable covariance matrix even without doing any simulations (more precisely, with a single simulation to compute \hbox{\boldmath{b}}_{l}, making sure the estimator is unbiased). This clearly deserves a more detailed investigation, in particular in the presence of a survey mask and other complications, but is beyond the goal of this paper.
6.1 Convergence and computational cost
Let us look at the computational costs of our method. As mentioned above the cost of a single optimization procedure is about 230 forward model calls, each being a few FFTs (to compute 2-LPT) of the mesh. On the data this is done once. For the simulations we have to do the same for the fiducial simulation realizations, and then for each of the Fisher matrix rows, where we inject a small amount of power into a single bandpower on top of the fiducial power. The latter would thus appear to dominate the computational cost, but because the Fisher matrix can be interpolated across the bandpowers (we only evaluate every fifth row), and because we can start the optimization from the fiducial model solution, which reduces significantly the number of iterations, in practice this entire operation is comparable in cost to the optimization of fiducial model simulation, or of the data. To achieve required precision one also needs to repeat the procedure until required precision is achieved. In our tests we have not found the need to average over many realizations.
Another possible acceleration is to extract the results without reaching the global minimum. As figure 3 shows the solutions between 100 iterations and final convergence are almost indistinguishable, since the optimizer is changing a large number of modes by a small amount that cannot change the solution significantly. This result is confirmed by the actual value of as a function of iteration number, which hardly changes at higher iterations even when the loss function is still changing significantly. In general, if we have modes and the optimizer has determined each to within , where is the error for individual mode, then for there is no loss relative to perfect optimization solution, in the sense that the total error is given by . As long as the error is stochastic the same will be true for the bandpower error, where we average several modes. But this additional noise adds to the overall loss function, and for large this can be a large number even if is small. For example only increases the overall error by 0.5%, suggesting that even the optimization solutions which are away from the global minimum by a large number may still be good enough. Another way to understand this argument is through the sampling perspective: when sampling each mode is on average away from MAP, and most of the sampling is done in a narrow range of loss function being above the global minimum. The sampler never descends down to the minimum, and its exact position is irrelevant for the posterior distribution of the modes, which is ultimately what determines the solution. Figure 5 suggests that we could stop the optimization after 50-100 steps and still satisfy the condition that . In practice we observe strong correlations in the modes of a single bandpower along the optimization path, so the argument above is too naive and one must develop a more robust criterion when to stop. We also observe this argument to fail when computing Fisher matrix elements, where small differences between solutions are relatively more important.
A related issue is the question of convergence to the global minimum. The loss function is only convex for the linear problem. For the nonlinear problem it becomes non-convex, meaning that there can be multiple local minima. In the mildly nonlinear regime this will not happen, but in deeply nonlinear regime this is very likely to happen. Changing the starting point of optimization, or using a stochastic descent which tries different paths down the hill are typical approaches to this problem, but a general solution is difficult. However, as argued above a local minimum is only problematic if it is very far from the global minimum. A situation where there are many high quality local minima at the bottom of the loss function is not of much concern, since all those solutions are in some sense equivalent, specially if they all give the same transfer function.
Another numerical issue are saddle points. These are numerically problematic, since they are attractors for Newton’s method optimizers: they may give an impression that convergence has been reached, only for the optimizer to suddenly start reducing its loss function again once it has evolved past the saddle point. We have observed this in the configurations with high noise. Adding noise or momentum to the gradient, taking absolute value of the Hessian eigenvalues or some other way to change the Hessian make these issues less problematic [34]. More work is needed to develop optimal optimization methods for our problem.
In any case, the large number of forward model evaluations, in the hundreds to thousands, and the need to compute the gradients puts a high demand on its computational cost, and approximate methods like 2-LPT or fast simulations like FastPM [27] are crucial for the success of this program. We note that FastPM with 10 time steps is 7 times slower than 2-LPT, while its transfer function with full N- body simulations at is 1.02, relative to 1.2 for 2-LPT, so its computational cost is not prohibitive, while dramatically improving the accuracy. However, one also needs to demonstrate that its gradients can be implemented accurately and can be evaluated. In practice the choice of forward model will depend on the level of noise and resolution in the data, with low resolution/high noise data requiring only 2-LPT and high resolution/low noise data requiring N-body simulations with low number of time steps. It remains to be seen whether this will change with a more realistic survey geometry, where window functions will be less compact. It may be possible to reduce this number further if the convergence of the optimizer is accelerated. We note that in our tests a simple conjugate gradient optimizer was a factor of 2-3 times slower than L-BFGS, which is the reason we only present the latter in these tests. We plan to investigate some of these acceleration issues further in the future.
7 Conclusions
The purpose of this paper is to develop a formalism to optimally extract cosmological information in the nonlinear regime, and to develop numerical methods that are computationally feasible and preserve the near optimality of the analysis, and allow an implementation of this formalism both in the linear and nonlinear regime. Traditional CMB/LSS analyses start with the power spectrum or correlation function of the observed field (such as CMB temperature and polarization, weak lensing shear or galaxy positions), and often end there. Even in this case the two-point function analysis is typically not optimal, because the data are not inverse covariance matrix weighted.
In the nonlinear regime the LSS evolution creates non-gaussian signatures. This is a generic term that contains a lot of different manifestations: they range from the higher order correlations such as bispectrum, to existence of special objects like the density peaks (e.g. clusters), to non-trivial topologies such as sheets, filaments and voids. The problem with these descriptions is that they are not exhaustive, i.e. it is not possible to prove that any specific combination extracts all the information. Indeed, often it is not even clear that these additional statistics increase the amount of information relative to the power spectrum at all. When combining different statistics one must also obtain their joint covariance matrix, which is a challenge in itself, since it typically requires a large number of mock data realizations (parametric bootstraps).
In this paper we present an (approximately) Bayesian analysis that attempts to nearly optimally extract the information contained in the data given a nonlinear structure formation models of LSS. By writing down an exact probabilistic model the optimality of the analysis is guaranteed as long as all the steps are solved exactly. In our approach we do not attempt to produce an exact solution, which is prohibitively expensive, but instead we develop a series of controlled approximations that we argue approximately preserve the optimality of the solution. Should these approximations turn out to be insufficient one can attempt to generalize the approach used here.
Specifically, our approach is to expand the likelihood around the maximum a posterior (MAP) solution for the initial modes, and integrate out the modes around this solution. In our specific implementation we assume existence of a single global MAP solution. This assumption is guaranteed in the linear model, and we argue that it is likely to be the valid at least in the mildly nonlinear regime, which is our primary focus. We do not analytically compute the integration volumes, which correspond to the product of all variances of the modes, but instead rely on simulations to compute their derivative with respect to the parameters. We do this using simulations, and the realization dependence of the variance is lost. We have argued that at the linear level this is a valid and lossless approach. We also expect any corrections to be small in the low noise limit, and in general we expect this to be a very good approximation. In this context it should be stated that in the presence of multiple local minima the solutions should not be expanded around the global MAP, but around the one that maximizes the loss function times the log of posterior volume (equation 4.10), which differs from what MAP maximizes by \det\tilde{\hbox{\boldmath{D}}}^{-1/2}. While we have argued that it is likely that one does not need to average over multiple local maxima, this does not imply that finding the one that maximizes MAP (or volume corrected MAP) is easy. A proliferation of lower quality (higher loss function) saddle points can still make convergence slow, and there is no guarantee that the optimization will find a global minimum, as opposed to a local minimum far from the best solution. We have argued that if these local maxima are close to the global MAP and all have the same transfer function it does not affect the unbiasedness of the solution.
Similarly, the optimality of the method is only guaranteed if the forward model is exact, even if it is unknown (exact in the sense that the nuisance parameters always cover the truth). Formally one adds more and more forward model parameters, and marginalizes over them. While it is easy to achieve this on large scales, it becomes increasingly more difficult to do so on small scales where many uncertain physical processes can affect the data, and one needs to include all possible variations consistent with the allowed physical processes, including physical priors on these parameters.
We argue that the bandpowers of optimally reconstructed initial power spectrum are lossless summary statistics, at least for the simplified case where the forward model connecting initial modes to the data is unique. We provide a derivation of the bandpower estimators, the corresponding Fisher matrix and noise bias, which fully quantify the information content of LSS data. Note that our approach does not involve any sampling of the modes, and we believe that current implementations of sampling procedure to be too slow to be of practical use in these applications.
We then generalize the analysis allowing also for parameters that control the mapping from the initial modes to the final data, both cosmology parameters such as matter density or massive neutrinos, and astrophysical modeling parameters such as galaxy bias, shear bias etc. We expect that these parameters in general are not degenerate with the power spectrum amplitude parameters, even if they may be at the linear level: in the nonlinear regime this degeneracy can be broken, since nonlinear gravitational effects cannot be exactly mimicked by the nuisance parameters. This allows for possibility of marginalizing over the nuisance parameters internally, without the need for external constraints. We show that one can treat primordial non-gaussianity as a forward model parameter and present an optimal analysis method to extract it. Our method involves only a single sum over the data, multiplied by a response function numerically derived at each data point.
For linear models, such as the CMB or the LSS on large scales, where the relation between the data and the modes is a simple Fourier transform, the MAP is guaranteed to be at the global maximum with a gaussian posterior and this guarantees optimality [7]. Even in this case we expect the numerical methods developed here, optimization and fast Fisher and noise matrix evaluations, to be useful in applying these estimators to the real data. We will present our results of applying these methods to weak lensing and CMB applications elsewhere.
This paper focuses on the theoretical formalism, deriving the summary statistics in the context of nonlinear structure formation models. We implement the procedure on a simplified, yet nontrivial test case, where we use 2-LPT or N-body simulations as a forward model to generate a nonlinear density field, while we use simplified 2-LPT times the transfer function for gradient evaluations during the optimization. We show that the resulting solution recovers the linear baryonic acoustic oscillations, achieving a near optimal BAO reconstruction. Since our gradients are not exact we do not claim that the reconstruction is optimal, although this is very likely to be the case for BAO (since 2-LPT is expected to fully model nonlinear BAO evolution). Our reconstructed power spectrum is unbiased and comes with the bandpower window matrix and the (disconnected) covariance matrix that describes well the solution.
Some of the more realistic effects that remain to be included are incomplete data coverage, nuisance parameters such as the biasing of galaxies, and baryonic effects on the matter distribution. Equally importantly, while for the weak lensing the noise model is simple, for galaxy data noise probability distribution needs to be understood better, and sub Poisson noise should be achievable [35]. When applying to weak lensing one additional complication is the elongated radial projection window, which does not allow for a full 3-d reconstruction even when tomographic reconstruction is used, and will lead to significant degeneracies in the radial direction.
Ultimately we expect the methods developed here to superseed the traditional data analysis methods in cosmology. When this will be achieved depends mostly on the remaining technical and computational challenges that need to be addressed. Among these are fast evaluations of derivatives of N-body simulations, where we expect FastPM code [27], which achieves percent level agreement of halo clustering against high resolution simulations with as few as five time steps, to be particularly suitable for this purpose. Another challenge is the ability to perform gradients of galaxy density field, as well as gradients of forward model parameters evaluated on MAP. We plan to investigate some of these issues in the future. Given the higher computational cost of this approach, and various numerical challenges in the nonlinear regime, there may still be a long road ahead before we achieve this goal, but the payoff will be an increased precision of cosmological constraints and (near) optimality of the analysis.
Acknowledgments
We thank Emanuele Castorina, Simone Ferraro, Ryan Giordano, Patrick McDonald, Martin White and Matias Zaldarriaga for useful discussions. We acknowledge support of NASA grant NNX15AL17G. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.
Appendix A Appendix: Beyond the gaussian posterior approximation
The power spectrum bandpower procedure described in the main text is explicitly unbiased, and provides bandpowers and their covariance matrix under the gaussian approximation for the posterior. In a typical application to LSS the number of modes in a bandpower is large, so that by the central limit theorem the gaussian approximation to the likelihood in equation 4.17 suffices. When the number of modes inside the bandpower is small this is no longer the case: a model with a large amount of power relative to the best estimate also has a larger sampling variance error, making this model more likely than the gaussian approximation would suggest. Since this situation is most likely to occur on large scales, where the modes are linear, we will consider linear case only. For the simplest linear model with complete coverage one can write the likelihood exactly in terms of () as (inverse) Wishart distribution [36]. For our purpose we need ,
[TABLE]
where
[TABLE]
and , where is the number of modes in the bin. For convenience we have normalized the log-likelihood to zero when .
A simple generalization for the optimal quadratic estimator is to use equation A.1 together with the uncorrelated version of the estimator in equation 4.34. Since the product of inverse square root of Fisher matrix with the Fisher matrix is the square root of the Fisher matrix, this requires its evaluation. The best option is for matrix {\hbox{\boldmath{F}}}^{-1/2} to be symmetric, which can be evaluated by diagonalizing {\hbox{\boldmath{F}}}={\hbox{\boldmath{ZDZ}}^{{\dagger}}}, and then evaluating the square root of the diagonal elements to give {\hbox{\boldmath{F}}}^{-1/2}=\hbox{\boldmath{ZD}}^{-1/2}\hbox{\boldmath{Z}}^{{\dagger}}. This will make the window function narrower (for narrow windows, by Taylor expansion, the off-diagonal terms will be approximately halved). Once diagonalized, we can use
[TABLE]
One can interpret this term as giving an effective number of modes measured in bandpower (divided by 2), given the noisy and incomplete data. The effective number of modes is not only a function of volume coverage, but also of noise and fiducial power: pixels that are noisy relative to the fiducial power are downweighted, reducing the effective volume. Note also that the signal to noise in bandpower is lower than , because of noise bias ,
[TABLE]
Appendix B Appendix: Alternative approaches to
Fisher matrix and noise bias
Fisher matrix is defined as the response of the modes in a bandpower bin to modes in a bandpower bin . To evaluate this response we can compute the cross-correlation of one mode of bin to the modes in bin . We label this mode with index , and there are modes inside the bin . We define a transfer function between initial and reconstructed modes as a cross-correlation divided by auto-correlation
[TABLE]
When maximally correlated and without noise this equals unity for . By defining this as a transfer function and dividing with the actual power we reduce the sampling variance. We can also choose the actual power for each mode to be fixed to the fiducial power, reducing fluctuations further. This is different from the typical mock (Monte Carlo) simulations, which require more realizations to converge. For the same reason our noise fluctuations are forced to have the same variance via random realizations. If necessary the procedure is repeated on several realizations, until we achieve the required convergence.
Fisher matrix is then given by
[TABLE]
To evaluate the noise bias we use the difference between the auto-correlation, which contains noise contribution, and cross-correlation, which does not,
[TABLE]
where
[TABLE]
From this we finally obtain the unbiased estimator by adding over all modes within the bandpower
[TABLE]
This is our minimum variance estimator in the sense of equation 4.32, i.e. with definition. It is manifestly unbiased, in the sense that (equation 4.36). It is quadratic in the optimal field reconstruction , which is in turn nonlinearly related to the data . The answer is given by the optimal field reconstruction, corrected for the noise bias and properly normalized.
The procedure is useful to clarify the nature of noise bias and Fisher matrix, but its convergence is very slow. The reason is that we have to compute cross-correlations between all the modes, then square it, and this process will be biased positive even if there are no correlations between the modes, and will only slowly converge to the correct value with the number of simulations. Moreover, the calculation of the Fisher matrix using B.2 is quadratic in the total number of modes, which can become prohibitive for or more modes. In the main text we developed an alternative approach that converges faster.
A second alternative method takes advantage of the Fisher matrix interpretation as the covariance matrix. Suppose that the Fisher matrix is diagonal. This is true for uniform noise and periodic box simulations used in this paper (figure 7). In this case we have from equation 4.29
[TABLE]
where is the variance of the quantity, ie . Assuming are gaussian distributed we can use Wick’s theorem . This, together with equation refminvar gives equation 4.31. As shown in figure 9 this works remarkably well. It is computed entirely from the data, so one can argue that this is the actual curvature matrix needed in a Bayesian analysis, as opposed to its ensemble average, the Fisher matrix. This equation needs to be generalized for realistic surveys with a survey mask, gaps, variable noise etc. to its matrix form, but we will not pursue it here.
Appendix C Appendix: Analysis in the presence of external additive nuisance parameters
Minimum variance power spectrum estimator in the form we presented in the main text has to be modified when the data are contaminated or do not contain information for some modes [4]. For example, in LSS galaxy surveys the mean density is unknown and obtained from the measured data itself, which means that low modes cannot be measured independently. The mean density mode needs to be assigned a large variance, which in turn forces its reconstructed value to zero. Because of the finite survey size a uniform density across the survey can be created from \hbox{\boldmath{k}}>0 modes: the window function of the \hbox{\boldmath{k}}=0 determines this mode mixing. Other contaminants exist as well, such as dust extinction, star-galaxy separation etc. Assuming we know their spatial structure, we can model them as a known data vector \hbox{\boldmath{L}}(\hbox{\boldmath{q}}), with an unknown amplitude . One wants to make sure that the estimator remains unbiased and that the contaminated modes do not affect the power spectrum estimator.
In some instances we want to remove these components from the data without knowing their actual values, in other cases we may be interested in their best reconstructed values. From the point of view of this paper we want to marginalize over these modes by integrating over their probability distribution. Let us discuss the case of simultaneous reconstruction and external parameter determination. Following [4] and generalizing to the nonlinear case we can write the data vector in the presence of external parameters as \hbox{\boldmath{d}}=\hbox{\boldmath{f(s)}}+\hbox{\boldmath{L}}(\hbox{\boldmath{q}})+\hbox{\boldmath{n}}, where is the nonlinear mapping and is the number of external parameters. For example, when determining the mean density one can write \hbox{\boldmath{L}}(q)=(1,1,...,1)q and if the data express galaxy numbers on a pixelized grid corresponds to the mean number of galaxies per pixel. If the number density is evolving with redshift then one must determine the relative evolution from the model of the galaxy selection, so that only the overall normalization needs to be determined.
We assume a gaussian prior for
[TABLE]
with the prior covariance matrix of the nuisance parameters. Often we have no prior knowledge of these nuisance parameters, in which case the prior is simply flat. The full posterior is
[TABLE]
We can now simultaneously maximize the posterior probability for and . This gives \hat{\hbox{\boldmath{q}}} and \hat{\hbox{\boldmath{s}}}, and the corresponding curvature matrix now includes properly the correlations between the nuisance parameters and the modes . If we want to know the probability distribution of the data independent of external parameters and the underlying field we marginalize over these parameters by integrating P_{q}(\hbox{\boldmath{q}})P_{s}(\hbox{\boldmath{s}})P_{n}[\hbox{\boldmath{d}}-(\hbox{\boldmath{f(s)}}+\hbox{\boldmath{L(q)}})] over d^{M_{q}}\hbox{\boldmath{q}}d^{M}\hbox{\boldmath{s}}, as we did in equation 4.10. The resulting posterior as a function of power spectrum will properly account for all the correlations between the nuisance parameters and the power spectrum. Since our estimator is built out of a sequence of optimization steps, there is no difference to the equations we wrote above: we are simply minimizing for a few more parameters, with no significant cost increase given the large number of modes.
When this is applied to the linear model one can write an explicit form of the covariance matrix, and with the use of Sherman-Morrison-Woodbury formula one can show that this corresponds to addition of S_{q}\hbox{\boldmath{LL}}^{{\dagger}} for each parameter to the covariance matrix of the data [4]. However, one does not need to restrict to linear nuisance parameter models, nor does one need to write down the covariance matrix: since our procedure never requires an explicit form of the curvature matrix or its inverse, there is no need to know its form in the presence of these contaminated modes: optimization procedure and the corresponding construction of the noise bias and Fisher matrix automatically give the desired result.
Appendix D Appendix: 2-LPT Derivative
The optimization involves the evaluation of the derivative (equation 3.5) term at each step
[TABLE]
where we have ignored the derivative of the prior term for the time being and is defined as the gradient of the forward model at every position with respect to every initial mode.
In the main text, we show the results using 2-LPT as the forward model. In this appendix we present the evaluation of its derivatives in this subsection. The corresponding derivatives at 1-LPT level have been presented in [13] and at 2-LPT level in [37], albeit in a different form from the one here. We use 2-LPT because it significantly improves the cross-correlation coefficient between the full nonlinear field and the approximate 2-LPT field [30]: in the limit of cross-correlation coefficient being unity the resulting solution multiplying 2-LPT with the transfer function becomes exact, and the error is given by the deviation of the cross-correlation coefficient from unity. As we discuss in the text, all of our 2-LPT fields are multiplied with the transfer function, which we will omit here for the sake of clarity. To make it more illustrative, we replace the generic notation used in the main text with more conventional notation for our specific model, i.e. we replace the density modes of interest with \hbox{\boldmath{\delta}}_{j}({\bf k}) where the subscript is [math] or for the real and complex parts of a Hermitian field respectively. We write the noise matrix as N_{ij}=\sigma^{2}(\hbox{\boldmath{r}}_{i})\delta^{D}_{ij}. Thus the derivative explicitly summed over all the positions is
[TABLE]
where \rho_{d}({\bf k}_{1})=\mathcal{F}\frac{[\hbox{\boldmath{d}}({\bf x})-\hbox{\boldmath{f_{\delta}}}({{\bf x}})]}{\sigma^{2}({\bf x})} and is the Fourier transform (). The forward model with CIC interpolation of particles from their Eulerian position is,
[TABLE]
where is the smoothing kernel and is the CIC kernel assignment. The 2 LPT forward model calculates the Eulerian positions of the particles to second order in Lagrangian displacement, . These displacements in Fourier domain are
[TABLE]
This second order displacement is simply the conventional definition
, written as convolutions in the Fourier space.
Thus, the complete derivative can be expressed as independent derivative of first and second order displacements
[TABLE]
where, we have defined
[TABLE]
which can be evaluated using simple Fourier transforms and re-griddings. Then, to first order (1-LPT), this leads to
[TABLE]
which leads to the final expression for the most modes
[TABLE]
while the factor of drops out for the self-conjugating (zero and the nyquist) modes of the hermitian field.
At second order (2-LPT), we have
[TABLE]
where in the first expression we use result from first order derivative and in the second step we are using the symmetry between and . As done above for , we have:
[TABLE]
where simplifications have been made using symmetry between different dummy variables and Hermitian properties of and .
The final expression is a convolution that can be evaluated as a product in the configuration space after defining 6 functions and a function as:
[TABLE]
leading to the second order terms
[TABLE]
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] A. J. S. Hamilton, MNRAS 289 , 285 (Aug. 1997).
- 2[2] M. Tegmark, Phys. Rev. D 55 , 5895 (May 1997).
- 3[3] J. R. Bond, A. H. Jaffe, and L. Knox, Phys. Rev. D 57 , 2117 (Feb. 1998).
- 4[4] G. B. Rybicki and W. H. Press, Ap J 398 , 169 (Oct. 1992).
- 5[5] N. Padmanabhan, U. Seljak, and U. L. Pen, New Astronomy 8 , 581 (2003).
- 6[6] F. Elsner and B. D. Wandelt, A&A 549 , A 111, A 111 (Jan. 2013), 1210.4931.
- 7[7] U. Seljak, Ap J 503 , 492 (Aug. 1998).
- 8[8] B. D. Wandelt, D. L. Larson, and A. Lakshminarayanan, Phys. Rev. D 70 , 083511 (2004), astro-ph/0310080.