Fast cosmic web simulations with generative adversarial networks
Andres C. Rodriguez Tomasz Kacprzak Aurelien Lucchi Adam Amara Raphael Sgier Janis Fluri Thomas Hofmann Alexandre Réfrégier
astro-ph.CO stat.ML

Gamow-Teller transitions from high-spin isomers in N=Z nuclei
    M. Honma
    December 30, 2023
=============================================================







§ INTRODUCTION


The large scale distribution of matter in the universe takes the form of a complicated network called the cosmic web <cit.>.

The properties of this distribution contain important cosmological information used to study the nature of dark matter, dark energy, and the laws of gravity <cit.>, as different cosmological models give rise to dark matter distributions with different properties.
Simulations of these cosmic structures <cit.> play a fundamental role in understanding cosmological measurements <cit.>.
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 large computational resources for a long period of time such as days or weeks <cit.>.
Furthermore, reliable measurements of cosmological parameters typically require a large number of simulations of various cosmological models <cit.>.
This creates a strong need for fast, approximate approaches for generating simulations of cosmic web <cit.>.


Here we demonstrate the possibility of using deep generative models to synthesize samples of the cosmic web.
Deep generative models <cit.> 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 uncorrelated to the training examples.
A 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 can potentially enable 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 <cit.>.
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 <cit.> or to recover certain features out of noisy astrophysical images <cit.>.
GANs were recently applied to generating samples of projected 2D mass distribution, called convergence <cit.>.
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.
For cosmological applications, it may be more efficient to work with the full 3D matter distributions generated by a GAN, rather then 2D convergence maps.
For gravitational lensing, the convergence map depends on the input distribution of background galaxies <cit.>; the 3D matter distribution is projected onto the sky plane by integrating the mass in radial direction against a lensing kernel, which depends on distribution n(z) of redshifts z of background galaxies.
For most lensing studies, the uncertainty on n(z) is large and is effectively marginalised over.
If the 3D distributions are simulated, then the projection can be done analytically, for a given n(z) <cit.>.
For a 2D generative model, a separate GAN would have to be trained for each n(z) distribution.
This may be particularly important for analyses beyond the power spectrum, such as peak statistics <cit.> or deep learning <cit.>, which use simulations to predict both the signal and its uncertainty.
In this paper we demonstrate the feasibility of GAN-based methods for capturing the type of matter distributions characteristic for in N-body simulations.
As the development of 3D generative methods for N-body data is likely to be a very challenging due to scalability issues and memory requirements, we consider this to be an important step in asserting that this approach is worth pursuing further.


In learning the cosmic web structures, which are more feature-rich than projected convergence maps, we encountered and addressed several important challenges.
The first was to handle data with very large dynamic range of the data; the density in the images created from slices of N-body simulations span several orders of magnitude.
Secondly, we explored how mode collapse, a feature of GANs causing the model to focuses on a single local minimum, affects the quality of results <cit.>.
As mode collapse is expected to depend on the degree of homogeneity between samples, we tested the performance of GANs for both large and small cosmological volumes, of size 500 and 100 Mpc;
the matter density distributions in large boxes are considerably more homogeneous than in small boxes.

Finally, expanding on the work of <cit.>, we additionally evaluate the cross-correlations of the GAN-generated data with itself and the training set.
A high cross-correlation would be an indication of lack of independence between the generated samples, a feature which we would judge to be undesirable in this task.



The paper is organised as follows.
In Section <ref> we describe the Generative Adversarial Networks.
Section <ref> contains the information on N-body simulations used.
Our implementation of the algorithm is described in Section <ref> and diagnostics used to evaluate its performance are detailed in Section <ref>.
We present the results in Section <ref> and conclude in Section <ref>.



§ 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: ↦ [0; 1] tries to probabilistically classify a sample  as being real or fake.
On the other hand, the generator G: ↦ 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  drawn from a given prior p_prior() which is typically a Gaussian distribution.
Formally, the two networks D and G play the following two-player minimax game:




    min_G max_D [ V(D, G) ]  
    
    V(D, G)     _x∼ p_data(x)[log D(x)]  + _z∼ p_prior(z)[log (1 - D(G(z)))],

where  is the expectation function.
The standard GAN approach <cit.> aims at finding a Nash Equilibrium of this objective by using gradient-based techniques in an alternating fashion, sometimes coupled with stabilization techniques <cit.>. As shown in <cit.>, for the Bayes-optimal discriminator D(x), the objective in Eq. <ref> reduces to the Jensen-Shannon divergence between p_data and the distribution induced by the generator. The work of <cit.> later generalized this to a more general class of f-divergences. An alternative formulation proposed in <cit.> 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 <cit.> to create 10 independent simulation boxes for both box sizes.
The cosmological model used was ΛCDM (Cold Dark Matter) 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:

    s(x) = 2x/(x+a)  - 1

where a 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 super-massive 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 a allows to control the median value of the images, and was fixed to a=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.

In this paper we used L-PICOLA: a faster, but approximate simulator.
For a real application of our method a more precise simulator would be used, such as GADGET-2 <cit.> or PkdGrav3 <cit.>.
Nevertheless, for the purpose of demonstration of performance of GANs, we consider L-PICOLA simulations to be sufficient.
We do not expect the results to differ much if GANs were trained on simulations generated using more precise codes.




§ IMPLEMENTATION AND TRAINING


We use a slightly modified version of the standard DCGAN architecture <cit.>, which was shown to achieve good results on natural images, including various datasets such as LSUN-Bedrooms (3 million indoor bedrooms images) <cit.> or the celebrity face dataset (CelebA, 200000 28x28 pixel celebrity faces) <cit.>.

Table <ref> 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 <cit.>, 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 generator might converge to them, leaving out parts of regions of the target distribution. Wasserstein-1 loss, has some empirical evidence to prevent mode collapse but still suffers from it.

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 <ref>.
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 <ref> 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.


Matter density distribution can be described as dimensionless over-density field in space δ( x) = ( ρ( x)-ρ̅) / ρ̅,
where ρ( x) is the matter density at position x and ρ̅ is the mean density in the universe.
The cross power spectrum P_δ_1δ_2 of the matter over-densities is calculated as follows

    ⟨δ̃_1(ℓ) δ̃_2^*(ℓ) ⟩ =
    (2π)^2 δ_D(ℓ-ℓ') P_×(ℓ).


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 <cit.>.
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.

We trained the model parameters using ADAM, a gradient based optimizer <cit.>, 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.
A quantitative assessment of the results was performed based on summary statistics commonly used in cosmology, described in Section <ref>.
The angular power spectrum is a standard measure used for describing the matter distribution <cit.>.
Another important statistic used for cosmological measurements is the distribution of maxima in the density distribution, often called “peak statistics” <cit.>.
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.
Another statistic we used was the histogram of pixel values of N-body and GAN-generated images.
Finally, we calculated the covariance between the power spectrum values at different k.





 §.§ Large images of size 500 Mpc


Figure <ref> presents the original images (top) and synthesized images (bottom), for the 500 Mpc simulations.
The plotted images were transformed using Equation <ref> to make it easier to assess the difference in the texture.
The cosmic web structure produced by the GAN model is visually very difficult to distinguish from the originals, even for human experts.
The GAN can capture the prominent features in this data, including halos and filaments.

Figure <ref> shows the summary statistics for the original (blue lines) and GAN-synthesized (red lines) samples for 500 Mpc images.
Mass density histograms, shown in the top left panel, agree well throughout most of the range, except for very large densities.
Peak statistics, shown in top right panel, also agree well, although slightly worse than the density histograms, especially for higher mass ranges, where the error can reach ∼10%.
The power spectrum is shown in the bottom left panel.
We focused on correlations at angular scales larger than a few Mpc,
as the current N-body simulations do not agree well in their predictions for smaller scales <cit.>.
We find that between k=0.06 and k=0.4 the agreement is 1-2%, while the rest of the range agrees within 15%, and for large k>0.8, the error starts to increase dramatically.
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 evident 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 uncorrelated realisations of cosmic web.
While the lack of cross correlation does not strictly imply independence, it assures that local structures are not consistently the same, and the data is not simply memorised and “pasted” in the same locations.
Finally, the correlation matrix of power spectra at different values of k is shown in the top panels of Figure <ref>.
	For the 500 Mpc images, the structure of the correlation matrix for GAN is similar to N-body: more correlation is observed at small scales. The numerical agreement, however, is good only for small scales, with ∼5% differences. For large scales, the errors reach 20%.












 §.§ Small images of size 100 Mpc


The example density distributions from 100 Mpc data is shown in Figure <ref>.
These images 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.
The agreement between the real and GAN-generated images is still good, although it is possible to distinguish them visually.
Notably, the filaments do not look as distinct as in the real images.
Even thought the images are not homogeneous, the network can still capture this variability: it does generate images full of structures, as well as rather empty ones.
However, the proportions of these types of distributions among the generated samples may differ between real and GAN data.
These differences will manifest themselves in the quantitative comparison.

Figure <ref> shows the summary statistics for 100 Mpc images.
The agreement between mass density histograms is good.
The difference in terms of peak statistics is on average small, although with deviations of ∼10 % in the middle of the mass range.
The error on the power spectrum is much larger: between k=0.13 and k=4 there is a 20% disagreement, and reaches 30% outside that range.
From from k>5 the error becomes large.
Similarly to 500 Mpc images, we do not observe large discrepancies in the cross power spectrum between pairs of GAN generated images, as well as between GAN-real pairs.
The agreement between the correlation matrix of the power spectra between N-body and GAN-generated is much better for 100 Mpc images.
	The differences are smaller than 5% for most of the correlation matrix, as shown in the bottom panels in Figure <ref>.








§ 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 produced from N-body simulations, and used it to generate a synthetic cosmic web.
Our GAN-generated images are visually very similar to the ones from N-body simulations: the generative model managed to capture the complex structures of halos, filaments and voids.
We compared the GAN-generated images to the N-body originals using several summary statistics and found a good agreement.
Most notably, for 500 Mpc, the power agreement on power spectrum was very good: between k=0.06 and k=0.4 the level of 1-2% is close to the requirements for precision cosmology <cit.>.
The correlation matrices of P(k) values had similar structures and agreed to around 5% at small scales, but the GANs did not reproduce the large scale correlations well, with ∼20% difference.
While more work would be needed to improve this agreement further, this result is  promising for using GANs as emulators of mass density distributions for practical applications.

For 100 Mpc images the error on the power spectrum was larger, reaching 20%.
We attribute this feature to the fact that images in the 100 Mpc sample are much more inhomogeneous than the 500 Mpc sample: some images contain dense regions with halos, and some relatively empty regions with few features.
We have seen empirically that this tends to induce a known phenomenon in GANs called mode collapse, where the training algorithm focuses on a subset of the target distribution .
This results in the model generating a few specific types of images, for example the ones with empty regions, more often than others.
We conclude that the application of GANs is suitable for large, homogeneous datasets.
For the type of inhomogeneous distributions appearing in the 100 Mpc sample, some techniques addressing mode-collapse <cit.> might be required if high-quality statistics are required.

An important advantage of the approach we presented here is that, once trained, it generates new samples in a fraction of a second on a modern Graphics Processing Unit (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 <cit.> or analyses using weak lensing shear peak statistics <cit.>.
Generative methods may become even more important in the future; the need for fast N-body simulations is anticipated to grow in the era of large cosmological datasets obtained by the Euclid[<www.euclid-ec.org>] and LSST[<www.lsst.org>] projects.
The need for fast simulations will be amplified further by the emergence of new analysis methods, which can be based on advanced statistics <cit.> or deep learning <cit.>.
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 <cit.> for estimating cosmological parameters from dark matter simulations.

Finally, it would be interesting to explore how many simulations are needed to train a GAN model for a given precision requirement.
Another future direction would be to further explore the agreement between the original and GAN-generated images in terms of advanced statistics, such as for example 3-pt functions or Minkowski functionals.
Going beyond the cross-correlations to further tests for independence of the GAN-generated samples could also be of interest.
We leave this exploration to future work.




This work was support in part by grant number 200021_169130 from the Swiss National Science Foundation.
We acknowledge the support of IT services of the Leonhard and Euler clusters at ETH Zurich.




unsrt





§ NETWORK ARCHITECTURE AND PARAMETERS


The architecture used in this work is shown in Table <ref>. See Section <ref> for more details.
Table <ref> contains the list of hyperparameters used.








