Fast Cosmic Web Simulations with Generative Adversarial Networks

. Dark matter in the universe evolves through gravity to form a complex network of halos, ﬁlaments, sheets and voids, that is known as the cosmic web. Computational models of the underlying physical processes, such as classical N-body simulations, are extremely resource intensive, as they track the action of gravity in an expanding universe using billions of particles as tracers of the cosmic matter distribution. Therefore, upcoming cosmology experiments will face a computational bottleneck that may limit the exploitation of their full scientiﬁc potential. To address this challenge, we demonstrate the application of a machine learning technique called Generative Adversarial Networks (GAN) to learn models that can eﬃciently generate new, physically realistic realizations of the cosmic web. Our training set is a small, representative sample of 2D image snapshots from N-body simulations of size 500 and 100 Mpc. We show that the GAN-produced results


Introduction
The large scale distribution of matter in the universe takes the form of a complicated network called the cosmic web [1][2][3][4].The properties of this distribution contain important cosmological information used to study the nature of dark matter, dark energy, and the laws of gravity [5][6][7], as different cosmological models give rise to dark matter distributions with different properties.Simulations of these cosmic structures [8,9] play a fundamental role in understanding cosmological measurements [10,11].These simulations are commonly computed using N-body techniques, which represent the matter distribution as a set of particles that evolve throughout cosmic time according to the underlying cosmological model and the laws of gravity.Creating a single N-body simulation requires the use of significant computational resources for a long period of time such as days or weeks [12,13].Furthermore, reliable measurements of cosmological parameters typically require a large number of simulations of various cosmological models [14,15].This creates a strong need for fast, approximate approaches for generating simulations of cosmic web [16][17][18].
Here we demonstrate the possibility of using deep generative models to synthesize samples of the cosmic web.Deep generative models [19,20] are able to learn complex distributions from a given set of data, and then generate new, statistically consistent data samples.Such a deep generative model can be trained on a set of N-body simulations.Once the training is complete, the generative model can create new, random dark matter distributions that are statistically independent from the training examples.A significant practical advantage of using a generative model is that the generation process is extremely fast, thus giving us the ability to generate a virtually unlimited number of samples of the cosmic web.Having access to such a large amount of simulations would yield to more reliable scientific studies and would therefore enhance our ability to understand the physics of the Universe.
In the last decade, Deep Learning approaches have achieved outstanding results in many fields, especially for computer vision tasks such as image segmentation or object detection [21].Deep convolutional neural networks (DCNN) have also recently been used as data generating mechanisms.Here a latent random vector, typically a high-dimensional Gaussian, is passed through a DCNN in order to output images.Generative Adversarial Networks (GAN) create such a model by adopting an adversarial game setting between two DCNN players, a generator and a discriminator.The goal of the generator is to produce samples resembling the originals while the discriminator aims at distinguishing the originals from the fake samples produced by the generator.The training process ends when a Nash equilibrium is reached, that is when no player can do better by unilaterally changing his strategy.
The rise of deep generative models has sparked a strong interest in the field of astronomy.Deep generative models have been used to generate astronomical images of galaxies [22][23][24] or to recover certain features out of noisy astrophysical images [24].GANs were recently applied to generating samples of projected 2D mass distribution, called convergence [25].This approach can generate random samples of convergence maps, which are consistent with the original simulated maps according to several summary statistics.The projection process, however, washes out the complex network structures present in the dark matter distribution.Here, we instead focus on generating the structure of the cosmic web without projection, therefore preserving the ability of the generative model to create halos, filaments, and sheets.We accomplish our goal by synthesizing thin slices of dark matter distribution which have been pixelised to create 2D images that serve as training data for a GAN model.A demonstration of this method on 2D slices presents a case for the development of deep learning methods able to generate full, 3D dark matter distributions.
The paper is organised as follows.In Section 2 we describe the Generative Adversarial Networks.Section 3 contains the information on N-body simulations used.Our implementation of the algorithm is described in Section 4 and diagnostics used to evaluate its performance are detailed in Section 5. We present the results in Section 6 and conclude in Section 7.

Generative Adversarial Networks
The basic idea behind GANs consists in pairing-up two neural networks: a generator network G and a discriminator network D. These networks are trained in an adversarial game setting.The discriminator D : x → [0; 1] tries to probabilistically classify a sample x as being real or fake.On the other hand, the generator G : z → x tries to generate samples that look like they were drawn from the true data distribution p data .This generator makes use of a random variable z drawn from a given prior p prior (z) which is typically a Gaussian distribution.Formally, the two networks D and G play the following two-player minimax game: where E is the expectation function.The standard GAN approach [20] aims at finding a Nash Equilibrium of this objective by using gradient-based techniques in an alternating fashion, sometimes coupled with stabilization techniques [26,27].As shown in [20], for the Bayes-optimal discriminator D(x), the objective in Eq. 2.1 reduces to the Jensen-Shannon divergence between p data and the distribution induced by the generator.The work of [28] later generalized this to a more general class of f-divergences.An alternative formulation proposed in [29] uses the Wasserstein-1 distance to measure how different the real and fake samples are.In this work we experimented with both the standard GAN approach as well as Wasserstein GAN.We found both approaches to produce similar results and here present the results for 500 Mpc using Wasserstein-1 distance and 100 Mpc using the standard GAN approach.

N-body simulations data
We created N-body simulations of cosmic structures in boxes of size 100 Mpc and 500 Mpc with 512 3 and 1,024 3 particles respectively.We used L-PICOLA [18] to create 10 independent simulation boxes for both box sizes.The cosmological model used was ΛCDM with Hubble constant H 0 = 100, h = 70 km s −1 Mpc −1 , dark energy density Ω Λ = 0.72 and matter density Ω m = 0.28.We used the particle distribution at redshift z = 0. We cut the boxes into thin slices to create grayscale, two-dimensional images of the cosmic web.This is accomplished by dividing the x-coordinates into uniform intervals to create 1,000 segments.We then selected 500 non-consecutive slices and repeated this process for the y and z axes, which gave us 1,500 samples from each of the 10 realizations, yielding a total of 15, 000 samples as our training dataset.We pixelised these slices into 256 × 256 pixel images.The value at each pixel corresponded to its particle count.After the pixelisation, the images are smoothed with a Gaussian kernel with standard deviation of one pixel.This step is done to decrease the particle shot noise.
Most existing GAN architectures are designed for natural images and therefore require an RGB representation with 3 channels and integer values between 0 and 255.We adapted the DCNN architecture to work on our grayscale, floating-point images.We scaled the image values to lie in the interval [−1, 1] as we empirically found this transformation to improve performance.Once we have trained our GAN model, newly generated samples are transformed back to the original range using an inverse transformation.The transformation between the original, smoothed image x and the scaled image s was chosen to be where k is a free parameter.This transformation is non-linear, and similar in nature to a logarithm function.This choice was motivated by the fact, that the cosmic web has a high dynamic range between empty regions of space (voids with no particles) and supermassive halos (with many, concentrated particles).This non-linear transformation enhances the contrast on features of interest, namely the network structure of filaments, sheets and halos.The parameter k allows to control the median value of the images, and was fixed to k = 4 throughout the experimental section.Immediately after the generation of a new, synthetic image, we apply the inverse function s −1 (x) to transform it to the original space.

Implementation and training
We use a slightly modified version of the standard DCGAN architecture [30], which was shown to achieve good results on natural images, including various datasets such as LSUN-Bedrooms (3 million indoor bedrooms images) [31] or the celebrity face dataset (CelebA, 200000 28x28 pixel celebrity faces) [32].Table 1 in Appendix A presents the details of the architecture used for our experiments.We used similar architectures for both the discriminator and the generator, consisting of five convolutional layers.The total number of trainable parameters in both networks is 3.2 • 10 7 .We trained the networks until we achieved convergence in terms of the discriminant score for the standard version and a stable distance between the generated and real images for Wasserstein-1.
A commonly faced problem when training GANs is a phenomenon called "mode collapse" [33][34][35], where the network focuses on a subset of the modes of the underlying data distribution.In these regions where the generator is fooling the discriminator well, the gradient signal becomes weak and the discriminator might be unable to properly lead the generator to the right target distribution.We addressed this problem by doing early stopping, effectively selecting the network parameters during the training process by choosing the network that displayed the best agreement in terms of the power spectrum statistics described in Section 5.This happened after 17 and 21 epochs (one epoch consists of one full training cycle over the training set) for the 500 and 100 Mpc images respectively, which took 16.1 and 7 hours on a single GPU Nvidia GTX 1080 with 8GB.Table 2 in Appendix A presents the set of hyperparameters used in our results.

Diagnostics
The diagnostic measures used in this work are: average histogram of pixel values in the images, average histogram of values of maxima ("peaks"), average auto power spectrum and the average cross-power spectrum of pairs of images within the sample.
(5.1) where δ1 ( ) and δ2 ( ) are the Fourier transforms of two over-density maps at each logarithmically spaced Fourier bin , and δ D is the Dirac delta function.To compute the auto power spectrum, we set δ1 ( ) = δ2 ( ).We compute both auto and cross power spectrum from 2D images using a discrete Fourier transform, followed by averaging over angles.One of the popular alternatives to power spectrum for analysing matter density distribution is the peak statistics.These statistics capture non-Gaussian features present in the cosmic web and are commonly used on weak lensing data [15,36].A "peak" is a pixel in the density map that is higher than all its immediate 24 neighbours.The peaks are then counted as a function of their height.

Results
We focused our study on two simulation regimes: large-scale distribution, simulated in boxes of size 500 Mpc, and small-scale distribution, with boxes of size 100 Mpc.For both configurations we ran 10 independent simulations.From these boxes, we cut out a total of 15,000 thin, 2D slices for each box size.We design a GAN model where both the discriminator and generator are deep convolutional neural networks.These networks consists of 5 layers, with 4 convolutional layers using filter sizes of 5 × 5 pixels.The total number of trainable parameters in both networks is 3.2 • 10 7 .We trained the model parameters using stochastic gradient descent, which yields a model that can generate new, random cosmic web images.
We assessed the performance of the generative model in several ways.First, we performed a visual comparison of the original and synthetic images.Figure 1 presents the original images (top) and synthesized images (bottom), for the 500 Mpc simulations.The cosmic web structure produced by the GAN model is visually very difficult to distinguish from the originals, even for human experts.Similar results were observed for the 100 Mpc simulations shown in Figure 3.
A quantitative assessment of the results was performed based on summary statistics commonly used in cosmology.The angular power spectrum is a standard measure used for describing the matter distribution [37].We focused on correlations at angular scales larger than a few Mpc, as the current of N-body simulations do not agree well in their predictions for smaller scales [38].Another important statistic used for cosmological measurements is the distribution of maxima in the density distribution, often called "peak statistics" [15,39].This statistic compares the number of maxima in the maps as a function of their values.We also assessed the statistical independence of samples generated by GANs, as real cosmic structures are expected to be independent due to isotropy and homogeneity of the universe, unless they are physically close to each other.To assess the independence of generated cosmic web distributions, we compare the cross-correlations between pairs of images.Finally, we compared the histograms of pixel values of N-body and GAN-generated images.
Figures 2 and 4 show the summary statistics for the original (blue lines) and GANsynthesized (red lines) samples, for 500 Mpc and 100 Mpc images, respectively.The top left panels show the dark matter density histograms, as calculated from the pixelised images.The top right panels show the density histograms of peaks in the matter distribution.Average power spectra are shown in the bottom right panels.These results demonstrate that the agreement between real and generated samples is generally very good across all summary N-body simulation samples Generative Adversarial Network (GAN) samples statistics.The exception is a small difference in the power spectrum for 100 Mpc images, on which we comment below.Finally, the bottom right panels show the average cross power spectra, with the coloured bands corresponding to the standard deviation calculated using all available image pairs.As expected, the cross power spectrum of the original images is close to zero.We do not find significant discrepancies in the cross power spectrum between pairs consisting of N-body-and GAN-generated image, as well as between pairs of GAN-generated images.This indicates that the generated images can be also considered as independent realisations of cosmic web.The cross power spectrum is calculated between pairs N-body images (blue points), between pairs of GAN images (red points), and between pairs consisting of one GAN and one N-body image (cyan points).The power spectra are shown in units of h −1 Mpc, where h = H 0 /100 corresponds to the Hubble parameter.
The small discrepancy between GAN-and N-body-generated images for 100 Mpc can be attributed to mode collapse (introduced in Section 4).Images of size 100 Mpc are less homogeneous than the ones of size 500 Mpc.The structures present in smaller images can vary from image to image: some may contain only empty space while some might be large structures.This may cause the GAN algorithm to focus preferentially on one type of images (for example, the ones with fewer large structures), therefore resulting in a slight mismatch in the power spectrum.Images of size 500 Mpc are much more homogeneous within the sample, which alleviates the mode collapse problem.

Conclusion
We demonstrated the ability of Generative Adversarial Models to learn the distribution describing the complex structures of the cosmic web.We implemented a generative model based on deep convolutional neural networks, trained it on 2D images of cosmic web from N-body simulations, and used it to generate synthetic cosmic web.Our GAN-generated images are visually very similar to the ones from N-body simulation: the generative model managed to capture the complex structures of halos, filaments and voids.We compared the GANgenerated images to the N-body originals using several summary statistics and found a good agreement.As the independence of the realisations is preserved, the GAN model can be used  to very efficiently sample from that distribution.We found that the GAN method may be susceptible to mode collapse when the cosmic web images are not homogeneous within the sample, as it was the case for 100 Mpc images.This results in a slight mismatch in average statistics between N-body-and GAN-generated images.
A significant advantage of the approach we presented here is that, once trained, it generates new samples in a fraction of a second on a modern GPU.Compared to a classical N-body technique, this constitutes a gain of several orders of magnitude in terms of simulation time.The availability of this approach has the potential to dramatically reduce the computational burden required to acquire the data needed for most cosmological analyses.Examples of such analyses include the computation of covariance matrices for cosmology with large scale structure [14] or analyses using weak lensing shear peak statistics [39].Generative methods may become even more important in the future; the need for fast N-body simulations is anticipated to grow significantly in the era of large cosmological datasets obtained by the Euclid1 and LSST2 projects.The need for fast simulations will be amplified further by the emergence of new analysis methods, which can be based on advanced statistics [40] or deep learning [41].These methods aim to extract more information from cosmological data and often use large simulation datasets.While we demonstrated the performance of GANs for 2D images using training on a single GPU, this approach can naturally be extended to generate 3D mass distributions [42]  Step size for the Adam optimizer Gradient penalty -1,000 Gradient penalty applied for Wasserstein-1 Parameter in s(x) to obtain the scaled images Table 2: Hyper-parameters used in our GAN implementations.Adam [43] is the algorithm used to estimate the gradient in our models.

Figure 1 :
Figure 1: Samples from N-body simulations (top two rows) and from our GAN model (bottom two rows) for a box size of 500 Mpc.Note that the transformation in Equation 3.1 with k = 20 was applied to the images shown above for better readability.
x N-body GAN x GAN N-body x GAN

Figure 2 :
Figure 2: Comparison of summary between N-body and GAN simulations, for box size of 500 Mpc.The statistics are: mass density histogram (upper left), peak count (upper right), power spectrum of 2D images (lower left) and cross power spectrum (lower right).The cross power spectrum is calculated between pairs N-body images (blue points), between pairs of GAN images (red points), and between pairs consisting of one GAN and one N-body image (cyan points).The power spectra are shown in units of h −1 Mpc, where h = H 0 /100 corresponds to the Hubble parameter.

Figure 3 :
Figure 3: Samples from N-body simulation (top two rows) and from GAN (bottom two rows) for the box size of 100 Mpc.In this figure, transformation S1 with k = 7 was applied.
x N-body GAN x GAN N-body x GAN

Figure 4 :
Figure 4: Comparison of summary statistics between N-body and GAN simulations, for a box size of 100 Mpc.The statistics are the same as in Figure 4.

Table 1 :
for estimating cosmological parameters from dark matter simulations.Architecture used in the discriminator and generator networks.We used a batch size of m = 16 samples.The neural network has ∼32 million trainable parameters.Parameters for our Wasserstein-1 distance implementation are shown in brackets.