Learning to predict the cosmological structure formation.

Matter evolved under the influence of gravity from minuscule density fluctuations. Nonperturbative structure formed hierarchically over all scales and developed non-Gaussian features in the Universe, known as the cosmic web. To fully understand the structure formation of the Universe is one of the holy grails of modern astrophysics. Astrophysicists survey large volumes of the Universe and use a large ensemble of computer simulations to compare with the observed data to extract the full information of our own Universe. However, to evolve billions of particles over billions of years, even with the simplest physics, is a daunting task. We build a deep neural network, the Deep Density Displacement Model ([Formula: see text]), which learns from a set of prerun numerical simulations, to predict the nonlinear large-scale structure of the Universe with the Zel'dovich Approximation (ZA), an analytical approximation based on perturbation theory, as the input. Our extensive analysis demonstrates that [Formula: see text] outperforms the second-order perturbation theory (2LPT), the commonly used fast-approximate simulation method, in predicting cosmic structure in the nonlinear regime. We also show that [Formula: see text] is able to accurately extrapolate far beyond its training data and predict structure formation for significantly different cosmological parameters. Our study proves that deep learning is a practical and accurate alternative to approximate 3D simulations of the gravitational structure formation of the Universe.

Astrophysicists require a large amount of simulations to extract the information from observations (1–8). At its core, modeling structure formation of the Universe is a computationally challenging task; it involves evolving billions of particles with the correct physical model over a large volume over billions of years (9–11). To simplify this task, we either simulate a large volume with simpler physics or a smaller volume with more complex physics. To produce the cosmic web (12) in large volume, we select gravity, the most important component of the theory, to simulate at large scales. A gravity-only -body simulation is the most popular and effective numerical method to predict the full 6D phase-space distribution of a large number of massive particles whose position and velocity evolve over time in the Universe (13). Nonetheless, -body simulations are relatively computationally expensive, thus making the comparison of the -body–simulated large-scale structure (of different underlying cosmological parameters) with the observed Universe a challenging task. We propose to use a deep model that predicts the structure formation as an alternative to -body simulations.

Deep learning (14) is a fast-growing branch of machine learning, where recent advances have led to models that reach and sometimes exceed human performance across diverse areas, from analysis and synthesis of images (15–17), sound (18, 19), text (20, 21), and videos (22, 23) to complex control and planning tasks as they appear in robotics and game play (24–26). This new paradigm is also significantly impacting a variety of domains in the sciences, from biology (27, 28) to chemistry (29, 30) and physics (31, 32). In particular, in astronomy and cosmology, a growing number of recent studies are using deep learning for a variety of tasks, ranging from analysis of cosmic microwave background (33–35), large-scale structure (36, 37), and gravitational lensing effects (38, 39) to classification of different light sources (40–42).

The ability of these models to learn complex functions has motivated many to use them to understand the physics of interacting objects, leveraging image, video, and relational data (43–53). However, modeling the dynamics of billions of particles in N-body simulations poses a distinct challenge.

In this work, we show that a variation on the architecture of a well-known deep-learning model (54) can efficiently transform the first-order approximations of the displacement field and approximate the exact solutions, thereby producing accurate estimates of the large-scale structure. Our key objective is to prove that this approach is an accurate and computationally efficient alternative to expensive cosmological simulations, and, to this end, we provide an extensive analysis of the results in the following section.

The outcome of a typical N-body simulation depends on both the initial conditions and on cosmological parameters which affect the evolution equations. A striking discovery is that the Deep Density Displacement Model (), trained by using a single set of cosmological parameters, generalizes to new sets of significantly different parameters, minimizing the need for training data on a diverse range of cosmological parameters.

Setup

We build a deep neural network, , with similar input and output of an -body simulation. The input to our is the displacement field from the Zel’dovich Approximation (ZA) (55). A displacement vector is the difference of a particle position at target redshift —i.e., the present time—and its Lagrangian position on a uniform grid. ZA evolves the particles on linear trajectories along their initial displacements. It is accurate when the displacement is small; therefore, ZA is frequently used to construct the initial conditions of -body simulations (56). As for the ground truth, the target displacement field is produced by using FastPM (57), a recent approximate N-body–simulation scheme that is based on a particle-mesh (PM) solver. FastPM quickly approaches a full N-body simulation with high accuracy and provides a viable alternative to direct N-body simulations for the purpose of our study.

A significantly faster approximation of N-body simulations is produced by second-order Lagrangian perturbation theory (2LPT), which bends each particle’s trajectory with a quadratic correction (58). The 2LPT is used in many cosmological analyses to generate a large number of cosmological simulations for comparison of the astronomical dataset against the physical model (59, 60) or to compute the covariance of the dataset (61–63). We regard 2LPT as an effective way to efficiently generate a relatively accurate description of the large-scale structure, and therefore we select 2LPT as the reference model for comparison with .

We generate 10,000 pairs of ZAs as input and accurate FastPM approximations as the target. We use simulations of -body particles in a volume of (600 million light years, where is the Hubble parameter). The particles have a mean separation of per dimension.

An important choice in our approach is training with a displacement field rather than a density field. Displacement field and density field are two ways of describing the same distribution of particles. And an equivalent way to describe a density field is the overdensity field, defined as , with denoting the mean density. The displacement field and overdensity field are related by Eq. .When the displacement field is small and has zero curl, the choice of overdensity vs. displacement field for the output of the model is irrelevant, as there is a bijective map between these two representations, described by the equation:However, as the displacements grow into the nonlinear regime of structure formation, different displacement fields can produce identical density fields (e.g., ref. 64). Therefore, providing the model with the target displacement field during the training eliminates the ambiguity associated with the density field. Our inability to produce comparable results when using the density field as our input and target attests that relevant information resides in the displacement field ().

Results and Analysis

Fig. 1 shows the displacement vector field as predicted by (Left) and the associated point-cloud representation of the structure formation (Right). It is possible to identify structures such as clusters, filaments, and voids in this point-cloud representation. We proceed to compare the accuracy of and 2LPT compared with ground truth.

Fig. 1.

The displacement vector field (Left) and the resulting density field (Right) produced by . The vectors in Left are uniformly scaled down for better visualization.

Point-Wise Comparison.

Let denote the displacement field, where is the number of spatial-resolution elements in each dimension (). A natural measure of error is the relative error , where is the true displacement field (FastPM), and is the prediction from 2LPT or . Fig. 2 compares this error for different approximations in a 2D slice of a single simulation. We observe that predictions are very close to the ground truth, with a maximum relative error of 1.10 over all 1,000 simulations. For 2LPT, this number is significantly higher at 4.23. In average, the result of comes with a 2.8% relative error, while for 2LPT, it equals 9.3%.

Fig. 2.

The columns show 2D slices of full-particle distribution (Upper) and displacement vector (Lower) by various models: FastPM, the target ground truth, a recent approximate N-body simulation scheme that is based on a PM solver (A); ZA, a simple linear model that evolves particle along the initial velocity vector (B); 2LPT, a commonly used analytical approximation (C); and deep-learning model () as presented in this work (D). While FastPM (A) served as our ground truth, B–D include color for the points or vectors. The color indicates the relative difference between the target location (A) or displacement vector and predicted distributions by various methods (B–D). The error bar shows that denser regions have a higher error for all methods, which suggests that it is harder to predict highly nonlinear region correctly for all models: , 2LPT, and ZA. Our model has the smallest differences between predictions and ground truth among the above models B–D.

Two-Point Correlation Comparison.

As suggested by Fig. 2, the denser regions seem to have a higher error for all methods—that is, more nonlinearity in structure formation creates larger errors for both and 2LPT. The dependence of error on scale is computed with two- and three-point correlation analysis.

Cosmologists often use compressed summary statistics of the density field in their studies. The most widely used of these statistics are the two-point correlation function (2PCF) and its Fourier transform, the power spectrum :where the ensemble average is taken over all possible realizations of the Universe. Our Universe is observed to be both homogeneous and isotropic on large scales—i.e., without any special location or direction. This allows one to drop the dependencies on and on the direction of , leaving only the amplitude in the final definition of . In the second equation, is simply the Fourier transform of and captures the dispersion of the plane-wave amplitudes at different scales in the Fourier space. is the 3D wavevector of the plane wave, and its amplitude (the wavenumber) is related to the wavelength by . Due to isotropy of the Universe, we drop the vector form of and .

Because FastPM, 2LPT, and take the displacement field as input and output, we also study the two-point statistics for the displacement field. The displacement power spectrum is defined as:

We focus on the Fourier-space representation of the two-point correlation. Because the matter and the displacement power spectrum take the same form, in what follows, we drop the subscript for matter and displacement field and use to stand for both matter and displacement power spectrum. We use the transfer function and the correlation coefficient as metrics to quantify the model performance against the ground truth (FastPM) in the two-point correlation. We define the transfer function as the square root of the ratio of two power spectra,where is the density or displacement power spectrum as predicted by 2LPT or , and, analogously, is the ground truth predicted by FastPM. The correlation coefficient r(k) is a form of normalized cross-power spectrum,where is the cross-power spectrum between 2LPT or predictions and the ground-truth (FastPM) simulation result. The transfer function captures the discrepancy between amplitudes, while the correlation coefficient can indicate the discrepancy between phases as functions of scales. For a perfectly accurate prediction, and are both 1. In particular, describes stochasticity, the fraction of the variance in the prediction that cannot be explained by the true model.

Fig. 3 shows the average power spectrum, transfer function , and stochasticity of the displacement field and the density field over 1,000 simulations. The transfer function of density from 2LPT predictions is 2% smaller than that of FastPM on large scales (). This is expected since 2LPT performs accurately on very large scales (). The displacement transfer function of 2LPT increases >1 at and then drops sharply. The increase of the 2LPT displacement transfer function is because 2LPT overestimates the displacement power at small scales (e.g., ref. 65). There is a sharp drop of power near the voxel scale because smoothing over voxel scales in our predictions automatically erases power at scales smaller than the voxel size.

Fig. 3.

(A) Displacement and density power-spectrum of FastPM (orange), 2LPT (blue), and c (green) (Top); transfer function—i.e., the square root of the ratio of the predicted power-spectrum to the ground truth (Middle); and 1 – , where is the correlation coefficient between the predicted fields and the true fields (Bottom). Results are the averaged values of 1,000 test simulations. The transfer function and correlation coefficient of the predictions are nearly perfect from large to intermediate scales and outperform our benchmark 2LPT significantly. (B) The ratios of the multipole coefficients () (to the target) of the two 3PCFs for several triangle configurations. The results are averaged over 10 test simulations. The error bars (padded regions) are the SDs derived from 10 test simulations. The ratio shows that the 3PCF of is closer than 2LPT to our target FastPM with lower variance.

Now, we turn to the predictions: Both the density and displacement transfer functions of the differ from 1 by a mere 0.4% at scale , and this discrepancy only increases to 2% and 4% for density field and displacement field, respectively, as increases to the Nyquist frequency at ?. The stochasticity hovers at ∼ and for most scales. In other words, for both the density and displacement fields, the correlation coefficient between the predictions and FastPM simulations, all the way down to small scales of , is >90%. The transfer function and correlation coefficient of the predictions show that it can reproduce the structure formation of the Universe from large to seminonlinear scales. significantly outperforms our benchmark model 2LPT in the two-point function analysis. only starts to deviate from the ground truth at fairly small scales. This is not surprising, as the deeply nonlinear evolution at these scales is more difficult to simulate accurately and appears to be intractable by current analytical theories (66).

Three-Point Correlation Comparison.

The three-point correlation function (3PCF) expresses the correlation of the field of interest among three locations in the configuration space, which is equivalently defined as bispectrum in Fourier space. Here, we concentrate on the 3PCF for computational convenience:where = and = . Translation invariance guarantees that is independent of . Rotational symmetry further eliminates all direction dependence except dependence on , the angle between and . The multipole moments of , , where is the Legendre polynomial of degree , can be efficiently estimated with pair counting (67). While the input (computed by ZA) do not contain significant correlations beyond the second order (power spectrum level), we expect to generate densities with a 3PCF that mimics that of ground truth.

We compare the 3PCF calculated from FastPM, 2LPT, and by analyzing the 3PCF through its multipole moments . Fig. 3 shows the ratio of the binned multipole coefficients of the two 3PCFs for several triangle configurations, , where can be the 3PCF for M or 2LPT and is the 3PCF for FastPM. We used 10 radial bins with . The results are averaged over 10 test simulations, and the error bars are the SD. The ratio shows that the 3PCF of M is closer to FastPM than 2LPT, with smaller error bars. To further quantify our comparison, we calculate the relative 3PCF residual defined bywhere is the number of () bins. The mean relative 3PCF residual of the and 2LPT predictions compared with FastPM are and , respectively. The accuracy on 3PCF is also an order of magnitude better than 2LPT, which indicates that the is far better at capturing the non-Gaussian structure formation.

Generalizing to New Cosmological Parameters

So far, we train our model using a “single” choice of cosmological parameters (hereafter ) and (68). is the primordial amplitude of the scalar perturbation from cosmic inflation, and is the fraction of the total energy density that is matter at the present time, and we will call it matter density parameter for short. The true exact value of these parameters is unknown, and different choices of these parameters change the large-scale structure of the Universe; Fig. 4.

Fig. 4.

We show the differences of particle distributions and displacement fields when we change the cosmological parameters and . (A) The error bar shows the difference of particle distribution (Upper) and displacement fields (Lower) between and the two extremes for A = 0.2 A0 (Center) and (Right). (B) A similar comparison showing the difference of the particle distributions (Upper) and displacement fields (Lower) for smaller and larger values of with regard to , which was used for training. While the difference for the smaller value of () is larger, the displacement for the larger () is more nonlinear. This nonlinearity is due to concentration of mass and makes the prediction more difficult.

Here, we report an interesting observation: The trained on a single set of parameters in conjunction with ZA (which depends on and ) as input can predict the structure formation for widely different choices of and . From a computational point of view, this suggests a possibility of producing simulations for a diverse range of parameters, with minimal training data.

Varying Primordial Amplitude of Scalar Perturbations .

After training the using , we change in the input of our test set by nearly one order of magnitude: and . Again, we use 1,000 simulations for analysis of each test case. The average relative displacement error of remains < per voxel (compared with when train and test data have the same parameters). This is still well below the error for 2LPT, which has relative errors of and for larger and smaller values of , respectively.

Fig. 5 shows the transfer function and correlation coefficient for both and 2LPT. The performs much better than 2LPT for . For small , 2LPT does a better job than predicting the density transfer function and correlation coefficient at the largest scales; otherwise, predictions are more accurate than 2LPT at scales larger than . We observe a similar trend with 3PCF analysis: The 3PCFs of predictions are notably better than 2LPT ones for larger , compared with smaller , where it is only slightly better. These results confirm our expectation that increasing increases the nonlinearity of the structure-formation process. While 2LPT can predict fairly well in linear regimes, compared with , its performance deteriorates with increased nonlinearity. It is interesting to note that the prediction maintains its advantage, despite being trained on data from more linear regimes.

Fig. 5.

Similar plots as in Fig. 3, except that we test the two-point statistics when we vary the cosmological parameters without changing the training set (which has different cosmological parameters) or the trained model. We show predictions from and 2LPT when tested on different (A) and (B). We show the transfer function—i.e., the square root of the ratio of the predicted power spectrum to the ground truth (Upper)—and 1 – , where is the correlation coefficient between the predicted fields and the true fields (Lower). The prediction outperforms 2LPT prediction at all scales except in the largest scales, as the perturbation theory works well in linear regime (large scales).

Varying Matter Density Parameter .

We repeat the same experiments, this time changing to 0.5 and 0.1, while the model is trained on , which is quite far from both of the test sets. For , the relative residual displacement errors of the and 2LPT averaged over 1,000 simulations are 3.8% and 15.2%, and for , they are 2.5% and 4.3%. Fig. 5 show the two-point statistics for density field predicted by using different values of . For , the results show that the outperforms 2LPT at all scales, while for smaller , outperforms 2LPT on smaller scales (). As for the 3PCF of simulations with different values of , the mean relative 3PCF residual of the for and are 1.7% and 1.2%, respectively, and for 2LPT, they are 7.6% and 1.7%, respectively. The prediction performs better at than . This is again because the Universe is much more nonlinear at than . The learns more nonlinearity than is encoded in the formalism of 2LPT.

Conclusions

To summarize, our deep model can accurately predict the large-scale structure of the Universe as represented by FastPM simulations, at all scales, as seen in Table 1. Furthermore, learns to predict cosmic structure in the nonlinear regime more accurately than our benchmark model 2LPT. Finally, our model generalizes well to test simulations with cosmological parameters ( and ) significantly different from the training set. This suggests that our deep-learning model can potentially be deployed for a range of simulations beyond the parameter space covered by the training data (Table 1). Our results demonstrate that the successfully learns the nonlinear mapping from first-order perturbation theory to FastPM simulation beyond what higher-order perturbation theories currently achieve.

Table 1.

A summary of our analysis

Data	Point-wise	T(k)k = 0.11	r(k)k=0.11	T(k)k=0.50	r(k)k=0.50	3PCF
Test phase
2LPT density	N/A	0.96	1.00	0.74	0.94	0.0782
D^3M density	N/A	1.00	1.00	0.99	1.00	0.0079
2LPT displacement	0.093	0.96	1.00	1.04	0.90	N/A
D^3M displacement	0.028	1.00	1.00	0.99	1.00	N/A
A_s=1.8A_0
2LPT density	N/A	0.93	1.00	0.49	0.78	0.243
D^3M density	N/A	1.00	1.00	0.98	1.00	0.039
2LPT displacement	0.155	0.97	1.00	1.07	0.73	N/A
D^3M displacement	0.039	1.00	1.00	0.97	0.99	N/A
A_s=0.2A_0
2LPT density	N/A	0.99	1.00	0.98	0.99	0.024
D^3M density	N/A	1.00	1.00	1.03	1.00	0.022
2LPT displacement	0.063	0.99	1.00	0.95	0.98	N/A
D^3M displacement	0.036	1.00	1.00	1.01	1.00	N/A
Ω_m=0.5
2LPT density	N/A	0.94	1.00	0.58	0.87	0.076
D^3M density	N/A	1.00	1.00	1.00	1.00	0.017
2LPT displacement	0.152	0.97	1.00	1.10	0.80	N/A
D^3M displacement	0.038	1.00	1.00	0.98	0.99	N/A
Ω_m=0.1
2LPT density	N/A	0.97	1.00	0.96	0.99	0.017
D^3M density	N/A	0.99	1.00	1.04	1.00	0.012
2LPT displacement	0.043	0.97	1.00	0.97	0.98	N/A
D^3M displacement	0.025	0.99	1.00	1.02	1.00	N/A

The unit of k is hMpc−1. N/A, not applicable.

Looking forward, we expect that replacing FastPM with exact N-body simulations would improve the performance of our method. As the complexity of our model is linear in the number of voxels, we expect to be able to further improve our results if we replace the FastPM simulations with higher-resolution simulations. Our work suggests that deep learning is a practical and accurate alternative to the traditional way of generating approximate simulations of the structure formation of the Universe.

Materials and Methods

Dataset.

The full simulation data consists of 10,000 simulations of boxes with ZA and FastPM as input–output pairs, with an effective volume of 20 (Gpc/h)3 (), comparable to the volume of a large spectroscopic sky survey like Dark Energy Spectroscopic Instrument or EUCLID. We split the full simulation dataset into 80%, 10%, and 10% for training, validation, and test, respectively. We also generated 1,000 simulations for 2LPT for each set of tested cosmological parameters.

Model and Training.

The adopts the U-Net architecture (54) with 15 convolution or deconvolution layers and ∼ trainable parameters. Our generalizes the standard U-Net architecture to work with 3D data (69–71). The details of the architecture are described in the following sections, and a schematic figure of the architecture is shown in . In the training phase, we use the Adam Optimizer (72) with a learning rate of 0.0001, and first- and second-moment exponential decay rates equal to 0.9 and 0.999, respectively. We use the mean-squared error as the loss function (Loss Function) and regularization with regularization coefficient 0.0001.

Details of the Architecture.

The contracting path follows the typical architecture of a convolution network. It consists of two blocks, each of which consists of two successive convolutions of stride 1 and a down-sampling convolution with stride 2. The convolution layers use 33 filters with a periodic padding of size 1 (Padding and Periodic Boundary) on both sides of each dimension. Notice that at each of the two down-sampling steps, we double the number of feature channels. At the bottom of the , another two successive convolutions with stride 1 and the same periodic padding as above are applied. The expansive path of our is an inverted version of the contracting path of the network. (It includes two repeated applications of the expansion block, each of which consists of one up-sampling–transposed convolution with stride 1/2 and two successive convolutions of stride 1. The transposed convolution and the convolution are constructed with 33 filters.)

We take special care in the padding and cropping procedure to preserve the shifting and rotation symmetry in the up-sampling layer in expansive path. Before the transposed convolution, we apply a periodic padding of length 1 on the right, down, and back sides of the box [padding = (0,1,0,1,0,1) in pytorch], and after the transposed convolution, we discard one column on the left, up, and front sides of the box and two columns on the right, down, and back sides [crop = (1,2,1,2,1,2)].

A special feature of the is the concatenation procedure, where the up-sampling layer halves the feature channels and then concatenates them with the corresponding feature channels on the contracting path, doubling the number of feature channels.

The expansive building block then follows a 11 convolution without padding, which converts the 64 features to the final 3D displacement field. All convolutions in the network except the last one are followed by a rectified linear unit activation and batch normalization.

Padding and Periodic Boundary.

It is common to use constant or reflective padding in deep models for image processing. However, these approaches are not suitable for our setting. The physical model we are learning is constructed on a spatial volume with a periodic boundary condition. This is sometimes also referred to as a torus geometry, where the boundaries of the simulation box are topologically connected—that is, , where is the index of the spatial location, and is the periodicity (size of box). Constant or reflective padding strategies break the connection between the physically nearby points separated across the box, which not only loses information but also introduces noise during the convolution, further aggravated with an increased number of layers.

We find that the periodic padding strategy significantly improves the performance and expedites the convergence of our model, comparing to the same network using a constant padding strategy. This is not surprising, as one expects that it is easier to train a model that can explain the data than to train a model that does not.

Loss Function.

We train the to minimize the mean square error on particle displacementswhere labels the particles and N is the total number of particles. This loss function is proportional to the integrated squared error, and by using a Fourier transform and Parseval’s theorem, it can be rewritten aswhere is the Lagrangian space position, and its corresponding wavevector. is the transfer function defined in Eq. , and is the correlation coefficient defined in Eq. , which characterize the similarity between the predicted and true fields, in amplitude and phase, respectively. Eq. shows that our simple loss function jointly captures both of these measures: As and approach 1, the loss function approaches 0.

Data Availability.

The source code of our implementation is available at https://github.com/siyucosmo/ML-Recon. The code to generate the training data is also available at https://github.com/rainwoodman/fastpm.