Machine-learning Cosmology from Void Properties

In this paper, we train our models on data from GIGANTES, a void catalog suite containing over 1 billion cosmic voids created by running the void finder VIDE (Sutter et al. 2015) on the spatial distribution of halos from the QUIJOTE simulations (Villaescusa-Navarro et al. 2020). VIDE is a well-established public Voronoi-watershed void finder toolkit based on ZOBOV (Neyrinck 2008) and used in a variety of papers (see, e.g., Hamaus et al. 2015, 2016, 2017, 2022; Baldi & Villaescusa-Navarro 2016; Kreisch et al. 2022; Contarini et al. 2022; Douglass et al. 2023). It is the primary tool used to create the GIGANTES data set. While VIDE has a wide variety of features, one of the most important properties is its ability to capture shape, which has proven to be particularly effective in extracting cosmological information (e.g., compared to a spherical void finder; see Kreisch et al. 2022). Different techniques can be used to model voids (see, e.g., Baddeley et al. 2015; Okabe et al. 2000; Colberg et al. 2008; Cautun et al. 2018). Using the Voronoi tessellation and the watershed transform (Platen et al. 2007), VIDE is able to identify voids in dark matter or tracer distributions and provides void information including void position (R.A., decl., redshift, or x, y, z; depending on whether the catalog is from observations and light cones, or from a simulation box), void radius (calculated as

$\begin{eqnarray}&&R={\left(\displaystyle \frac{3}{4\pi }V\right)}^{1/3},\end{eqnarray} \tag{ 1 }$

where V is the total volume of the Voronoi cells that the void contains), void ellipticity (see Section 2.2), and void density contrast (see Section 2.2), among others.

GIGANTES is the first data set built for analyzing cosmic voids with machine-learning techniques. It provides void catalogs for over 7000 cosmologies, including fiducial and nonfiducial ones, and covers redshifts z = 0.0, 0.5, 1.0, and 2.0, both in real space and in redshift space. Cosmological parameters {Ω_m , Ω_b , h, n_s , σ₈, M_v , w} are varied throughout the simulations.

The data set for our paper is constituted by the void catalogs from the high-resolution QUIJOTE Latin hypercube simulations at z = 0.0. Thus, we have 2000 void catalogs covering a volume of $1\,{\left({h}^{-1}\mathrm{Gpc}\right)}^{3}$ each, from 2000 cosmological parameters arranged in a Latin hypercube with boundaries defined by

$\begin{eqnarray}&&{{\rm{\Omega }}}_{m}\in \,[0.1,0.5]\end{eqnarray} \tag{ 2 }$

$\begin{eqnarray}&&{{\rm{\Omega }}}_{b}\in \,[0.03,0.07]\end{eqnarray} \tag{ 3 }$

$\begin{eqnarray}&&h\in \,[0.5,0.9]\end{eqnarray} \tag{ 4 }$

$\begin{eqnarray}&&{n}_{s}\in \,[0.8,1.2]\end{eqnarray} \tag{ 5 }$

$\begin{eqnarray}&&{\sigma }_{8}\in \,[0.6,1.0].\end{eqnarray} \tag{ 6 }$

On average, each catalog contains about 11,000 voids. The large number of catalogs and voids allows us to train models to perform likelihood-free inference. In this work, we quantify how much information is embedded into (1) the distribution of void ellipticities, (2) the distribution of void density contrasts, and (3) void catalogs that contain radius, ellipticity, and density contrast for each void. We now describe how the ellipticity and density contrast are computed for each void.

The ellipticity of a void is computed by VIDE from the eigenvalues and eigenvectors of the inertia tensor (Sutter et al. 2015) and defined as

$\begin{eqnarray}&&\epsilon =1-{\left(\displaystyle \frac{{J}_{1}}{{J}_{3}}\right)}^{1/4},\end{eqnarray} \tag{ 7 }$

where J₁ and J₃ are the smallest and largest eigenvalues of the inertia tensor, which is calculated as

$\begin{eqnarray}\begin{array}{rcl}{M}_{{xx}} & = & \displaystyle \sum _{i=1}^{{N}_{p}}({y}_{i}^{2}+{z}_{i}^{2})\\ {M}_{{xy}} & = & -\displaystyle \sum _{i=1}^{{N}_{p}}{x}_{i}{y}_{i}\end{array}\end{eqnarray} \tag{ 8 }$

and complemented with cyclic permutations. Here, Npis the number of particles in the void while x_i , y_i , and z_i represent the cartesian coordinates of the particle i with respect to the void macrocenter. The macrocenter in the VIDE algorithm is defined as

$\begin{eqnarray}&&{{\boldsymbol{X}}}_{v}=\displaystyle \frac{1}{\displaystyle \sum _{i}{V}_{i}}\displaystyle \sum _{i}{{\boldsymbol{x}}}_{i}{V}_{i},\end{eqnarray} \tag{ 9 }$

where x _i and V_i are the positions and Voronoi cell volumes of each particle i in the void (Sutter et al. 2015).

The void density contrast, ⁵ Δ, is defined ⁶ as the ratio of the minimum density along the ridge of the void, Δ_ridge, versus the minimum density in the void, ${{\rm{\Delta }}}_{\min }$:

$\begin{eqnarray}&&{\rm{\Delta }}=\displaystyle \frac{{{\rm{\Delta }}}_{\mathrm{ridge}}}{{{\rm{\Delta }}}_{\min }}.\end{eqnarray} \tag{ 10 }$

We show in Figure 1 the ellipticity and density contrast as a function of void radii. The void radii are binned from 10.5 to 61.5 Mpc h⁻¹, and the error bars are the standard deviations of the mean for each bin. For comparison purposes, in Figure 1 we use the same void radius range as in the GIGANTES analysis (Kreisch et al. 2022), that is from 10.5 to 61.5 Mpc h⁻¹. As Figure 1 shows, when the void radius is under 30 Mpc h⁻¹, void ellipticity decreases as the void radius increases. In this range, smaller voids are more elliptical than larger voids. This is also consistent with voids found in the Magneticum simulation (Schuster et al. 2023) and voids in the Sloan Digital Sky Survey Data Release 9 (Sutter et al. 2014b), and with the expectation that small voids are less evolved structures (Sheth & van de Weygaert 2004). For voids with radii larger than 30 Mpc h⁻¹, we notice that an increase in radius leads to a mild increase in the mean values of ellipticity. However, this increase is accompanied by a notable expansion of the error bars, which can be attributed to the scarcity of voids within that range. The bottom panel shows that the void density contrast increases with the void radius. This is likely related to the locations of the detected voids since smaller voids are more commonly found in denser regions or by the edges of the larger voids, corresponding to local overdensities. Therefore, there might be smaller changes in matter density across smaller voids. On the other side, larger voids have usually been evolving for a longer time than small voids (Sheth & van de Weygaert 2004). Therefore, even if dark energy generically has a late-time impact, larger voids were dominated by dark energy earlier than smaller voids. Consequently, dark energy within large voids has had more time to make those voids larger and consequently emptier, causing a greater density contrast. Figures 2 and 3 show examples from our data set of voids with different ellipticities and density contrasts (respectively).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In our analysis, we do not impose any cut to the voids properties (e.g., only considering voids with ellipticities smaller than a given value) since machine-learning methods can learn to marginalize over features that are not correlated with cosmology (Villaescusa-Navarro et al. 2021).

The large GIGANTES data set enables cosmic void exploration with machine-learning methodologies. In this work, we train neural networks to perform likelihood-free inference on the value of five cosmological parameters: {Ω_m , Ω_b , h, n_s , σ₈}. We make use of two different architectures that perform likelihood-free inference but whose input is different:

• 1.  
  Fully connected layers. The inputs to these models are one-dimensional arrays. In our case, the data can be either the distribution of void ellipticities or the distribution of void density contrasts.
• 2.  
  Deep sets. The inputs to these models are void catalogs. The catalogs will contain a set of voids (each catalog can have a different number of voids) where each void will be characterized by three properties: (1) ellipticity, (2) density contrast, and (3) radius.

In Figure 4, we provide a visual explanation of these two architectures.

Zoom In Zoom Out Reset image size

The output of the two models is a one-dimensional array with 10 values (two values per parameter): five numbers for the posterior mean and five numbers for the posterior standard deviation. In order to train the networks to output these quantities, we use the loss function of moment networks (Jeffrey & Wandelt 2020). To avoid problems that arise from different terms having different amplitudes, we use the change described in Villaescusa-Navarro et al. (2022). Thus, the loss function we employ is

$\begin{eqnarray}\begin{array}{rcl}L & = & \displaystyle \frac{1}{| B| }\left[\mathrm{log}\left(\displaystyle \sum _{j\in B}{\left({\theta }_{j}-F({x}_{j})\right)}^{2}\right)\right.\\ & & +\,\left.\mathrm{log}\left(\displaystyle \sum _{j\in B}{\left({\left({\theta }_{j}-F({x}_{j})\right)}^{2}-G{\left({x}_{j}\right)}^{2}\right)}^{2}\right)\right]\,,\end{array}\end{eqnarray} \tag{ 11 }$

where F(x) and G(x) are functions that will output the marginal posterior mean and standard deviation when input x in batch B for the true value of parameters θ. The sum runs over all elements in a given batch and the functions are parameterized as fully connected layers or deep sets as described above. In the next section, we provide further details on the architecture and training procedure.

In order to evaluate the accuracy and precision of the model, we make use of four different statistics: (1) the coefficient of determination(R²), (2) the root mean squared error (RMSE), (3) the mean magnitude of the relative error (MMRE), and (4) the normalized chi-squared, defined as

$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}{\left(\displaystyle \frac{\mathrm{true}-\mathrm{mean}}{\mathrm{standard}\,\mathrm{deviation}}\right)}^{2}\,,\end{eqnarray} \tag{ 12 }$

where mean and standard deviation represent the posterior mean and standard deviation predicted by the model, corresponding to F(x) and G(x) in Equation (11). This measurement is directly related to the second term that we want to minimize in the loss function. The normalized chi-squared is useful to determine whether the error bars are accurately determined: a value larger/smaller than 1 indicates that the errors are underestimated/overestimated. On the other hand, we mostly rely on the R² values to evaluate the linear correlations between our models' predictions and the true values.

We start by training fully connected layers on ellipticity histograms that are constructed as follows. First, all the voids in a given simulation are selected. Second, the void ellipticities are assigned to 18 equally spaced bins from = 0.0 to = 1.0. Next, the number count in each bin is divided by the total number of voids in the considered simulations. In this way, we construct 2000 ellipticity histograms, one per simulation. We further preprocess the data so that each bin has a mean of 0.0 and a standard deviation of 1.0.

The architecture of the model consists of a series of blocks that contain a fully connected layer with a LeakyRelu nonlinear activation function and a dropout layer. We note that the last layer only contains a fully connected layer. To find the hyperparameters for our model, we use the TPESampler in Optuna package (Akiba et al. 2019) to perform Bayesian optimization over the following ranges of each hyperparameter:

• 1.  
  Number of hidden layers ∈ [1, 5].
• 2.  
  Number of neurons in each hidden layer ∈ [4, 1000].
• 3.  
  Learning rate ∈ [0.1, 1.0 × 10⁻⁵].
• 4.  
  Weight decay ∈ [1.0, 1.0 × 10⁻⁸].
• 5.  
  Dropout rate ∈ [0.2, 0.8].

The optimization is carried out by minimizing the validation loss for a total of 1000 trials. For each trial, the training is performed through gradient descent using the AdamW optimizer (Loshchilov & Hutter 2017) for 1000 epochs and a batch size of 128.

We find that the model trained on ellipticity histograms can only infer the value of Ω_m , and we show the results of inferring the value of this parameter from the ellipticity histograms of the test set in the top row of Figure 5. We find that the model is able to infer the value of Ω_m with an RMSE value of 0.06 and a mean relative error of 19%. We report the values of the four accuracy metrics in Table 1.

Zoom In Zoom Out Reset image size

We also train fully connected layers on density contrast histograms, which are constructed and standardized in the same way as the ellipticity histograms. They contain 18 bins equally spaced from Δ = 1.0 to Δ = 3.0. The model architecture, training procedure, and hyperparameter tuning are the same as the ones outlined above.

In this case, we find that the model is able to infer both Ω_m and σ₈. We show the results of this model in the middle row of Figure 5. For Ω_m , the accuracy and precision of the model are comparable to the one achieved by the model trained on ellipticity histograms. For σ₈, the model can infer its value with a mean relative error of 9%. From the values of the χ², we can conclude that the magnitudes of the error bars are well estimated. Looking at the value of the R², we can see that the estimate of Ω_m is more accurate than the estimate of σ₈. We report all the accuracy metrics in Table 1.

We have also trained neural networks using both ellipticity and density contrast histograms. In this case, the input to the network is a one-dimensional vector with 36 dimensions: 18 from the ellipticity histogram and 18 from the density contrast histogram. Using this setup, we find similar results as the ones we obtain using the density contrast histograms, possibly indicating that the information in the two statistics may be somewhat redundant.

We now illustrate a different architecture to extract information from the void catalogs directly. In this work, we are not quantifying the cosmological information encoded in the clustering of voids, but just on their intrinsic properties. For this reason, we make use of deep sets (Zaheer et al. 2017): an architecture that is able to work with sets of objects whose number may vary from set to set. Each object in the set has a collection of properties associated with it, p . Our goal is to train a model on void catalogs, to perform likelihood-free inference on the value of the cosmological parameters θ .

Our model will take all the voids in a given simulation as input and output the posterior mean ($\bar{{\theta }_{i}}$) and standard deviation (δ θ_i ) of the cosmological parameter i. Every void is characterized by three numbers: (1) its ellipticity, (2) its density contrast, and (3) its radius. We relate the input and the output through

$\begin{eqnarray}&&(\bar{{\theta }_{i}},\delta {\theta }_{i})=f\left(\bigoplus g({\boldsymbol{p}})\right),\end{eqnarray} \tag{ 13 }$

where f and g are neural networks (in this case fully connected layers) and $\bigoplus$ is a permutation invariant operation that applies to all the elements in the deep sets. In our case, we use four different permutation invariant operations: (1) the sum, (2) the mean, (3) the maximum, and (4) the minimum. The result of these four operations is concatenated before being passed to the last fully connected layer f.

The neural networks g and f are parameterized as three and four fully connected layers where the number of neurons in the latent space is considered as a hyperparameter. We use the LeakyRelu activation functions with slope −0.2 in all layers but the last. We train the model for 2000 epochs with a batch size of six using the same loss function in Equation (11) in order to perform likelihood-free inference.

Similar to the hyperparameter optimization process of the fully connected layers, we use Optuna to perform Bayesian optimization for the deep sets model. The ranges of the hyperparameters are defined below:

• 1.  
  Number of neurons in each hidden layer ∈ [10, 1000].
• 2.  
  Learning rate ∈ [1.0 × 10⁻⁴, 1.0 × 10⁻⁸].
• 3.  
  Weight decay ∈ [1.0 × 10⁻³, 1.0 × 10⁻¹²].

Here, we only perform 50 trials to optimize hyperparameters because of the longer training time. Optuna is run to minimize the value of the validation loss. While we train the model to predict the value of the five cosmological parameters we consider, we find that we are only able to infer the value of Ω_m , n_s , and σ₈. We show the results in Figure 5.

Our model is able to infer the value of Ω_m , n_s , and σ₈ with a mean relative error of approximately 10%, 3%, and 4%, respectively. This represents a significant improvement over the results we obtained using histograms of the void ellipticities or density contrasts, and is also reflected in the values of the R² statistics. This is of course expected, given that with a sufficiently complex model, this method will not throw away information when we input the data to the networks. ⁷ We note that for n_s and σ₈, our model achieves a slightly high value of χ². We have checked that this is mainly driven by a single void catalog that has a slightly extreme cosmology. If we had removed this catalog, the obtained χ² would be much closer to 1. Thus, we think our error bars are not strongly underestimated.

The reason for unsuccessful predictions could simply be the lack of information in the data sets we are using. However, we note that the parameter Ω_b will only change the amplitude and shape of the initial power spectrum in N-body simulations. Other works have shown that this effect is smaller than the ones induced by other parameters, such as Ω_m and σ₈ (see, e.g., Uhlemann et al. 2020).

We note that these results could be further improved if more void properties were to be used, or if the network would exploit clustering information (Makinen et al. 2022; Shao et al. 2022; Villanueva-Domingo & Villaescusa-Navarro 2022).