# Cosmic Voids in Generated N-body Simulations

Olivia Curtis,1 Tereasa G. Brainerd,1 Anthony Hernandez2
1Department of Astronomy & Institute for Astrophysical Research, Boston
University, Boston, MA 02215, USA
2Department of Computer Science and Engineering, University of South Florida,
Tampa, FL 33620, USA
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

A Generative Adversarial Network (GAN) was used to investigate the statistics
and properties of voids in $\Lambda$CDM universes. A total of 15,000 2D
images, extracted from N-body simulations, were used to train the GAN, which
was then used to generate 15,000 novel 2D images. The total number of voids
and the distribution of void sizes is similar in both sets of images. However,
the generated images yield somewhat fewer small voids than do the simulated
images. In addition, the generated images yield far fewer voids with central
density contrast $\sim-1$. Because the generated images yield fewer of the
emptiest voids, the distribution of mean interior density contrast is
systematically higher for the generated voids than it is for the simulated
voids. The mean radial underdensity profiles of the largest voids are similar
in both sets of images, but systematic differences are apparent. On small
scales ($r<0.5r_{v}$), the underdensity profiles of the voids in the generated
images exceed those of the voids in the simulated images. On large scales
($r>0.5r_{v}$), the underdensity profiles of the voids in the generated images
exceed those of the voids in the simulated images. The discrepancies between
the void properties in the two sets of images are likely attributable to
systematic difficulties faced by deep learning techniques, such as neural
networks struggling to capture absolute patterns in the data.

###### keywords:

methods: numerical – methods: statistical – dark matter – large-scale
structure of Universe

††pubyear: 2021††pagerange: Cosmic Voids in Generated N-body Simulations–A

##  1 Introduction

Galaxies in the universe are distributed within an interconnected, large-scale
network of walls and filaments that stretch for hundreds of megaparsecs.
Between these structures lie vast, underdense regions of space known as voids.
Voids can reach up to $100h^{-1}$Mpc in diameter (see, e.g., Giovanelli &
Haynes 1991 and references therein), and they have the potential to serve as
excellent laboratories for testing the popular $\Lambda$ Cold Dark Matter
($\Lambda$CDM) model of structure formation. Due to their underdense nature,
voids are dominated by dark energy and they are only weakly influenced by the
non-linear effects of gravity. Because of this, the shapes and distributions
of voids provide tests of modified theories of gravity, as well as constraints
on the dark energy equation of state, inflationary models, and the sum of the
neutrino masses (see, e.g., Li et al. 2012, Clampitt et al. 2013, Sahlén 2019,
Cai et al. 2015, Falck et al. 2018).

Existing void catalogs (e.g., Hoyle & Vogeley 2004, Sánchez et al. 2016, Mao
et al. 2017) suffer from small number statistics since they are insufficiently
deep to fully sample cosmic structure on the largest scales. For example,
Sánchez et al. (2016) identified only 87 voids in the first 139 sq. degrees of
the Dark Energy Survey (e.g., Dark Energy Survey Collaboration, & Flaugher
2005). The relatively small number of voids in current observational catalogs
limits our ability to place strong observational constraints on fundamental
statistics such as void frequency and mean void radius. However, upcoming
surveys from the Vera C. Rubin Observatory, as well as the Nancy Grace Roman
and Euclid satellites, are expected yield exceptionally deep maps that will
sample structure in the universe on gigaparsec scales. This should yield
significant improvements in our understanding of the properties of voids in
the observed universe.

From the standpoint of large-scale structure theory, gigaparsec-scale
simulations are required in order to determine whether $\Lambda$CDM is able to
successfully reproduce the largest structures in the universe. While Gpc-scale
N-body mock catalogs do exist (e.g., Kim et al. 2011, Falck et al. 2021), the
largest voids are so rare that only a few dozen of these objects are found
within a single Gpc-scale simulation. N-body simulations of this size are
computationally expensive to run, making it challenging to use these types of
mock catalogs to fully investigate the statistics of the largest voids in
$\Lambda$CDM universes. In order to make progress on the theoretical
properties of the largest voids, it would be helpful to be able to use a
modeling technique that is capable of quickly producing independent, novel
catalogs of large-scale structure. One technique that may prove advantageous
for such applications is the use of Generative Adversarial Networks (GANs).

GANs originate within the deep learning community and they have already begun
to demonstrate their usefulness for astrophysical applications. For example,
Rodríguez et al. (2018) used GANs to generate novel 2D images of large-scale
structure up to scales of $500h^{-1}$Mpc. Similarly, Feder et al. (2020) and
Curtis & Brainerd (2020) showed that GANs are capable of generating novel 3D
density maps of the large-scale structure on megaparsec scales. Kodi Ramanah
et al. (2020) experimented with the optimization of a network architecture
that could generate 3D density maps in boxes with sidelength $500h^{-1}$Mpc.
In addition, Mustafa et al. (2019) used GANs to produce novel weak lensing
convergence maps, Schawinski et al. (2017) used GANs to recover features in
Hubble Space Telescope images beyond the deconvolution limit, Gan et al.
(2021) demonstrated the ability to eliminate the effects of seeing on images
obtained with the Subaru telescope, and Zingales & Waldmann (2018) showed that
GANs can be used to retrieve the spectra of exoplanet atmospheres (allowing
computationally expensive Bayesian analyses to be bypassed).

In order for GANs to develop into computationally inexpensive alternatives to
enormous, Gpc-scale N-body simulations, they must first demonstrate the
ability to reproduce the large-scale structure that is consistently produced
by smaller, more typical N-body simulations. In this paper we focus on the
statistics and properties of voids in $\Lambda$CDM universes, and we compare
results obtained from N-body simulations to results obtained from a suitably
trained GAN. The paper is organized as follows. In §2, we discuss our
computational methods. These include the N-body simulations that were used to
train the GAN and the details of the particular GAN that we adopted. We also
discuss the network architecture and optimisation of the network weights. In
§3 we discuss the underdensity probability function, the voidfinder algorithm
and various void properties. A discussion of our results is presented in §4.
Throughout, we adopt the following cosmological parameters: $h=0.7$,
$\Omega_{\Lambda}=0.7$, $\Omega_{m0}=0.3$, and $\sigma_{8}=0.9$.

##  2 Methods

###  2.1 N-body Simulations

Before a GAN can be used to generate novel density maps of $\Lambda$CDM
universes, it first must be taught the properties that are expected for the
large-scale structure. Teaching the GAN these expected properties is a process
known as “training”. Training the GAN requires a set of independent N-body
simulations, the results of which allow the GAN to learn the expected
properties of the large-scale structure and then to extrapolate from its
training in order to produce new images. Here, the training set consisted of
ten $\Lambda$CDM simulations. Each simulation adopted a cubical box with
periodic boundary conditions, a box length of $L=512h^{-1}$Mpc, and a total of
$512^{3}$ dark matter particles of mass $m_{p}=8.3\times
10^{10}h^{-1}M_{\odot}$. A unique set of initial conditions for each
simulation was generated using
N-GenIC111https://www.h-its.org/2014/11/05/ngenic-code/ (Springel et al.,
2005) and the well-tested Gadget-2222https://wwwmpa.mpa-
garching.mpg.de/gadget/ code (Springel, 2005) was used to evolve the
simulations from redshift $z=50$ to the present epoch. All simulations were
run on the Shared Computing Cluster at the Massachusetts Green High-
Performance Computing Center.

To create a training set for the GAN, 15,000 2D images were extracted from the
N-body simulations using the method adopted by Rodríguez et al. (2018). First,
each simulation was divided each into 1,000 slices along each of the three
primary axes. Each slice was then pixellated using a nearest neighbor mass
assignment scheme and 500 non-consecutive slices along each axis were
selected. Lastly, each image was smoothed with a Gaussian filter with a width
of 1 pixel.

Before being passed to the GAN for training, the density maps from the N-body
simulations must be normalized. This is necessary because normalizing reduces
the cost of transforming the latent space distribution into the space expected
of the true data distribution. Following Rodríguez et al. (2018), the density
maps from our N-body simulations were normalized into the range [-1,1) using

$s(x)=\frac{2x}{x+a}~{}~{}.$ (1)

Here, $a$ is an integer that controls the median pixel value of the normalized
map and the properties of $s(x)$ resemble those of a logarithm function. For
consistency with Rodríguez et al. (2018) we adopt $a=4$ throughout. Examples
of normalized 2D mass density images obtained from the N-body simulations are
shown in the top panels of Figure 1.

Figure 1: Examples of normalized matter density maps at redshift $z=0$. Each
image has a sidelength of $512h^{-1}$ Mpc and a thickness of $0.512h^{-1}$Mpc.
Top: Images obtained from N-body simulations. Bottom: Images obtained from a
trained Wasserstein GAN (WGAN). Here the local matter density has been
normalized to the range [-1,1) using Equation 1 (see text). A full colour
version of this figure is available in the online journal.

###  2.2 Generative Adversarial Networks

A GAN is a game between two deep convolutional neural networks: a
discriminator network ($D$) and a generator network ($G$). The discriminator,
$D:(\vv{\bm{x}},\theta_{D})\rightarrow[0;1]$, attempts to label a sample,
$\vv{\bm{x}}$, as being either ‘real’ or ‘fake.’ Here, ‘real’ means
$\vv{\bm{x}}$ is drawn from the true data distribution ($p_{data}$) and ‘fake’
means $\vv{\bm{x}}$ is a generated image drawn from the set of generated data
($p_{gen}$). The variable $\theta_{D}$ represents the hyper-parameters that
characterize $D$. The generator network,
$G:(\vv{\bm{z}},\theta_{G})\rightarrow\vv{\bm{x}}$, attempts to map a random
variable, $\vv{\bm{z}}$, to a sample $\vv{\bm{x}}$ that appears to be drawn
from $p_{data}$. The random variable is drawn from a latent prior space,
$p_{prior}(\vv{\bm{z}})$, which is generally taken to be a Gaussian.

The discriminator network is trained to maximize the probability of correctly
labelling samples that are drawn from the training set and samples that are
drawn from the generator network. Meanwhile, the generator network is trained
to minimize $\log_{10}(1-D(G(\vv{\bm{z}})))$; i.e., the probability of the
generator network producing a sample that the discriminator network mislabels
as ‘real’. Formally, $D$ and $G$ play what is known as the two-player minimax
game:

$\begin{array}[]{c c}min&max\\\
\textsuperscript{G}&\textsuperscript{D}\end{array}[L(G_{\theta_{G}},D_{\theta_{D}})]$
(2)

with loss function

$L(G_{\theta_{G}},D_{\theta_{D}}):=\mathbb{E}_{\vv{\bm{x}}\sim
p_{data}(\vv{\bm{x}})}[\log_{10}D_{\theta_{D}}(\vv{\bm{x}})]+\\\
\mathbb{E}_{\vv{\bm{z}}\sim
p_{prior}(\vv{\bm{z}})}[\log_{10}(1-D_{\theta_{D}}(G_{\theta_{G}}(\vv{\bm{z}})))]~{},$
(3)

where $\mathbb{E}$ is the expectation function. (Note: Equation 3 reduces to
the Jensen-Shannon divergence between $p_{data}$ and $p_{gen}$; see Goodfellow
et al. 2020.)

####  2.2.1 Wasserstein GANs

For our work, we adopted a particular type of GAN known as a Wasserstein GAN.
The motivation for this choice stems from the phenomenon of mode collapse.
Mode collapse is said to have occurred in a generative network when the
generator produces only one particular output. This happens when the generator
discovers one output that consistently deceives the discriminator, at which
point the generator stops learning. For example, suppose a GAN is being
trained to generate pictures of rooms that would typically be found in a house
(see, e.g., Gulrajani et al. 2017). If the generator stops learning after it
discovers that one particular generated image of a bathroom consistently
deceives the discriminator, mode collapse occurs and the GAN will simply
continue to generate that same image of the bathroom as its sole output.

Wasserstein GANs mitigate the problem of mode collapse by adopting what is
known as an ‘Earth mover’s distance’ loss function (see, e.g., Arjovsky et al.
2017). In general, the loss function is a function that accumulates error, and
the network attempts to optimize the loss function during the training
process. The error is distributed to the network’s weights, allowing training
to occur. An intuitive explanation of the Earth mover’s distance is found in
its name. Imagine trying to relocate a pile of dirt from one location to
another. The Earth mover’s distance is then the metric that provides the most
cost effective way to transport the dirt.

When applied to GANs, the Earth mover’s distance allows the generator to find
the most efficient way of transforming $p_{gen}$ into $p_{data}$. In return,
the discriminator network now scores the ‘realness’ and ‘fakeness’ of an image
(i.e., it scores how well $p_{gen}$ resembles $p_{data}$). For this reason,
Arjovsky et al. (2017) refer to the discriminator network as a ‘critic
network’, a terminology that we also adopt in this paper. WGANs address the
problem of mode collapse by not rewarding the generator network for its
ability to generate a single result that consistently deceives the
discriminator network. Instead, the generator network attempts to optimize its
ability to produce results that receive a high score of ‘realness’ from the
critic network.

Additional discussion of Wasserstein GANs and the Earth Mover’s distance is
given in Appendix A.

####  2.2.2 Adam Optimizer

In order to carry out back propagation, we use a standard Adaptive Moment
Optimizer known as Adam (see Kingma & Ba 2015). Adam is a stochastic gradient
descent algorithm that updates the network weights using running averages of
the first and second moments of the gradients of the loss function. The Adam
optimizer is a popular choice for performing stochastic gradient descent
because it is computationally efficient and has low memory requirements.

For given a set of parameters, $w^{(t)}$, and a loss function, $L^{(t)}$, at
training iteration, $t$, the weights at the $t+1$ training iteration are found
by first calculating the biased first and second moments using Equations 4 and
5, respectively:

$m_{w}^{t+1}\leftarrow\beta_{1}m_{w}^{t}+(1-\beta_{1})\nabla_{w}L^{t}$ (4)
$v_{w}^{t+1}\leftarrow\beta_{2}v_{w}^{t}+(1-\beta_{2})(\nabla_{w}L^{t})^{2}.$
(5)

Here, $\beta_{1}$, $\beta_{2}$ $\in$ [0,1) are the exponential decay rates for
the moment estimates, $m_{w}^{(t)}$ are the biased first moment estimates of
the weights at step t, $v_{w}^{(t)}$ are the biased second moment estimates,
and $\nabla_{w}$ is the gradient operator with respect to the network
parameters.

Typical values for $\beta_{1}$ and $\beta_{2}$ are 0.9 and 0.999, so if
$m^{0}_{w}$ and $v^{0}_{w}$ are initialized to zero, a significant bias
towards small values in future updates will occur. This bias can be eliminated
(making the moments less sensitive to $\beta_{1}$ and $\beta_{2}$ at early
timesteps) by normalizing the $\beta$’s such that
$\sum_{i=0}^{t}\beta^{i}=\frac{1-\beta^{t}}{1-\beta}$. This normalization can
then be removed by dividing the first and second moments by $1-\beta^{t+1}$ in
Equations 6 and 7 below:

$\hat{m}_{w}=\frac{m_{w}^{t+1}}{1-\beta_{1}^{t+1}}$ (6)
$\hat{v}_{w}=\frac{v_{w}^{t+1}}{1-\beta_{2}^{t+1}}~{}.$ (7)

Finally, the weights at step $t+1$ are calculated from the unbiased moments
using

$w^{t+1}\leftarrow
w^{t}-\alpha\frac{\hat{m}_{w}}{\sqrt{\hat{v}_{w}}+\epsilon},$ (8)

where $\alpha$ is the algorithm’s learning rate and $\epsilon$ is a small
number that is used to prevent division by zero. Table 1 lists the values of
the hyperparameters that we adopted for our work.

Table 1: Hyperparmeters used in the Adam optimizer, as well as other important network parameters (see Table 2). Hyperparameter | Value | Description
---|---|---
Batch Size | 16 | Number of training samples fed into the network before back propagation
z dimension | 200 | Dimension of the random Gaussian vector
$\alpha_{D}$ | 1e-6 | Learning rate for the discriminator’s (D) Adam optimizer
$\beta_{1,D}$ | 0.5 | First exponential decay rate for D’s Adam optimizer
$\beta_{2,D}$ | 0.999 | Second exponential decay rate for D’s Adam optimizer
$\alpha_{G}$ | 5e-5 | Learning rate for the generator’s (G) Adam optimizer
$\beta_{1,G}$ | 0.5 | First exponential decay rate for G’s Adam optimizer
$\beta_{2,G}$ | 0.999 | Second exponential decay rate for G’s Adam optimizer
$\alpha$ | 0.8 | Momentum used in batch normalization
$n_{disc}$ | 5 | Number of times D is updated for every G update
c | 0.01 | Weight clipping parameter used to uphold Lipschitz continuity in WGAN architectures
$\alpha_{LeakyReLU}$ | 0.2 | Negative slope coefficient for LeakyReLU activation functions

####  2.2.3 Network Architecture

Throughout, we used a standard WGAN architecture with deep convolutional
layers. Details of the architecture are summarized in Table 2. The critic
takes a monochromatic image of size $256\times 256$ pixels as its input. Four
2D convolutional layers then down-sample the image. The output of the critic
is a score in the range $(-1,1)$, which describes how ‘real’ ($-1$) or ‘fake’
($+1$) the critic thinks the image is. For the convolutional layers, a
$4\times 4$ pixel kernel with a stride of size $(2,2)$ is used. In addition,
the inputs of the network are batch normalized before each mini-batch update.
The first four convolutional layers have a Leaky Rectified Linear Unit (Leaky-
ReLU) activation function, and the last layer has a linear activation
function. The first hidden layer down-samples the monochromatic image into a
tensor of shape $128\times 128\times 64$, where the last dimension is the
number of channels. Each consecutive hidden layer then doubles the number of
channels to 128, 256, and 512 respectively.

The generator takes a random 1D Gaussian vector of length 200 as its input.
The first hidden layer is a linear layer that reshapes the latent space vector
into a tensor of shape $16\times 16\times 512$ (i.e., 512 low resolution
images of size $16\times 16$ pixels). Four 2D convolutional transpose layers
up-scale the image while reducing the number of channels. The first four
layers use a ReLU activation function and have their inputs batch normalized
before each mini-batch update. The output of the generator is a monochromatic
image of size $256\times 256$ pixels where each pixel is cast into the range
$(-1,1)$ using a hyperbolic tangent function. To initialize the weights of the
critic and generator, a Glorot weight initialization (see, e.g., Glorot &
Bengio 2010) was used.

The networks were implemented using Python’s
TensorFlow333https://www.tensorflow.org/ package (e.g., Oliphant 2007; Abadi
et al. 2016) and were trained on two NVIDIA Tesla P100 graphics processing
units with compute capability 6.0. The critic and generator networks had 2.89
million and 16.00 million trainable parameters, respectively. Training was
carried out with Python version 3.7.7, TensorFlow version 1.15.0, CUDA version
10.0, and NVIDIA cuDNN version 7.6, and took approximately four hours for 20
training epochs. The implementation of the WGAN described in this paper is
available on the first author’s GitHub page.444https://github.com/o-curtis/

Table 2: WGAN network architecture used in this paper. Layer | Operation | Output | Dimension
---|---|---|---
Critic |  |  |
X |  |  | 256 x 256 x 1
$D_{0}$ | Conv2D | LeakyReLU-BatchNorm | 128 x 128 x 64
$D_{1}$ | Conv2D | LeakyReLU-BatchNorm | 64 x 64 x 128
$D_{2}$ | Conv2D | LeakyReLU-BatchNorm | 32 x 32 x 256
$D_{3}$ | Conv2D | LeakyReLU-BatchNorm | 16 x 16 x 512
$D_{4}$ | Linear | Linear | 1
Generator |  |  |
Z |  |  | 200 x 1
$G_{0}$ | Linear | ReLU-BatchNorm | 16 x 16 x 512
$G_{1}$ | Conv2DTranspose | ReLU-BatchNorm | 32 x 32 x 256
$G_{2}$ | Conv2DTranspose | ReLU-BatchNorm | 64 x 64 x 128
$G_{3}$ | Conv2DTranspose | ReLU-BatchNorm | 128 x 128 x 64
$G_{4}$ | Conv2DTranspose | Tanh | 256 x 256 x 1

##  3 Results

The WGAN was trained for 20 epochs and the network parameters were saved after
each training epoch. After 20 epochs, the outputs of each saved network were
analyzed to determine which training epoch produced matter power spectra that
were most similar to the training set. In the case of our WGAN, we found this
to occur after four training epochs. This network was then used to generate a
set of 15,000 images which were used for the analyses we present below. Two of
these generated images are shown in the bottom panels of Figure 1.

Figure 2 shows the results for the matter power of spectra of our simulated
and generated datasets. The power spectra were calculated in the standard way
(see, e.g., Davis et al. 1985) and they encompass distance scales that are
typical of void sizes (i.e., $10-100h^{-1}\rm{Mpc})$. Points in Figure 2 show
the mean power spectra, computed using all simulated and generated images, and
error bars were obtained from 10,000 bootstrap resamplings of the data. From
Figure 2, $P(k)$ from the generated images agrees with $P(k)$ from the
simulated images on scales $<20h^{-1}\rm{Mpc}$. On larger scales, however, the
generated images yield systematically lower values of $P(k)$ than do the
simulated images.

Figure 2: Mean matter power spectra, $P(k)$, computed from 15,000 simulated
images (circles) and 15,000 images generated with a WGAN (diamonds). Error
bars were computed using 10,000 bootstrap resamplings of the data.

###  3.1 Underdensity Probability Function

In the observed universe, the Void Probability Function (VPF) is commonly used
to determine whether a randomly selected region in an image is devoid of
galaxies (e.g., White 1979; Lachieze-Rey et al. 1992; Vogeley et al. 1994).
The VPF depends on the galaxy $N$-point correlation function via

$P_{0}(N,A)=\\\
exp\left[\sum_{N=1}^{\infty}\frac{(-n)^{N}}{N!}\int_{A}w_{N}(\vv{\bm{x}}_{1},...,\vv{\bm{x}}_{N})d^{2}x_{1}...d^{2}x_{N}\right]$
(9)

(e.g., White 1979), where $n$ is the average number density of galaxies in the
field, $w_{N}$ are the $N$-point correlation functions, and $\vv{\bf{x}}_{i}$
are the i-th galaxy positions in the area $A$.

Because our GAN outputs consist of images of a smoothed mass density
distribution (i.e., not the locations of individual galaxies or individual
particles) we do not compute the VPF. Instead, we compute a similar statistic
known as the Underdensity Probability Function (UPF). The UPF measures the
probability that, within an area of radius $r$, a randomly selected region of
space is less dense than some particular threshold density. To identify
underdense regions in our simulated and generated images, we define a
normalized density contrast $\Delta\equiv\frac{\delta\rho}{\rho}$ and we adopt
a threshold for the mass density contrast of $\Delta<-0.6$. Here $\rho$ is the
local matter density and $\delta\rho$ is the difference between the local
matter density and mean matter density of the image. That is, the UPF is
defined as $U(r)=P(\Delta<\Delta_{\rm crit};r)$ where $\Delta_{\rm
crit}=-0.6$.

The UPF for the simulated and generated images was determined by computing the
mean mass density contrast within randomly-placed circles of radius $r$. If
the average density contrast within a circle of radius $r$ is less than
$\Delta_{\rm crit}$, then that particular region qualifies as being
sufficiently underdense to be included in the calculation of the UPF. A total
of $N_{\rm test}=10^{4}$ random circles of radius $r$ were placed within each
image and the total number of underdense regions, $N_{0}$, was then computed.
The UPF was then defined as $U(r)=N_{(0;~{}\Delta<\Delta_{\rm crit})}/N_{\rm
test}$.

Figure 3 shows the mean UPF obtained from 15,000 simulated images and 15,000
generated images. The UPFs were calculated following the steps described above
and error bars were calculated by bootstrap resampling of both sets of
distributions. Note that, since smaller circles in an image are encompassed by
larger circles in the same image, the data points in Figure 3 are not
independent of one another. Hence, the point-to-point scatter is considerably
smaller than the formal, bootstrapped error bars. Similar to the power
spectra, there is good agreement between between the simulated and WGAN images
on small scales. On larger scales, there is a significantly lower chance of
finding underdense regions in the generated data than in the simulated data,
with the difference increasing for regions $\geq 60h^{-1}$ Mpc.

Figure 3: Mean underdensity probability function for images obtained from
N-body simulations (circles) and WGAN-generated images (diamonds). Error bars
were computed from 10,000 bootstrap resamplings of the data and are omitted
when they are comparable to or smaller than the data points.

###  3.2 Void Identification

To identify voids in both the simulated and generated density fields, we
adopted a voidfinder algorithm that is similar to the algorithm used by El-Ad
& Piran (1997) and Hoyle & Vogeley (2004) to identify voids in catalogs of
galaxies. When applied to an observational dataset, the voidfinder algorithm
identifies voids using the observed locations of galaxies. In our case there
are no galaxies in the images. Instead, the images consist of pixels with
known values of the mass density. Our version of the voidfinder algorithm,
therefore, scans images of the large-scale structure as defined by the
pixellated mass density. The algorithm yields a set of individual void centres
and a list of which pixels in a given image belong to which voids.

Our voidfinder algorithm can be summarized in four steps: [1] each pixel is
classified as being either a ‘wall pixel’ or a ‘void pixel’, [2] the distance
between every void pixel and the nearest wall pixel is calculated, [3] void
centres are defined to be the locations of the most isolated void pixels
(i.e., the void pixels that are farthest from a wall pixel), and [4] as part
of the process of identifying underdense regions that are interconnected, the
areas of the larger voids are ‘enhanced’ to include overlapping regions of low
mass density. The first step is performed by scanning over every pixel in the
image using a circular aperture of radius $d=3$ pixels that takes the periodic
boundary into account. If the pixel that the circular aperture is centred on
contains at least $N=9$ pixels with a normalised density contrast $>-0.6$, the
pixel is classified as a ‘wall pixel’. That is, if $\sim d^{2}$ of the pixels
in the neighbourhood around the central pixel are sufficiently dense compared
to the mean, the central pixel is classified as a wall pixel. Any pixel that
is not classified as a wall pixel is classified as a ’void pixel.’ The second
step is to iterate over each void pixel and determine its distance to the
nearest wall pixel. This is done by finding the circle of maximum radius that
can be centred around each void pixel such that the circle contains no wall
pixels. Following the terminology of Hoyle & Vogeley (2004), we refer to these
maximal circles as ‘holes’.

Once all of the pixels have been classified as being either a void pixel or a
wall pixel, the void centres are then identified. To do this, we begin by
sorting all holes in a given image from largest to smallest in terms of their
radii. We then define the centre of the largest hole in the image to be the
centre of the first void. Having found this first void centre, we then iterate
over all holes. If the next hole in the sorted list overlaps with any
previously-identified void by more than 10% of its own area, that hole is not
considered to be an individual void and, at this point, is temporarily
ignored. If the next hole in the sorted list does not overlap with any
previously-identified void (or if the overlap is $\leq 10$% of its own area),
the hole is classified as distinct void, with the centre of the hole being the
void centre. This process continues until all holes with radius $>3$ pixels
have been considered.

In order to account for the fact that underdense regions of space are
interconnected, the final step in the process of identifying the voids is to
consider the degree of overlap between holes and voids. The process is
intended to examine the underdense regions near the edges of the holes and
either assign those underdense regions to appropriate voids or reject them as
being part of a void. If a hole overlaps a single void by $>50$% of its own
area, then we consider the hole to be part of the void. If a hole overlaps
multiple voids by $>50$% of its own area, we consider it to be part of the
largest overlapping void. If a hole overlaps a void by $\leq 10$% of its own
area, then the hole is identified as a distinct void, but the area of overlap
is assigned to the larger void. If the overlap between a hole and a void is
$>10$% of the hole’s area but $<50$%, the pixels in the hole are not assigned
to any void. By this process, the total areas of the larger voids are enhanced
relative to the areas of the maximal circles that are centred on the void
centres.

Figure 4 shows the results of the voidfinder algorithm as applied to a single
image. The different panels of Figure 4 show: a) the normalised density
contrast, b) the locations of wall pixels (white) and void pixels (black), c)
the distance to the nearest wall pixel at each location (i.e., the maximal
radius for each pixel), d) the largest holes (i.e., the circles that define
the void centres) overlaid on the density contrast, e) the largest holes
overlaid on the wall/void pixel image, f) the locations of the largest holes
overlaid on the distance to the nearest wall pixel, g) the final void areas
overlaid on the wall/void pixel image, and h) same as g), but with the largest
holes indicated. The 15,000 images from the N-body simulations yielded a total
of 2.5 million voids and the 15,000 generated images yielded a total of 2.8
million voids.

###  3.3 Void Properties

Below we will adopt the following terminology when discussing the void
properties. The ‘radius’ of a void is simply the radius of the void’s maximal
circle (i.e., the largest circle, centred on the void centre, with radius
equal to the distance to the nearest wall pixel). The ‘effective radius’ of a
void corresponds to the radius of a circle that, after the area of the void
has been adjusted for any overlaps, is equal to the final area assigned to the
void. Since some of the smaller voids lose area to neighbouring large voids,
their effective radii will be less than the radii of their maximal circles.
Conversely, the larger voids with areas that have been enhanced by the
addition of pixels from neighbouring small voids (or significantly overlapping
holes) will have effective radii that are larger than the radius of their
maximal circle.

Figure 4: Illustration of the voidfinder algorithm. a) Normalised mass
density. b) Locations of wall pixels (white) and void pixels (black). c)
Maximal radius for each void pixel. d) – f) Largest hole contained within each
void, overlaid on the panels directly above (red circles). g) Final void areas
overlaid on the wall/void pixel distribution. h) Same as (g), but with the
largest hole contained within each void superposed on the image. A full colour
version of this figure is available in the on-line journal.

The top panels of Figure 5 show normalised probability distributions for: a)
the void radii and b) the effective void radii. Overall there is good
agreement between the results for the simulated and generated images. Compared
to the simulated images, however, the generated images produce fewer of the
smallest voids. The median void radius in the generated images is identical to
that of the simulated images ($10.0h^{-1}$ Mpc), but the median effective void
radius in the generated images ($10.0h^{-1}$ Mpc) is somewhat larger than that
in the simulated images ($9.6h^{-1}$ Mpc). Figure 5c) shows the normalised
probability distributions for the density contrast at the centres of the
voids. Here, again, the distributions are similar and the median value of the
central density contrast is nearly identical for the simulated and generated
images ($-0.83$ in the simulated images vs. $-0.81$ in the generated images).
However, compared to the simulated images, the generated images produce far
fewer voids for which the central density contrast is $\sim-1$ (i.e., the most
underdense void centres). Conversely, the generated images produce somewhat
more voids with central density contrasts between $\sim-0.9$ and $\sim-0.5$
than do the simulated images. Figure 5d) shows the normalised probability
distributions for the mean density contrast within the voids. While the median
values of the distributions in Figure 5d) are similar ($-0.79$ for the
simulated images vs. $-0.76$ for the generated images), there is a clear
offset between the two distributions such that the mean density contrast of
the generated voids is systematically higher than that of the simulated voids.

Figure 5: Normalized probability distributions of various void properties. a)
Void radii. b) Void effective radii. c) Value of the density contrast at the
centres of the voids. d) Mean interior density contrast, computed using the
final areas of the voids. Error bars are omitted because they are comparable
to or smaller than the sizes of the data points.

For the remainder of our analyses we will focus on a comparison of the best-
resolved voids. These are the largest voids in the images, and here we will
restrict the analysis to only those voids with radii $\geq 10$ pixels
(corresponding to a radius $\geq 20h^{-1}$ Mpc). Because large voids are rare,
these represent only a small fraction of the total number of voids in the
samples (15,063 in the simulated sample and 16,165 in the generated sample).
However, unlike the smallest voids which are poorly resolved, these largest
voids are sufficiently well resolved to allow a computation of their radial
density profiles.

For comparison with our complete samples of voids shown in Figure 5, Figure 6
shows the same normalised probability distributions, but using only the
largest voids. Compared to the results in Figure 5, there is a much greater
difference between the distributions of void radii and void effective radii in
Figure 6. This is due to a combination of two effects that go into assigning
each void its final area (from which the effective radius is determined).
First, the larger holes have a larger number of underdense pixels just outside
their circumferences than do the smaller holes. Therefore, a large hole that
has been identified as being a void has a higher probability of having its
final area enhanced relative to the area of the hole (i.e., by the addition of
pixels from smaller holes that overlap the large hole by a significant
amount). This necessarily increases the effective radius of a void with an
enhanced area. In addition, because of their relatively large size, the
probability of the largest voids overlapping by a small amount is higher than
the probability of the smallest voids overlapping by a similarly small amount.
In the case of these small overlaps ($\leq 10$% of the area of the smallest of
the two associated holes), the overlap area is assigned to the larger void,
resulting in the smaller of the two voids having a final area that is reduced
compared to the size of its hole (and hence, having an effective radius
smaller than the radius of its associated hole). From panels a) and b) of
Figure 6, the median values of the void radii ($20.0h^{-1}$ Mpc) and void
effective radii ($21.4h^{-1}$ Mpc) are identical in the simulated and
generated images.

Figure 6c) shows trends that are similar to those in Figure 5c). While the
median central densities of the largest voids are similar in both the
simulated and the generated images ($-0.82$ for the simulated images vs.
$-0.81$ for the generated images), the generated images produce far fewer of
the largest voids with central densities $\sim-1$. Like the complete void
sample, the generated images also produce somewhat more voids with central
densities larger than $\sim-1$, but the range over which this occurs
($\sim-0.9$ to $\sim-0.75$) is smaller than in the full sample. Figure 6d)
also shows trends that are similar to those in Figure 5d). The median values
of the distributions in Figure 6d) are similar ($-0.80$ for the simulated
images vs. $-0.76$ for the generated images), and there is a clear offset
between the two distributions such that the mean density contrast of the
largest generated voids is systematically higher than that of the largest
simulated voids.

Figure 6: Same as Figure 5, but here only voids with radius $\geq 20h^{-1}$
Mpc are used.

Lastly, Figure 7 shows the mean radial underdensity profiles for the largest
voids, plotted in terms of dimensionless distance relative to the void radius,
($r/r_{v}$). As expected, the central regions of the voids are extremely
underdense and the density contrast increases with radius from the void
centres. While the radial underdensity profiles are similar in both the
simulated and generated images, there are clear differences. For distances
$r\lesssim 0.6r_{v}$, the underdensities of the generated voids are somewhat
higher than the underdensities of the simulated voids. For distances $r\gtrsim
0.6r_{v}$, the sense of the of disagreement reverses, with the underdensities
of the generated voids being somewhat lower than the underdensities of the
simulated voids.

Figure 7: Mean radial underdensity profiles of largest voids, plotted as a
function of dimensionless distance relative to the void centres, ($r/r_{v}$).
Error bars are omitted because they are comparable to or smaller than the data
points.

##  4 Discussion

We have explored the use of a standard deep convolutional Wasserstein GAN
(WGAN) for the purposes of generating novel 2D images of the smoothed mass
density field of $\Lambda$CDM universes. The WGAN was trained using 15,000 2D
images of the large-scale structure that were obtained from N-body
simulations. The trained WGAN was used to generate 15,000 2D images and a
voidfinder algorithm, modified to identify regions of low mass density (as
opposed to regions of low galaxy density), was then used to create void
catalogs from the simulated and generated images. The simulated and generated
images yield a similar number of voids (2.5 million voids in the simulated
images vs. 2.8 million voids in the generated images) and the properties of
the voids in the simulated and generated images are also similar. Systematic
differences between the properties of the simulated and generated voids are,
however, apparent. Compared to the simulated images, the generated images
yield somewhat fewer of the smallest voids and significantly fewer voids with
central density contrast $\sim-1$. Overall, the generated images yield fewer
of the emptiest voids and, as a result, the distribution of mean interior
density contrast for the generated voids is offset systematically from the
distribution of mean interior density contrast for the simulated voids.

Neural networks function by learning the average trends in the target
distribution but they struggle to capture absolute patterns in the data (see,
e.g., Li et al. 2018 for a detailed discussion). That is, our deep
convolutional WGAN is clearly able to detect trends in the pixel-to-pixel
variations of the smoothed density contrast, but it is not reliably capturing
the absolute highs and lows. This effect can be seen in Figures 5c) and 6c),
where the WGAN images produce significantly fewer voids with central density
constrast $\sim-1$ than do the simulated images. The failure of the WGAN to
reliably capture the absolute highs and lows of the smoothed density contrast
also explains the systematic offsets between the distribution of mean interior
density contrast for the generated and simulated images in Figures 5d) and
6d).

GANS present an opportunity to investigate the distribution of large-scale
structure on Gpc-scales across multiple epochs of the universe without the
need for computationally-expensive Gpc-scale N-body simulations. In order for
GANs to become a truly viable alternative to simulations, they need to be able
to fully reproduce the structure that is seen in typical N-body simulations.
While our WGAN produced voids with properties that were in broad general
agreement with the properties of voids in N-body simulations, the systematic
differences between the populations are an indication that improvements in the
WGAN approach are still necessary.

For future work, the network architecture will need to be optimised in such a
way that it allows the trained WGAN to better capture absolute trends in the
target images. Such optimisation may be dependent on the size of the
convolutional kernels and/or the values of the hyperparameters used in the
generator and critic networks. The use of inception blocks (see Szegedy et al.
2015), for example, would help the network to pick up features of varying
sizes by scanning the input layer with kernels of varying sizes. This alone
might resolve the issue of the centres of the generated voids being
insuffciently underdense. Finally, it will be interesting to see whether more
sophisticated network architectures, such as the one used by Kodi Ramanah et
al. (2020) to generate large-scale structure maps, yield improved results for
the statistics of generated voids.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Abadi et al. (2016) Abadi M., et al., 2016, in Keeton K., Roscoe T., eds, Proc. of the 12th USENIX Symposium on Operating Systems Design and Implementation. USENIX Association, Berkeley, California, USA, p. 265
  * Arjovsky et al. (2017) Arjovsky M., Chintala S., Bottou L., 2017, in Precup D., Teh Y. W., eds, Proc. of Machine Learning Research Vol. 70, Proc. of the 34th International Conference on Machine Learning. Journal of Machine Learning Research, Inc., Brookline, Massachusetts, USA, p. 214
  * Cai et al. (2015) Cai Y.-C., Padilla N., Li B., 2015, MNRAS, 451, 1036
  * Clampitt et al. (2013) Clampitt J., Cai Y.-C., Li B., 2013, MNRAS, 431, 749
  * Curtis & Brainerd (2020) Curtis O., Brainerd T. G., 2020, RNAAS, 4, 90
  * Dark Energy Survey Collaboration, & Flaugher (2005) Dark Energy Survey Collaboration, Flaugher B., 2005, International Journal of Modern Physics A, 20, 3121
  * Davis et al. (1985) Davis M., Efstathiou G., Frenk C. S., White S. D., 1985, ApJ, 292, 371
  * El-Ad & Piran (1997) El-Ad H., Piran T., 1997, ApJ, 491, 421
  * Falck et al. (2018) Falck B., Koyama K., Zhao G.-B., Cautun M., 2018, MNRAS, 475, 3262
  * Falck et al. (2021) Falck B., Wang J., Jenkins A., Lemson G., Medvedev D., Neyrinck M. C., Szalay A. S., 2021, preprint (arXiv:2101.03631)
  * Feder et al. (2020) Feder R. M., Berger P., Stein G., 2020, Phys. Rev. D, 102, 103504
  * Gan et al. (2021) Gan F. K., Bekki K., Hashemizadeh H., 2021, preprint (arXiv:2103.09711)
  * Giovanelli & Haynes (1991) Giovanelli R., Haynes M. P., 1991, ARA&A, 29, 499
  * Glorot & Bengio (2010) Glorot X., Bengio Y., 2010, in Teh Y. W., Titterington M., eds, Proc. of Machine Learning Research Vol. 9, Proc. of the Thirteenth International Conference on Artificial Intelligence and Statistics. Journal of Machine Learning Research Inc., Brookline, Massachusetts, USA, p. 249
  * Goodfellow et al. (2020) Goodfellow I., Pouget-Abadie J., Mirza M., Xu B., Warde-Farley D., Ozair S., Courville A., Bengio Y., 2020, Communications of the Association for Computing Machinery, 63, 139
  * Gulrajani et al. (2017) Gulrajani I., Ahmed F., Arjovsky M., Dumoulin V., Courville A., 2017, in von Luxburg U., Guyon I., Bengio S., Wallach H., Fergus R., eds, Proc. of the 31st International Conference on Neural Information Processing Systems. Curran Associates Inc., Red Hook, NY, USA, p. 5769
  * Hoyle & Vogeley (2004) Hoyle F., Vogeley M. S., 2004, ApJ, 607, 751
  * Kim et al. (2011) Kim J., Park C., Rossi G., Lee S. M., Gott III J. R., 2011, Journal of Korean Astronomical Society, 44, 217
  * Kingma & Ba (2015) Kingma D. P., Ba J., 2015, in Bengio Y., LeCun Y., eds, Proc. of the 3rd International Conference on Learning Representations. International Conference on Learning Representations, San Diego, California, USA, https://arxiv.org/pdf/1412.6980.pdf
  * Kodi Ramanah et al. (2020) Kodi Ramanah D., Charnock T., Villaescusa-Navarro F., Wandelt B. D., 2020, MNRAS, 495, 4227
  * Lachieze-Rey et al. (1992) Lachieze-Rey M., da Costa L. N., Maurogordato S., 1992, ApJ, 399, 10
  * Li et al. (2012) Li B., Zhao G.-B., Koyama K., 2012, MNRAS, 421, 3481
  * Li et al. (2018) Li Y., Yu R., Shahabi C., Liu Y., 2018, in Bengio Y., LeCun Y., eds, Proc. of the 6th International Conference on Learning Representations. International Conference on Learning Representations, San Diego, California, USA, https://openreview.net/forum?id=SJiHXGWAZ
  * Mao et al. (2017) Mao Q., et al., 2017, ApJ, 835, 161
  * Mustafa et al. (2019) Mustafa M., Bard D., Bhimji W., Lukić Z., Al-Rfou R., Kratochvil J. M., 2019, Computational Astrophysics and Cosmology, 6, 1
  * Oliphant (2007) Oliphant T. E., 2007, Computing in Science & Engineering, 9, 10
  * Rodríguez et al. (2018) Rodríguez A. C., Kacprzak T., Lucchi A., Amara A., Sgier R., Fluri J., Hofmann T., Réfrégier A., 2018, Computational Astrophysics and Cosmology, 5, 4
  * Sahlén (2019) Sahlén M., 2019, Phys. Rev. D, 99, 063525
  * Sánchez et al. (2016) Sánchez C., et al., 2016, MNRAS, p. stw2745
  * Schawinski et al. (2017) Schawinski K., Zhang C., Zhang H., Fowler L., Santhanam G. K., 2017, MNRAS, 467, L110
  * Springel (2005) Springel V., 2005, MNRAS, 364, 1105
  * Springel et al. (2005) Springel V., et al., 2005, Nature, 435, 629
  * Szegedy et al. (2015) Szegedy C., et al., 2015, in Baldwin R., ed., Proc. of the 2015 Institute of Electrical and Electronics Engineers’ conference on computer vision and pattern recognition. Institute of Electrical and Electronics Engineers, Los Alamitos, California, USA, p. 1
  * Vogeley et al. (1994) Vogeley M. S., Geller M. J., Park C., Huchra J. P., 1994, AJ, 108, 745
  * White (1979) White S. D. M., 1979, MNRAS, 189, 831
  * Zingales & Waldmann (2018) Zingales T., Waldmann I. P., 2018, AJ, 156, 268

##  Appendix A Wasserstein GAN

Arjovsky et al. (2017) suggest switching the standard GAN loss function to the
Wasserstein metric (a.k.a. the Kantorovich–Rubinstein metric or Earth mover’s
distance) so that the loss function becomes:

$W(\mathbb{P}_{data},\mathbb{P}_{gen}):=\inf_{\gamma(x,y)\in\Pi(\mathbb{P}_{data},\mathbb{P}_{gen})}\mathbb{E}_{(x,y)\sim\gamma}[||x-y||]~{}~{}.$
(10)

Here, $\Pi(\mathbb{P}_{data},\mathbb{P}_{gen})$ is the set of all couplings of
$\mathbb{P}_{data}$ and $\mathbb{P}_{gen}$. If each distribution is thought to
be a unit amount of soil piled over a metric space, then $\gamma(x,y)$
represents how much soil must be moved from $x$ to $y$ such that
$\mathbb{P}_{data}$ becomes $\mathbb{P}_{gen}$. In other words, Equation 10
reveals the cost of the optimal way to transfer soil.

In practice, solving 10 is intractable. However, Arjovsky et al. (2017) have
shown that the Wasserstein metric reduces to

$W(\mathbb{P}_{data},\mathbb{P}_{gen}):=-\mathbb{E}_{\vv{\bm{z}}\sim
p_{prior}(\vv{\bm{z}})}[f(G_{\theta_{G}}(\vv{\bm{z}}))]~{}~{},$ (11)

provided that $f$ exists and is K-Lipschitz. In the WGAN algorithm, $f$ is
analogous to the set of weights, $w$, that parametrise the critic network. In
other words, given weights, $w$, lying in a compact space, $W$, one back
propagates through $\mathbb{E}_{\vv{\bm{z}}\sim
p_{prior}(\vv{\bm{z}})}[f_{w}(G_{\theta_{G}}(\vv{\bm{z}}))]$ to update the
generator’s weights, just as one would do in the standard GAN architecture.
However, to ensure the weights lie in $W$, the weights are clipped to a box
$W=[-c,c]$, where $c$ is a clipping parameter.

The training algorithm for a WGAN is similar to the standard GAN algorithm.
Aside from the different loss functions and weight clipping, the only change
in the WGAN algorithm is that the critic network has its weights updated an
integer times more than the does the generator network (i.e., in the standard
GAN architecture, both $G$ and $D$ are updated evenly). This change possible
because, unlike the Jensen-Shannon Divergence, the Earth mover’s distance is
continuous and differentiable everywhere, and allows the critic to train until
optimality is reached.

Allowing the critic to train until optimality is reached makes mode collapse
improbable in the WGAN architecture. This is due to the fact that the critic
is no longer trained on how well it can discriminate between real and fake
images. Rather, the critic is trained on how it scores the realness or
fakeness of an image. of an image. Training the critic in this way improves
network stability by giving the generator better feedback on how to update its
weights.