Matter power spectrum using artificial neural networks

PkANN – I. Non-linear matter power spectrum interpolation through artificial neural networks

Shankar Agarwal, Filipe B. Abdalla, Hume A. Feldman, Ofer Lahav & Shaun A. Thomas Department of Physics & Astronomy, University of Kansas, Lawrence, KS 66045, USA. Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK. emails:,,,,

We investigate the interpolation of power spectra of matter fluctuations using Artificial Neural Network (PkANN). We present a new approach to confront small-scale non-linearities in the power spectrum of matter fluctuations. This ever-present and pernicious uncertainty is often the Achilles’ heel in cosmological studies and must be reduced if we are to see the advent of precision cosmology in the late-time Universe. We show that an optimally trained artificial neural network (ANN), when presented with a set of cosmological parameters ( and redshift ), can provide a worst-case error per cent (for ) fit to the non-linear matter power spectrum deduced through -body simulations, for modes up to . Our power spectrum interpolator is accurate over the parameter space. This is a significant improvement over some of the current matter power spectrum calculators. In this paper, we detail how an accurate interpolation of the matter power spectrum is achievable with only a sparsely sampled grid of cosmological parameters. Unlike large-scale -body simulations which are computationally expensive and/or infeasible, a well-trained ANN can be an extremely quick and reliable tool in interpreting cosmological observations and parameter estimation. This paper is the first in a series. In this method paper, we generate the non-linear matter power spectra using halofit and use them as mock observations to train the ANN. This work sets the foundation for Paper II, where a suite of -body simulations will be used to compute the non-linear matter power spectra at sub-per cent accuracy, in the quasi-non-linear regime . A trained ANN based on this -body suite will be released for the scientific community.

Key words: methods: numerical – cosmological parameters – cosmology: theory – large-scale structure of Universe.

1 Introduction

Studying the growth of structure and the distribution of galaxies in our Universe is a potent method for understanding both fundamental physics and the cosmological model. Measurements of the large scale structure are capable of testing our theory of gravity (Amendola, Kunz, & Sapone, 2008), distinguishing between dark energy models and probing the absolute neutrino mass scale (Thomas, Abdalla, & Lahav, 2010). It is also an immense complement to the incrementally cosmic variance limited cosmic microwave background. This complementarity explains the overwhelmingly large number of galaxy surveys on the horizon that promise to refine or even alter our understanding of the cosmos, e.g. DES (The Dark Energy Survey Collaboration, 2005), the Large Synoptic Survey Telescope (LSST) (Ivezic et al., 2008), and the Baryon Oscillation Spectroscopic Survey (BOSS) (Eisenstein et al., 2011). These surveys promise to achieve high-precision measurements of galaxy power spectrum amplitudes and offer a possibility to improve constraints on cosmological parameters including dark energy and neutrino masses. However, with this promise comes a great technical and systematic difficulty.

Arguably the most ubiquitous problem in both galaxy clustering and weak lensing surveys is that as structures collapse they evolve from being linear, for which one can solve analytically, to non-linear, for which one cannot. Using -body simulations (Heitmann et al., 2010; Agarwal & Feldman, 2011) and analytical studies inspired from perturbation theory (PT) (Scoccimarro, Zaldarriaga, & Hui, 1999; Saito, Takada, & Taruya, 2008), the non-linear effects have been shown to be significant compared to the precision levels of future surveys. A consequence of this is the uncertainty in calculating the theoretical power spectra over smaller scales and at low redshift. There is frequently a choice to either exclude – and therefore waste – the wealth of available and expensively obtained data, or to use an inaccurate procedure, which may bias and invalidate any measurement determined with anticipated precision.

At present there are several approaches to dealing with this fruitful yet frustrating regime. One is to use sophisticated -body simulations commonly produced with codes such as enzo (O’Shea et al., 2010) and gadget (Springel, 2005). The most popular non-linear prescription halofit (Smith et al., 2003) is a semi-analytical fit and has been fantastically successful. However, with larger and ever-improving state of the art of -body simulations, the non-linearities on smaller scales have been shown to be at levels higher than the ones that were used in calibrating halofit. As such, on smaller scales () the matter power spectra predicted by halofit do not match the high precision -body results well enough. If we are to perform precision cosmology it is imperative to go far beyond the levels of precision offered by current analytical approximations. An obstacle to further progress in obtaining accurate fits to underlying spectra is the vast computational demand from detailed -body simulations and a high dimensionality in the cosmological parameter space.

There have been attempts (see Bird, Viel, & Haehnelt 2012) to calibrate halofit using -body simulations to predict suppression of matter power spectrum for cosmological models with massive neutrinos. However, semi-analytical fits like halofit will themselves become obsolete with near-future surveys that promise to reach per cent level of precision. Moreover, implementing neutrinos as particles in numerical simulations is a topic of ongoing research, with results (see Brandbyge & Hannestad 2009; Viel, Haehnelt, & Springel 2010) contradictory at a level (factor of or higher) that can not be justified as due to (non)-inclusion of baryonic physics.

An alternative procedure to tackle small-scale non-linearities is to use higher order PT (e.g. Saito et al. 2008; Nishimichi et al. 2009; Saito, Takada, & Taruya 2009) to push further into the quasi-linear domain. Using high-resolution -body simulations as reference, *CarWhiPad09 have shown that although PT improves upon a linear description of the power spectrum at large scales (), it expectedly fails on smaller scales (). The range of scales where PT is reliable at the per cent level is both redshift and cosmological model dependent. For cosmologies close to WMAP best-fit parameters, Taruya et al. (2009) have shown that at redshift , the one-Loop standard perturbation sequence to the non-linear matter power spectrum is expected to converge with the -body simulation results to within per cent - only for scales . With the measurements from surveys expected to be at per cent level precision, these upcoming data sets create new challenges in analyses and need alternative ways to efficiently estimate cosmological parameters.

One might argue that a machine-learning approach to determine the non-linear response from varying parameter settings is a rather black-box approach that goes against the traditional approach to spectra: based on scientific understanding and physics. However, we view this direction as a pragmatic one: a new approach is urgent given the impending flood of new data from upcoming surveys, and in an age of supposed precision cosmology, we will be theory limited in this specific area. It is therefore crucial to strive towards per cent level precision in the determination of the non-linear power spectrum.

Machine-learning techniques – in the form of Gaussian processes – have already been used as cosmological non-linear emulators in Habib et al. (2007), Schneider et al. (2008), Heitmann et al. (2009), Lawrence et al. (2010) and *Schneider11. Gaussian process modeling (see MacKay 1997; Rasmussen & Williams 2006 for a basic introduction to Gaussian processes) is a non-linear interpolation scheme that, after optimal learning, is capable of making predictions when queried at a suitable input setting. There are several advantages and disadvantages when using neural networks and Gaussian processes to interpolate data. From a practical point of view, a neural network compresses data into a small number of weight parameters, so a large number of simulations could be fitted into a small number of files whereas a Gaussian process has to carry a large matrix which can be of the order of the number of points used for training the Gaussian process. Heitmann et al. (2009) dealt with large matrices by using principal component analysis (PCA) to reduce their sizes to ones easily manipulated. Again from a practical point of view, usually Gaussian processes can do better than neural networks in the case of a small number of training points given that a neural network could be flexible enough to be misused and misfit the data. From a theoretical point of view, the two methods should fare equally especially as there are certain kernels used in Gaussian processes which are equivalent to the interpolation and fit one would have with neural networks. Overall, given the implementation, we believe that the two methods should produce equivalent results especially if the artificial neural network (ANN) procedure is trained by a larger number of simulations.

While using any machine-learning technique in lieu of -body simulation output, it is critical that (i) the queried input setting not lie outside the input parameter ranges that are used during machine learning and (ii) the input parameter space must be sampled densely enough for the machine procedure to interpolate/predict accurately. Machine learning has been used in the fitting of cosmological functions (Auld et al. 2007; Fendt & Wandelt 2007; Auld, Bridges, & Hobson 2008) and photometric redshifts (Collister & Lahav, 2004).

halofit is accurate at the per cent level at best (see Heitmann et al. 2010). A far more accurate matter power spectrum calculator is the cosmic emulator (see Heitmann et al. 2009; Lawrence et al. 2010); although accurate at sub-per cent level, it makes predictions that are valid only for redshifts and does not include cosmological models with massive neutrinos. In order to extend the interpolation validity range to , as well as improve the accuracy levels, we work on a new technique to fit results from cosmological -body simulations using an ANN procedure with an improved Latin hypercube sampling of the cosmological parameter space. Using halofit spectra as mock -body spectra, we show that the ANN formalism enables a remarkable fit with a manageable number of simulations.

The outline of this paper is as follows. We discuss the concept of machine learning, in particular the details of neural networks, in Section 2. The improved Latin hypercube sampling of the underlying parameter space, which keeps the simulation number manageable and fitting accuracy high, is detailed in Section 3. Our fitting results are included in Section 4. Finally, we conclude in Section 5.

2 Machine Learning - The Neural Network

Machine learning is associated with a series of algorithms that allow a computational unit to evolve in its behavior, given access to empirical data. The major benefit is the potential to automatically learn complex patterns. As a subset of artificial intelligence, machine learning has been used in a variety of applications ranging from brain-machine interfaces to the analyses of stock market. There exist a range of techniques (see e.g. Nilsson 2005) including genetic algorithms, decision tree learning and gaussian processes (as referenced above). In this work we focus on the neural network technique.

An ANN is simply an interconnection of neurons or nodes analogous to the neural structure of the brain. This can take a more specific form whereby the nodes are arranged in a series of layers with each node in a layer connected, with a weighting, to all other nodes in adjacent layers. This is often referred to as a multi-layer perceptron (MLP). In this case one can impart values onto the nodes of the first layer (called the input layer), have a series of hidden layers and finally receive information from the last layer (called the output layer). The configuration of nodes is often called the network’s architecture and is specified from input to output as : : … : : . That is, a network with an architecture 7 : 49 : 50 has seven inputs, a single hidden layer with 49 nodes and finally 50 outputs. An extra node (called the bias node) is added to the input layer and connects to all the nodes in the hidden layers and the output layer (see Bishop 1995 for more details). The total number of weights for a generic architecture : : … : : can be calculated using the formula


where the summation index is over the hidden layers only. For a network with a single hidden layer, the second term on the right-hand side is absent. Specifically, the architecture 7 : 49 : 50 has a total of weights, collectively denoted by w.

Each node (except the input nodes) is a neuron with an , , taking as its argument


where the sum is over all nodes (of the previous layer) sending connections to the th node (of the current layer). The activation functions are typically taken to be sigmoid functions such as . The sigmoid functions impart some degree of non-linearity to the neural network models. A network becomes overly non-linear if the weights w deviate significantly from zero. This drives the activation of the nodes to saturation. The number and size of the hidden layers add to the complexity of ANNs. For a particular input vector, the output vector of the network is determined by progressing sequentially through the network layers, from inputs to outputs, calculating the activation of each node.

Adjusting the weights w to get the desired mapping is called the training of the network. For matter power spectrum estimation, we use a training set of -body simulations for which we have full information about the non-linear matter power spectra , as well as the underlying cosmological parameters: , where and are the present-day normalized Hubble parameter in units of 100 km s Mpc, the present-day matter and baryonic normalized energy densities, the primordial spectral index, the constant equation of state parameter for dark energy, the amplitude of fluctuation on an 8 Mpc scale and the total neutrino mass, respectively.

Given the training set, the network can be used to learn some parametrization to arbitrary accuracy by training the weights w. This is done by minimizing a suitable ,


with respect to the weights . The sum is over all the cosmologies in the training set. The sum is over all the nodes in the output layer. Note that each output node samples the matter power spectrum at some specific scale, . is the true non-linear matter power spectrum for the specific cosmology . In this paper, we use halofit’s approximation for . In Paper II -body simulations will be used to calculate . Given the weights w, is the ANN’s predicted power spectrum for the th cosmology. In our fitting procedure, we work with the ratio of the non-linear to linear power spectrum, namely , where is calculated using camb (Lewis, Challinor, & Lasenby, 2000). As such, weighting the numerator in Eq. 3 by gives,


The ratio is a flatter function and gives better performance, particularly at higher redshifts where the ratio tends to 1. Given the weights w, in Eq. 5 is the network’s prediction of the ratio for the specific cosmology . The predicted non-linear spectrum in Eq. 4 is recovered by multiplying by the corresponding linear spectrum .

We ran -body simulations over a range of cosmological parameters with the publicly available adaptive mesh refinement (AMR), grid-based hybrid (hydro+-body) code enzo111 (Norman et al. 2007; O’Shea et al. 2010). We include radiative cooling of baryons using an analytical approximation (Sarazin & White, 1987) for a fully ionized gas with a metallicity of . The cooling approximation is valid over the temperature range from K. Below K, the cooling rate is effectively zero. However, we do not account for metal-line cooling, supernova (SN) feedback or active galactic nucleus (AGN) feedback. It is worth mentioning here that van Daalen et al. (2011) have shown that the inclusion of AGN feedback can reproduce the optical and X-ray observations of groups of galaxies, and decrease the power relative to dark matter-only simulations at , ranging from per cent at to as much as per cent at . As such, understanding and including the effects of baryonic physics in numerical simulations will be critical to predicting the non-linear matter power spectrum at sub-per cent level. Further, the ANN prescription we are using in this paper could also be used for fitting these kinds of baryonic effects by introducing additional parameters beyond the cosmological ones, especially since gasdynamical runs are much more expensive than dark matter-only simulations.

Figure 1: Top: matter power spectrum evaluated at redshift (upper set) and (lower set). At each redshift, the spectrum is calculated using (i) linear theory (dot-dashed), (ii) -body (dots), (iii) halofit (dashed) and (iv) cosmic emulator (solid, see Lawrence et al. 2010). The cosmological parameters are: with . In this method paper, we use the linear matter power spectrum for scales . Due to the lack of force resolution on small scales, our -body power spectrum is accurate at per cent level for . Middle: The ratio of the -body spectrum to cosmic emulator’s prediction at . The error bars correspond to the scatter in the -body results. The horizontal dotted lines denote and per cent error. Bottom: The same as the middle panel, at .
Figure 2: Linear theory, halofit and -body spectra from Fig. 1, left-hand panel are re-plotted – with the only difference that on scales , the halofit and -body spectra are approximated by linear theory. The stitched spectra are then sampled at 50 -values between . In this paper, we use the halofit spectrum as for ANN training. In Paper, II the -body spectrum will be used along with the one-Loop standard PT.

In Fig. 1, top panel, we show the power spectrum for a cosmological model (0.13, 0.0224, 0.986, -1.23, 0.72, 0 eV), with . The spectrum is evaluated at redshift (upper set) and (lower set). At each redshift, the power spectrum is calculated using (i) linear theory (dot-dashed), (ii) -body (dots), (iii) halofit (dash) and (iv) cosmic emulator (solid). The vertical dash line at is the highest upto which our -body power spectrum is accurate at per cent level. We average over realizations of the initial power spectrum to suppress the scatter in the -body results. On smaller scales, our numerical simulations do not have the force resolution to give accurate results. The ratio of the -body spectrum to cosmic emulator’s prediction is shown in the middle () and bottom () panels. The error bars correspond to the scatter in the -body results.

In this method paper, we match the linear theory power spectrum to the non-linear power spectrum from simulations at . In Paper II we intend to use the one-Loop standard PT as implemented by Saito et al. (2008) for estimating the matter power spectrum upto and stitch it with the -body spectrum. The stitched spectrum will be sampled at 50 -values between . In Fig. 2, we show the stitched-and-sampled -body power spectrum (solid dots) which we will use as for ANN training in Paper II. The halofit spectrum, sampled at the same -values, is also shown (open circles) and is used as in this method paper.

In Eq. 5, optimizing the weights w in order to minimize generates an ANN that predicts the power spectrum very well for the specific cosmologies in the training set. However, such a network might not make accurate predictions for cosmologies included in the training set. This usually indicates (i) an overly simple network architecture (very few hidden layer nodes), (ii) very sparsely/poorly sampled parameter space and/or (iii) a highly complex non-linear mapping that actually over-fits to the noise on the training dataset. In order to generate smoother network mappings that generalize better when presented with new cosmologies that are not part of the training set, a penalty term is added to the cost function ,


where is the weight connecting the th node to the th node of the next layer. , usually a quadratic sum of the weights, prevents them from becoming too large during the training process, by penalizing in proportion to the sum. The hyperparameter controls the degree of regularization of the network’s non-linearity. After having initialized , its value itself is re-estimated during the training process iteratively. See Bishop (1995) for more details.

Thus, the overall cost function which is presented to the ANN for minimization with respect to the weights is,


To minimize the cost function w.r.t. the weights w, we use an iterative quasi-Newton algorithm which involves evaluating the inverse Hessian matrix


Specifically, we employ the Broyden–Fletcher–Goldfarb–Shanno (BFGS) method to approximate the inverse Hessian matrix (see Bishop 1995, for details). Starting with randomly assigned weights w, their values are re-estimated iteratively, assuring that each iteration proceeds in a direction that lowers the cost function . After each iteration to the weights, Eq. 7 is also calculated for what is known in neural network parlance as a validation set, in order to avoid over-fitting to the training set. The validation set for our application of neural networks, is a small set of simulations with known and . The final weights are chosen such as to give the best fit (minimum ) to the validation set. The network training is considered finished once is minimized with respect to the validation set. The trained network can now be used to predict for new cosmologies. In practice, a number of networks are trained that start with an alternative random configuration of weights. The trained networks are collectively called a committee of networks and subsequently give rise to better performance. The final output is usually given by averaging over the outputs of the committee members.

The ANN technique has been used successfully in empirical photometric redshift estimation with the annz (Collister & Lahav, 2004) package. annz learns an effective parameterization of redshift as a function of broad-band photometric colours by training on a representative set of galaxies that have both photometric and spectroscopic information. This has been shown to be more successful than template or synthetic-based methods (Abdalla et al. 2011; Thomas, Abdalla, & Lahav 2011). In this work we use an MLP similar to the original annz engine.

Our intention is to use this neural network technique to learn the non-linear matter power spectrum as a function of cosmological parameters by training on -body simulations. This natural fitting procedure removes both the effort and unnecessary potential bias that results from invoking ultimately imperfect sets of fitting equations such as the halofit. As can be seen in Fig. 2, the halofit predictions are in error by as much as per cent on small scales. As we will discuss in Section 4, the ANN technique is extremely fast and, more importantly, accurate.

3 Latin Hypercube Parameter Sampling

In order to fit a set of parameters optimally one strives to sample them as finely and as evenly as possible. However, a regularly spaced grid with sampling intervals along one dimension and parameters scales as . For a six-dimensional parameter space with only grid intervals, this quickly escalates to points. The problem is exacerbated because an -body simulation is a computationally expensive activity. To further compound this issue, each parameter configuration needs to be simulated over multiple realizations to beat down simulation (sample) variance. An alternate approach could be to interpolate the fitting function over a selection of randomly distributed points throughout the parameter space. However, this is prone to statistical clustering and will lead to a degradation of the machine-learning fit for parameters or regions affected by it. In order to circumvent these problems, we select parameters distributed on a Latin hypercube.

Figure 3: Left: An example of a Latin hypercube distribution. Every interval and is sampled; however each row and column are sampled only once. Right: an improved Latin hypercube where the distribution is more evenly spread through the space. Each subspace is equally sampled and there are no voids or clusters as in the left-hand panel (bottom left and right corners, respectively).

A square grid is said to be populated as a Latin square if, and only if, there is exactly one sample in each row and each column of the square. This is illustrated clearly in Fig. 3. A similar sampling scheme was developed first by Leonhard Euler who indexed the samples with Latin characters, motivating the name ‘Latin square’. A Latin hypercube is a generalization of Euler’s Latin square to a higher dimensional parameter space and is an example of a stratified sampling technique. This ensures that each and every segment/interval along a parameter axis is sampled with high resolution without a vast number of points. That is, one can sample a -dimensional space with simulations and have all parameters evaluated along every increment, where and are the upper and lower limits of the parameter, respectively. Therefore, it is independent of . However, a random implementation of a Latin hypercube algorithm can still lead to statistically under-sampled regions. An example of this can be seen in Fig. 3. Each panel shows a random implementation of Latin hypercube algorithm. In both panels, the square is partitioned into four subspaces. The left-hand panel has voids (and clusters) in two of its subspaces. The right-hand panel has each subspace equally sampled (while still obeying the Latin hypercube definition) and represents an improved Latin hypercube sampling. In this case the sample space is partitioned into equally probable subspaces and the variance in the pairwise separation of the sampled points is minimized.

Since the introduction of the Latin hypercube sampling technique (McKay et al., 1979), the procedure has become common in computer science, uncertainty analysis and engineering emulation (where simulation of complex machinery is overwhelmingly time consuming). Similarly, variations of the Latin hypercube sampling technique have been implemented in cosmological analyses before, e.g., Habib et al. (2007), Heitmann et al. (2009), Schneider et al. (2011) and references therein. In this paper, we use the improved Latin hypercube technique to set up the cosmological models to be used for ANN training.

3.1 Setting up an improved Latin hypercube for cosmological parameters

We varied six cosmological parameters between the limits specified in Table 1. The limits on this six-dimensional parameter space are chosen so as to include the WMAP 7-year+BAO+ (Komatsu et al., 2011) constraints (see Table 1).

Cosmological parameters Lower value Upper value WMAP 7-year+BAO+
0.110 0.165 0.1352 0.0036
0.021 0.024 0.02255 0.00054
0.85 1.05 0.968 0.012
-1.35 -0.65 -1.1 0.14
0.60 0.95 0.816 0.024
(eV) 0 1.1
Note. Komatsu et al. (2011); per cent CL for
Table 1: The six cosmological parameters and their ranges, used in generating the ANN training and validation sets. This six-dimensional parameter space is sampled using the improved Latin hypercube technique (see text for details). The last column shows the corresponding WMAP 7-year+BAO+ constraints at 68 per cent CL.

Throughout this paper, we only consider spatially flat models with the present-day CMB temperature . We also assume that all massive neutrino species are degenerate. The effective number of neutrino species is fixed at . We derive the Hubble parameter using the WMAP 7-year+BAO constraint on the acoustic scale , where is the distance to the last scattering surface and is the sound horizon at the redshift of last scattering. We derive as follows.

(i) For a given and , compute the redshift of the last scattering surface, , using the fit proposed by Hu & Sugiyama (1996):


(ii) For a given , and , choose a value for and compute its evolution, . Here we follow section 3.3 from Komatsu et al. (2011), which includes the effect of massive neutrinos on :



is the present-day neutrino temperature and is the present-day normalized photon energy density. Given , the function calculates the contribution of neutrinos to the radiation energy density at scale factor .

(iii) Using from step (ii), compute the comoving sound horizon at the last scattering redshift :


(iv) Using from step (iii), together with the WMAP 7-year+BAO constraint on the acoustic scale , compute the comoving distance to the last scattering surface, :


(v) Using from step (ii), compute the comoving distance to the surface of last scattering :


(vi) Compare results from steps (iv) and (v). Minimize the difference by varying in step (ii) and re-estimating steps (ii)-(v).

Using Table 1 as the parameter priors, we sampled this six-dimensional parameter space with an improved Latin hypercube technique. We generated 130 cosmologies to be used as the ANN training set and another 32 cosmologies for the validation set. We show the training set (upper triangle) and the validation set (lower triangle) in Fig. 4.

Figure 4: Upper triangle: ANN training set with 130 viable cosmologies, in a six-dimensional parameter space. Lower triangle: ANN validation set with 32 viable cosmologies, in a six-dimensional parameter space. See Table 1 for the parameter priors used to generate the training and validation sets.
Figure 5: Upper triangle: extending the ANN training set (upper triangle in Fig. 4) with 70 cosmologies with . Lower triangle: extending the ANN validation set (lower triangle in Fig. 4) with 18 cosmologies with .

As can be seen in Fig. 3, a major advantage of improved Latin hypercube sampling technique is the relatively uniform coverage it provides. This is, of course, highly useful for training a machine-learning algorithm. As with any interpolation mechanism, one hopes that the neural network can generalize from what it has learned to new and slightly different input data (in this case cosmological parameters). In reality, the response will be uncertain in poorly trained areas. Therefore, the caveat with our sampling will reside near the edges of the parameter hypercube. A parameter value that we might want emulated may not be encapsulated within the hypervolume of a simulated, and therefore trained, point. This can be understood with reference to Fig. 4. The performance of a neural network can severely degrade near the parameter boundaries. The solution is simply to choose prior ranges that are marginally wider than those of real interest. The allowance could easily be found empirically by projecting the hypercube realizations. The real problem in cosmology therefore arises when one has a parameter that is physically bounded, an example being the neutrino mass .

Adding several additional simulations at the parameter boundary may not be a computationally feasible solution to the problem due to the multi-dimensionality of the parameter space. Instead we propose to use a nested hypercube with dimensions. We fixed and varied the rest of the parameters over their aforementioned limits. We extended the ANN training and validation sets to include this five-dimensional hyperplane. Towards this, we generated 70 (for training) and 18 (for validation) cosmologies with . Fig. 5 shows the five-dimensional hyperplane.

4 Analysis and Results

We tested the precision with which our neural network can predict the non-linear matter power spectrum. We selected the combination 7 : : 50 as our network architecture, where (number of nodes in the hidden layer) was varied from 7 to 56, in steps of 7. The number of inputs were fixed at 7, corresponding to including redshift . Note that we do not sample the redshift in the Latin hypercube but instead evaluate , at 111 redshifts between and , using the existing prescription halofit coupled with the camb software. These halofit-generated spectra serve as mock -body spectra. We sampled at 50 points between . Since our training and validation sets have () and () cosmologies, respectively, we calculated for each cosmology, at 111 redshifts. These are scaled by their respective linear spectra , before being fed to the neural network. Thus, the overall size of the training set that we train our ANN with is . Likewise, we have patterns in the validation set. We trained a committee of 16 ANNs at each setting. The weights for each ANN were randomly initialized (the random configuration being different for each ANN). The weights are allowed to evolve until (see Eq. 7) is minimized with respect to the cosmologies in the validation set.

Figure 6: Left-hand column: percentage error at redshift , between the predicted non-linear power spectrum (using ANN) and the true underlying mock (using halofit) for training set cosmologies shown in the upper triangles of Figs. 4 and 5. The profile is continuous as the 50 output values have been spline interpolated over the functional range. The rows (from top to bottom) correspond to in increments of 7. The mean error over all 200 cosmologies is shown by a dashed line – an indicator of any bias in the ANN training scheme. Middle- column: the same as the left-hand column, at redshift . Right-hand column: the same as the left-hand column, at redshift .

In Fig. 6, we show the performance of the trained ANNs with varying units, when presented with each of the 200 cosmologies in the training set. Note that we average the predictions over all 16 ANN committee members. The rows correspond to (from top to bottom) in increments of 7. The columns (from left to right) correspond to . The mean error over all 200 cosmologies in the training set is shown by a dashed line in each panel, to get an idea about any systematics in our ANN training scheme. With , the ANN predictions at redshifts and , on scales, are within per cent of the halofit power spectra. Although we show results at and , we have checked that the predictions are per cent level for all . Predictions are at the per cent level even up to redshift for , after which the performance degrades to per cent. We have checked and confirmed that the worst-performing cosmologies correspond to the parameter settings in which at least four of the six cosmological parameters are at their boundary values.

As we mentioned earlier, this fitting procedure will be less accurate near the boundaries of the parameter ranges because some parameter configurations may not be encapsulated within the volume of a training point. This also explains why the ANN performance is better at – the mid-point of the redshift range. Looking at the bias (dashed line in Fig. 6), we see that the distribution of errors in the ANN predictions is centered on zero, indicating that our interpolations are not biased. A negligible bias, and the fact that for cosmological setting the parameter priors (see Table 1) the non-linear power spectrum at is correctly predicted within per cent up to , demonstrates the stability of our ANN strategy. This marks a remarkable improvement over the currently popular interpolation scheme – the cosmic emulator, which has a significant number ( per cent) of cosmological models with errors at per cent level. We note, however, that the cosmic emulator, based on Gaussian processes, is able to achieve sub-per cent accuracy with only 37 distinct cosmologies while in the ANN scheme we use a suite of around 200 cosmologies. Comparing Fig. 9 from Heitmann et al. (2009) with our Fig. 6, we see that the ANN implementation performs better on all scales and redshifts.

In order to check the performance of our trained ANNs over parameter configurations that do not touch the Latin hypercube, we generated a testing set of 330 cosmologies (of which 150 have ). See Table 2 for the parameter limits of the testing set.

Cosmological parameters Lower value Upper value WMAP 7-year+BAO+
0.120 0.150 0.1352 0.0036
0.022 0.023 0.02255 0.00054
0.90 1.00 0.968 0.012
-1.15 -0.85 -1.1 0.14
0.70 0.85 0.816 0.024
(eV) 0 0.50
Note. Komatsu et al. (2011); per cent CL for
Table 2: The six cosmological parameters and their ranges, used in generating the ANN testing set. This six-dimensional parameter space is sampled using the improved Latin hypercube technique (see the text for details). The parameter ranges are chosen so as to avoid the boundaries of the parameter space. See Table 1 for the parameter boundaries. Note that the lower bound on neutrino mass is still set at zero, since neutrinos are physically bound (). The last column shows the WMAP 7-year+BAO+ constraints at 68 per cent CL.
Figure 7: Upper triangle: ANN testing set with 180 cosmologies with . Lower triangle: extending the ANN testing set with 150 cosmologies with . See Table 2 for the parameter priors used to generate the testing set.

A testing set serves another crucial purpose. Increasing the number of nodes in the hidden layer increases the flexibility of a neural network. An increasingly complex network can make extremely accurate predictions on the training set. This is evident from Fig. 6, where the prediction over the training set becomes progressively better (from top to bottom) with increasing units. However, such complex networks can adversely affect their generalizing ability when presented with a new dataset. Measuring the performance of a neural network on an independent dataset (known as a testing set) as a function of units helps in controlling its complexity. We show the testing set in Fig. 7, with the lower triangle corresponding to the 150 cosmologies with .

The performance of the trained ANNs as a function of units, over the cosmologies in the testing set, is shown in Fig. 8. Increasing from 42 to 49 reduces the error marginally. We have checked that increasing beyond 49 does not contribute to any further error reduction on the testing set, indicating that saturates the generalizing ability of the network. With , the ANN prediction for every cosmology, on scales, at redshifts , is within per cent of the halofit power spectra. For , in the high range, the performance degrades slightly to per cent. As expected, the ANN performs exceedingly well away from the boundaries of the parameter ranges. It is quite remarkable that our ANN scheme is capable of making predictions at sub-per cent level, especially on the testing set that is a part of the ANN training process.

Figure 8: The same as Fig. 6, using 330 testing set cosmologies shown in Fig. 7.

In Fig. 9, left panel, we show the non-linear power spectrum at redshifts and predicted by our trained ANN (solid dots) for the same cosmology that was used to generate the halofit power spectrum of Fig. 2, namely (0.13, 0.0224, 0.986, -1.23, 0.72, 0 eV) with . is fixed at 49, as discussed above. For comparison, the halofit power spectra are re-plotted as starred points. The prediction errors at and are shown in the middle and the right-hand panels, respectively. The ANN predictions are well within per cent over the scales of interest.

We reiterate that this method for reconstructing the non-linear power spectrum will only function for the parameters and ranges that have been simulated and trained with PkANN. The intention of this study is to provide a technique for high precision fits in the concordance model for the oncoming generation of surveys. This should therefore act as a safety mechanism as it demonstrates that the range of validity has been breached, as often occurs with blind application of other fits. In Paper II, we will put our ANN formalism to test using matter power spectra calculated from -body simulations. On mildly non-linear scales () the power spectrum from -body simulations is expected to vary smoothly as a function of cosmological parameters (Heitmann et al., 2006, 2009; Agarwal & Feldman, 2011). As such, our ANN interpolation scheme that we have shown here to work well on halofit power spectra is expected to perform satisfactorily on -body spectra as well.

Figure 9: Left: ANN predictions (solid dots) at (upper) and (lower) for the halofit power spectra of Fig. 2 (re-plotted here as starred points – difficult to see due to the excellent ann predictions). The cosmological model is (0.13, 0.0224, 0.986, -1.23, 0.72, 0 eV) with . is fixed at 49, as discussed in the text. The ratio of ann predictions to the halofit spectra is shown in the other panels. Middle: percentage error in ANN predictions, at . Right: percentage error in ANN predictions, at . The ANN predictions are within per cent over the scales of interest.

5 Conclusions

The advent of the era of precision cosmology poses an immense challenge to theoretical physics. The upcoming generation of surveys has the potential to breach per cent level of accuracy. Such high-precision data will improve our constraints on cosmological parameters including dark energy and neutrino masses. Efficiently dealing with this impending flood of precise data on ever smaller scales and lower redshifts requires that we move on from linear theory as well as any imperfect sets of fitting equations. Although numerical simulations are capable of achieving the levels of precision required by the near-future surveys, the high dimensionality of the cosmological parameter space renders their brute force usage impractical.

We have introduced a unique approach to coping with non-linearities in the matter power spectra in cosmology. By employing a multi-layer perceptron neural network together with an improved Latin hypercube parameter sampling technique, we have demonstrated that the non-linear spectrum can be reconstructed from a full set of cold dark matter parameters to better than per cent over the parameter space spanning confidence level around the WMAP 7-year central values. Parameters that are likely to reside by some hard physical prior, such as the neutrino mass, can be successfully brought under the realms of ANNs by sprinkling of extra simulations in the corresponding (e.g. ) hyper-plane.

Looking forward, our ANN procedure can be readily employed for a variety of cosmological tasks such as fitting halo mass functions obtained through high resolution -body simulations. Moreover, mixed datasets such as the matter power spectra and the halo mass functions can be combined and presented to a neural network as the training set. An ANN trained with such a heterogeneous dataset would be capable of cosmological parameter estimation when presented with the combined observations of the matter power spectrum and the measured halo mass function. The implementation of our technique avoids complex calculations and, through the execution of only the neural network weights, is extremely fast (predictions take less than a second). An automated PkANN function will be released with our program of -body results in a companion paper shortly. Beyond this we hope that with our method a collaborative effort could reduce non-linear error to only uncertainty in the -body simulation’s baryon interactions.

6 Acknowledgements

We thank Salman Habib, Katrin Heitmann and Joop Schaye for their insightful comments and suggestions on the present paper. This work is supported by the National Science Foundation through TeraGrid resources provided by the NCSA. SAT acknowledges UCL’s Institute of Origins for a Post-doctoral Fellowship. FBA and OL acknowledge the support of the Royal Society via a Royal Society URF and a Royal Society Wolfson Research Merit Award, respectively.


  • Abdalla et al. (2011) Abdalla F. B., Banerji M., Lahav O., Rashkov V., 2011, MNRAS, 417, 1891
  • Agarwal & Feldman (2011) Agarwal S., Feldman H. A., 2011, MNRAS, 410, 1647
  • Amendola et al. (2008) Amendola L., Kunz M., Sapone D., 2008, JCAP, 4, 13
  • Auld et al. (2008) Auld T., Bridges M., Hobson M. P., 2008, MNRAS, 387, 1575
  • Auld et al. (2007) Auld T., Bridges M., Hobson M. P., Gull S. F., 2007, MNRAS, 376, L11
  • Bird et al. (2012) Bird S., Viel M., Haehnelt M. G., 2012, MNRAS, 2175
  • Bishop (1995) Bishop C. M., 1995, Neural Networks for Pattern Recognition. Oxford Univ. Press, New York
  • Brandbyge & Hannestad (2009) Brandbyge J., Hannestad S., 2009, JCAP, 5, 2
  • Carlson et al. (2009) Carlson J., White M., Padmanabhan N., 2009, PRD, 80, 043531
  • Collister & Lahav (2004) Collister A. A., Lahav O., 2004, PASP, 116, 345
  • Eisenstein et al. (2011) Eisenstein D. J., Weinberg D. H., Agol E., Aihara H., Allende Prieto C., Anderson S. F., Arns J. A., Aubourg É., Bailey S., Balbinot E., et al., 2011, AJ, 142, 72
  • Fendt & Wandelt (2007) Fendt W. A., Wandelt B. D., 2007, ApJ, 654, 2
  • Habib et al. (2007) Habib S., Heitmann K., Higdon D., Nakhleh C., Williams B., 2007, PRD, 76, 083503
  • Heitmann et al. (2006) Heitmann K., Higdon D., Nakhleh C., Habib S., 2006, ApJL, 646, L1
  • Heitmann et al. (2009) Heitmann K., Higdon D., White M., Habib S., Williams B. J., Lawrence E., Wagner C., 2009, ApJ, 705, 156
  • Heitmann et al. (2010) Heitmann K., White M., Wagner C., Habib S., Higdon D., 2010, ApJ, 715, 104
  • Hu & Sugiyama (1996) Hu W., Sugiyama N., 1996, ApJ, 471, 542
  • Ivezic et al. (2008) Ivezic Z., Tyson J. A., Acosta E., Allsman R., Anderson S. F., Andrew J., Angel R., Axelrod T., Barr J. D., Becker A. C., Becla J., Beldica C., Blandford R. D., Bloom J. S., Borne K., Brandt W. N., Brown M. E., Bullock J. S., Burke D. L., Chandrasekharan S., Chesley S., Claver C. F., Connolly A., Cook K. H., Cooray A., Covey K. R., Cribbs C., Cutri R., Daues G., Delgado F., Ferguson H., Gawiser E., Geary J. C., Gee P., Geha M., Gibson R. R., Gilmore D. K., Gressler W. J., Hogan C., Huffer M. E., Jacoby S. H., Jain B., Jernigan J. G., Jones R. L., Juric M., Kahn S. M., Kalirai J. S., Kantor J. P., Kessler R., Kirkby D., Knox L., Krabbendam V. L., Krughoff S., Kulkarni S., Lambert R., Levine D., Liang M., Lim K., Lupton R. H., Marshall P., Marshall S., May M., Miller M., Mills D. J., Monet D. G., Neill D. R., Nordby M., O’Connor P., Oliver J., Olivier S. S., Olsen K., Owen R. E., Peterson J. R., Petry C. E., Pierfederici F., Pietrowicz S., Pike R., Pinto P. A., Plante R., Radeka V., Rasmussen A., Ridgway S. T., Rosing W., Saha A., Schalk T. L., Schindler R. H., Schneider D. P., Schumacher G., Sebag J., Seppala L. G., Shipsey I., Silvestri N., Smith J. A., Smith R. C., Strauss M. A., Stubbs C. W., Sweeney D., Szalay A., Thaler J. J., Vanden Berk D., Walkowicz L., Warner M., Willman B., Wittman D., Wolff S. C., Wood-Vasey W. M., Yoachim P., Zhan H., for the LSST Collaboration, 2008, arXiv:0805.2366v2
  • Komatsu et al. (2011) Komatsu E., Smith K. M., Dunkley J., Bennett C. L., Gold B., Hinshaw G., Jarosik N., Larson D., Nolta M. R., Page L., Spergel D. N., Halpern M., Hill R. S., Kogut A., Limon M., Meyer S. S., Odegard N., Tucker G. S., Weiland J. L., Wollack E., Wright E. L., 2011, ApJS, 192, 18
  • Lawrence et al. (2010) Lawrence E., Heitmann K., White M., Higdon D., Wagner C., Habib S., Williams B., 2010, ApJ, 713, 1322
  • Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473
  • MacKay (1997) MacKay D. J. C., 1997, Gaussian processes – a replacement for supervised neural networks? Lecture notes for a tutorial at NIPS 1997
  • McKay et al. (1979) McKay M. D., Beckman R. J., Conover W. J., 1979, Technometrics, 21, 239
  • Nilsson (2005) Nilsson N. J., 2005, Machine learn., 56, 387
  • Nishimichi et al. (2009) Nishimichi T., Shirata A., Taruya A., Yahata K., Saito S., Suto Y., Takahashi R., Yoshida N., Matsubara T., Sugiyama N., Kayo I., Jing Y., Yoshikawa K., 2009, PASJ, 61, 321
  • Norman et al. (2007) Norman M. L., Bryan G. L., Harkness R., Bordner J., Reynolds D., O’Shea B., Wagner R., 2007, arXiv:0705.1556
  • O’Shea et al. (2010) O’Shea B. W., Bryan G., Bordner J., Norman M. L., Abel T., Harkness R., Kritsuk A., 2010, Astrophysics Source Code Library, 10072
  • Rasmussen & Williams (2006) Rasmussen C. E., Williams C., 2006, Gaussian Processes for Machine Learning. MIT Press, Cambridge, MA
  • Saito et al. (2008) Saito S., Takada M., Taruya A., 2008, Physical Review Letters, 100, 191301
  • Saito et al. (2009) —, 2009, PRD, 80, 083528
  • Sarazin & White (1987) Sarazin C. L., White III R. E., 1987, ApJ, 320, 32
  • Schneider et al. (2011) Schneider M. D., Holm Ó., Knox L., 2011, ApJ, 728, 137
  • Schneider et al. (2008) Schneider M. D., Knox L., Habib S., Heitmann K., Higdon D., Nakhleh C., 2008, PRD, 78, 063529
  • Scoccimarro et al. (1999) Scoccimarro R., Zaldarriaga M., Hui L., 1999, ApJ, 527, 1
  • Smith et al. (2003) Smith R. E., Peacock J. A., Jenkins A., White S. D. M., Frenk C. S., Pearce F. R., Thomas P. A., Efstathiou G., Couchman H. M. P., 2003, MNRAS, 341, 1311
  • Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  • Taruya et al. (2009) Taruya A., Nishimichi T., Saito S., Hiramatsu T., 2009, PRD, 80, 123503
  • The Dark Energy Survey Collaboration (2005) The Dark Energy Survey Collaboration, 2005, arXiv:astro-ph/0510346
  • Thomas et al. (2010) Thomas S. A., Abdalla F. B., Lahav O., 2010, Physical Review Letters, 105, 031301
  • Thomas et al. (2011) —, 2011, MNRAS, 412, 1669
  • van Daalen et al. (2011) van Daalen M. P., Schaye J., Booth C. M., Dalla Vecchia C., 2011, MNRAS, 415, 3649
  • Viel et al. (2010) Viel M., Haehnelt M. G., Springel V., 2010, Journal of Cosmology and Astro-Particle Physics, 6, 15
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description