# Sensitivity Analysis of Simulation-Based Inference for Galaxy Clustering

Chirag Modi, Shivam Pandey, Matthew Ho, Changhoon Hahn, Bruno

Régaldo-Saint Blancard, Benjamin Wandelt

► To cite this version:

Chirag Modi, Shivam Pandey, Matthew Ho, Changhoon Hahn, Bruno Régaldo-Saint Blancard, et al. Sensitivity Analysis of Simulation-Based Inference for Galaxy Clustering. Monthly Notices of the Royal Astronomical Society, 2025, 536 (1), pp.254-265. 10.1093/mnras/stae2473. hal-04246817

HAL Id: hal-04246817

https://hal.science/hal-04246817v1

Submitted on 17 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Sensitivity analysis of simulation-based inference for galaxy clustering

Chirag Modi<sup>1,2,*</sup> Shivam Pandey,<sup>3</sup> Matthew Ho<sup>1,4</sup>, ChangHoon Hahn<sup>1,5</sup>, Bruno Régaldo-Saint-Blancard<sup>1</sup> and Benjamin Wandelt<sup>2,4</sup>

Accepted 2024 October 28. Received 2024 October 28; in original form 2023 September 29

## ABSTRACT

Simulation-based inference (SBI) is a promising approach to leverage high-fidelity cosmological simulations and extract information from the non-Gaussian, non-linear scales that cannot be modelled analytically. However, scaling SBI to the next generation of cosmological surveys faces the computational challenge of requiring a large number of accurate simulations over a wide range of cosmologies, while simultaneously encompassing large cosmological volumes at high resolution. This challenge can potentially be mitigated by balancing the accuracy and computational cost for different components of the forward model while ensuring robust inference. To guide our steps in this, we perform a sensitivity analysis of SBI for galaxy clustering on various components of the cosmological simulations: gravity model, halo finder, and the galaxy–halo distribution models (halo-occupation distribution, HOD). We infer the $\sigma_8$ and $\Omega_m$ using galaxy power spectrum multipoles and the bispectrum monopole assuming a galaxy number density expected from the luminous red galaxies observed using the Dark Energy Spectroscopy Instrument. We find that SBI is insensitive to changing gravity model between N-body simulations and particle mesh simulations. However, changing the halo finder from friends of friends to Rockstar can lead to biased estimate of $\sigma_8$ based on the bispectrum. For galaxy models, training SBI on more complex HOD leads to consistent inference for less complex HOD models, but SBI trained on simpler HOD models fails when applied to analyse data from a more complex HOD model. Based on our results, we discuss the outlook on cosmological simulations with a focus on applying SBI approaches to future galaxy surveys.

**Key words:** cosmological parameters – large-scale structure of Universe.

## 1 INTRODUCTION

The three-dimensional distribution of galaxies provides a powerful means to characterize the nature of dark matter and dark energy, to measure sum of the neutrino masses and to test gravity theory on cosmological scales. This has been the focus of various existing, ongoing, and planned galaxy redshift surveys, including Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et al. 2013), Subaru Prime Focus Spectrograph (PFS; Takada et al. 2014; Tamura et al. 2016), Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration 2016a, b; Abareshi et al. 2022), the ESA *Euclid* satellite mission (Laureijs et al. 2011), and the NASA Nancy Grace Roman Space Telescope (Roman; Spergel et al. 2015; Wang et al. 2022). However, as galaxies are a complex and biased tracer of the underlying density field, the complicated process of galaxy formation limits the ease of extracting the cosmological information from the galaxy surveys. While the clustering amplitude of the galaxy density field can be measured to percent-level precision, it cannot

straightforwardly be related to the clustering amplitude of the matter density field. Traditional methods of cosmological analysis have also largely been based on using only two- or three-point clustering statistics and analytical models based on perturbation theory (PT; Chen, Vlah & White 2022; D'Amico et al. 2022; Philcox & Ivanov 2022). As a result, these can access only linear and quasi-linear scales and are unable to exploit the full information from galaxy redshift surveys.

Over the last few years, simulation-based inference (SBI), also called likelihood-free inference or implicit-likelihood inference, has emerged as a promising approach to overcome these limitations of traditional analysis (Alsing, Wandelt & Feeney 2018; Alsing et al. 2019; Jeffrey, Alsing & Lanusse 2021; Hahn et al. 2022). This approach uses high-fidelity cosmological simulations (or forward models<fnref n="1" />) to directly model the cosmological observables in full detail. The latest SBI methods combine these simulations with neural density estimation approaches to infer the cosmological

<footnote n="1"><sup>1</sup>Center for Computational Mathematics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA</footnote>

<footnote n="2"><sup>2</sup>Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA</footnote>

<footnote n="3"><sup>3</sup>Columbia Astrophysics Laboratory, Columbia University, 550 West 120th Street, New York, NY 10027, USA</footnote>

<footnote n="4"><sup>4</sup>CNRS & Sorbonne Université, Institut d'Astrophysique de Paris (IAP), UMR 7095, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="5"><sup>5</sup>Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA</footnote>

<footnote>*E-mail: modichirag@nyu.edu</footnote>

<footnote n="1">¹In this work, we will use ‘simulations’ and ‘forward models’ interchangeably.</footnote>
parameters efficiently. Using cosmological forward models allows us to use any higher-order summary statistics of the data such as bispectrum, wavelet scattering coefficients, $k$-nearest neighbours, or even machine-learnt optimal statistics that can be evaluated in the simulations (e.g. Banerjee & Abel 2021; Eickenberg et al. 2022; Naidoo, Massara & Lahav 2022; Valogiannis & Dvorkin 2022). It also enables us to push beyond quasi-linear scales while robustly accounting for observation systematics such as imaging, completeness, fibre collisions, etc., in our modelling (Hahn et al. 2017, 2023). Meanwhile, since we use neural density estimators, we do not need to assume a Gaussian distribution for the data likelihood but can instead learn the target distributions from the simulations themselves (Hahn et al. 2019). We refer the readers to Cranmer, Brehmer & Louppe (2020) for a review on SBI. This method has also recently been applied to analyse survey data for weak lensing in Jeffrey et al. (2021) and galaxy clustering data (Hahn et al. 2023).

However, scaling SBI approaches to the next generation of surveys is not straightforward. SBI uses numerical simulations to build a model for analysing data. Thus, the accuracy and robustness of inference with SBI depends to a large extent on (i) the accuracy of the simulators and (ii) the number of simulations used to train the SBI procedure. Accounting for both these criterion simultaneously can be challenging. If the underlying simulator does not accurately model the observed data, then the inference is not reliable (Cannon, Ward & Schmon 2022). This is known as model misspecification, and the only way to safeguard against it is by using the most accurate simulations for analysis. However, this makes these simulations increasingly computationally expensive and hence for a fixed computational budget, there is a trade-off between the accuracy and the number of these simulations. This challenge is further exacerbated with the increasing volumes of cosmological surveys, and probing observables like emission-line galaxies that increasingly reside in lower mass haloes, thus requiring higher resolution simulations. Both of these factors make the simulations more expensive for a given accuracy threshold. To put things in context, the largest simulation suite currently available for training SBI for galaxy clustering (Quijote simulations; Villaescusa-Navarro et al. 2020; Hahn et al. 2023) consists only of $1 (Gpc h^{-1})^3$ in volume, which is smaller than the SDSS-III BOSS survey, and has coarse resolution of $1 Mpc h^{-1}$. Given the current status, scaling SBI approaches to the scale and fidelity required in the future can be computationally prohibitive and requires strategic planning.

**Motivation:** We take first steps towards investigating the simulation requirements for scaling SBI approaches to the next generation of galaxy clustering surveys, and study the sensitivity of SBI to the different components of the forward models used in cosmological simulations. Our goal is to ensure the robustness of inference while balancing the component models to potentially ease the computational requirements. This is motivated by the following observation: the different stages (component models) of simulations have very different computational cost and accuracy. Specifically, for dark-matter only simulations for galaxy clustering, there are three stages in the forward model: (i) evolution of dark matter under gravity, (ii) finding dark matter haloes, and (iii) populating these haloes with observed galaxies. The gravity evolution is the most computationally expensive part of the simulation, but we are also the most confident in our understanding of the underlying physics. On the other hand, we are the most uncertain about the halo-galaxy connection models, having to infer and marginalize over its parameters during the analysis. This interplay leads us to ask the question, do we need the most accurate models of the gravity evolution if we are uncertain about other components of the model, such as how to populate galaxies in the haloes?

Does it bias our results if we do not use the most accurate model for all components of the forward models? A sensitivity analysis of SBI to the different components of the cosmological forward models will answer these questions.

Covering all aspects of this sensitivity analysis is beyond the scope of a single work as the number of cases to investigate increases combinatorially with different components of the forward model, summary statistics, and parameters considered. As a result, here we will focus only on the two traditional summary statistics of galaxy clustering, power spectrum multipoles and bispectrum, but push to smaller scales than the current PT-based analyses (Ivanov, Simonović & Zaldarriaga 2020; D'Amico et al. 2022; Philcox & Ivanov 2022). We will focus on only two cosmological parameters, $\Omega_m$ and $\sigma_8$, which are well constrained by these statistics. We will consider two-component models for each of the aforementioned three stages of these simulations, gravity evolution, halo finders, and galaxy occupation and study their impact on inference.

We begin in Section 2 by describing the different forward models we will consider for the sensitivity analysis. We describe the simulation data used for each of these models in Section 3 and outline our SBI methodology in Section 4. Finally, we present our results in Section 5 and discuss implications in Section 6.

## 2 FORWARD MODELS

In this section, we describe the different models that we will consider for each of the three stages of cosmological simulations. For every stage, we implement two different component models, a simple, often computationally cheap model, and a more complex, often computationally expensive model. Our end-to-end simulations will then consist of all possible combinations of these component models.

### 2.1 Gravity models

The first step in a cosmological simulation is to evolve dark matter particles under gravity from their initial conditions set at earlier times, to their final distribution at the time of observations. This evolution is generally the most computationally expensive part of the simulations. Here, we will consider two different gravity simulations commonly used in cosmology.

(i) **N-body simulations:** These are the most accurate simulations to evolve cold dark matter (CDM) particles under gravity, for e.g. Springel (2005) and Garrison et al. (2021). N-body simulations accurately estimate gravitational forces for particles on all scales, including the particle-particle interactions on the smallest scales at every time-step, and the evolution is simulated with very small (often adaptive) time-stepping for many hundreds of time-steps.

We will use the QUIJOTE N-body simulations (Villaescusa-Navarro et al. 2020) which simulate $1024^3$ CDM particles in a $1 \text{ Gpc } h^{-1}$ box with TREEPMGADGET-III code, initialized at $z = 127$ using second-order Lagrangian Perturbation Theory (2LPT) and gravitationally evolved until $z = 0.5$. Each of these simulations requires approximately 5000 CPU hours.

(ii) **Particle-mesh simulations:** Particle-mesh (PM) simulations trade-off accuracy for speed as compared to the N-body simulations. These estimate the gravitational forces by interpolating CDM particles on a uniform force grid. As a result, these lose information on scales smaller than the grid resolution but are able to solve the Poisson equations using highly efficient fast Fourier transforms. Thus, these simulations are accurate only on the large scales but can be more than $100 \times$ cheaper than the N-body simulations (e.g. Tassev, Zaldarriaga & Eisenstein 2013; Feng et al. 2016). Recent
GPU implementations of PM simulations further increase these computational gains (Modi, Lanusse & Seljak 2021; Li et al. 2022).

For this work, we will use FastPM PM scheme (Feng et al. 2016). In each simulation, we evolve 1024³ CDM particles on a force grid of 2048³ for 10 time-steps, starting from $z = 10$ until $z = 0$. Each simulation required 200 CPU hours, a factor of 10 less than the Quijote simulations.

## 2.2 Halo model

The next step in cosmology simulations is to find high-density regions called dark matter haloes, where the dark matter particles have self-collapsed under gravity. These regions serve as sites for galaxy formation. In this work, we will use two halo finders commonly used in the community (Knebe et al. 2011).

(i) **Friends of friends (FoF):** FoF is a cluster-finding algorithm, where the clusters represent haloes in this context. Operationally, FoF finds the clusters in the simulation as follows – if two particles, two clusters, or a particle and a cluster are separated by a distance smaller than a pre-defined distance (linking length), then they are merged to form a bigger cluster (halo). We use the 3D FoF halo finder implemented in NBodykit (Hand et al. 2018). By default, this uses a linking-length of $0.2l_p$ where $l_p$ is the mean inter-particle distance.<fnref n="2" />

(ii) **Rockstar:** Rockstar algorithm is a more sophisticated phase-space algorithm for finding haloes. We only give an intuition of the algorithm here and refer the reader to the original paper (Behroozi, Wechsler & Wu 2013) for further details. Briefly, the Rockstar halo finder starts by identifying FoF haloes in 3D position space with a large linking length. It then iteratively refines these clusters using both the positions and velocities of individual CDM particles by pruning those which are inconsistent with expected phase space distribution. These haloes are generally considered to be more realistic than FoF haloes. Rockstar halo finder also estimates physical properties of the halo such as its spin, concentration etc., which are not estimated by FoF haloes.

## 2.3 Galaxy models

In CDM simulations, dark matter haloes need to be populated with galaxies. This is usually done with a statistical framework called the halo-occupation distribution (HOD; Berlind & Weinberg 2002; Zheng, Coil & Zehavi 2007). HOD provides a prescription for determining the number of galaxies, as well as their positions and velocities within every halo. The flexibility and accuracy of this framework relates to the number of parameters in the HOD prescription, which need to be inferred and marginalized during analysis. Other approaches to populate galaxies in CDM simulations, such as sub-halo abundance matching (SHAM) and semi-analytical models (Somerville & Davé 2015) require additional information from the simulations such as sub-halo distribution and merger trees, but this makes the forward simulations significantly more expensive. Hence, here we will focus on using only the following two HOD models.

(i) **Zheng07 model:** The standard HOD model (Zheng et al. 2007) assumes that the galaxy occupation depends only on the halo mass, $M_h$. We take $M_h$ to correspond to the group mass i.e. sum of all particles assigned to the halo by the FoF or Rockstar halo finder. This

model has five free HOD parameters that determine the number of central and satellite galaxies: ($\log M_{min}$, $\sigma_{log M}$, $\log M_0$, $\log M_1$, $\alpha$). Central galaxies are placed at the centre of the haloes and assigned the velocity same as the halo. Satellite galaxies are placed according to positions and velocities sampled from an NFW profile (Navarro, Frenk & White 1997).

(ii) **Zheng07ex model:** Our second model extends the standard HOD model by including additional parameters to model assembly, concentration, and velocity biases, leading to a total of 10 free HOD parameters (Hahn et al. 2023). These are implemented using the decorated HOD prescription of Hearin et al. (2016). The assembly bias parameters ($A_c, A_s$) modify the number of galaxies based on halo concentration. The concentration bias ($\eta_{conc}$) modifies the positions of satellite galaxies to allow deviation from the NFW profile of their haloes. Lastly, the central and satellite velocity biases ($\eta_c, \eta_s$) re-scale the velocities of central and satellite galaxies with respect to the host halo. This HOD model was used for a recent analysis of a subset of BOSS galaxies in the South Galactic Cap with SBI in Hahn et al. (2022).

## 2.4 End-to-end forward models

We combine the aforementioned components of our simulations in all possible combinations to generate simulations with different end-to-end forward models to train SBI procedure. However, there are two caveats:

(1) Given the two gravity, halo-finding, and HOD models each, we can have a maximum of eight Latin-hypercubes (LHs) with different forward models. However, in practice, we use only six of these as the Rockstar halo finder is not compatible with the PM simulations in its default settings. Due to the missing small-scale forces in PM simulations, the CDM particles are less clustered in phase space and Rockstar with default configuration aggressively prunes these particles resulting in inaccurate halo mass function and clustering. While it may be possible to overcome this by modifying Rockstar, it is out of scope for this work.

(2) FoF halo finder does not estimate halo concentration accurately. Thus in our FoF catalogues, it is instead estimated using analytical mass-concentration formulas from Dutton & Macciò (2014). As a result, in the Zheng07ex model, the assembly bias parameter does not capture bias based on halo assembly but instead only results in a different dependence on halo mass than is included in the standard Zheng07 HOD model. However, this caveat should not affect our conclusions.

# 3 DATA

In this section, we combine the component models described in the previous section to generate training data sets for SBI.

## 3.1 Simulations

Our simulated data consist of galaxy catalogues in redshift space at $z = 0.5$. The average number density of galaxies is $\bar{n} = 4 \times 10^{-4}$ ($h \, \text{Mpc}^{-1}$)³ with an average satellite fraction of 20 per cent. We expect similar level of co-moving galaxy number density from the luminous red galaxies (LRG) observed using the DESI survey (Zhou et al. 2023), though our estimate of satellite fraction is approximately 5–10 per cent higher compared to expectations from DESI LRGs (Berti, Dawson & Dominguez 2023; Yuan et al. 2023) In SBI, we need a training data set to learn the relationship between the observed data and underlying cosmology parameters over a wide range. Thus,

<footnote n="2">²In 3D FoF, all the distances are measured only in the three-dimensional position space as opposed to a 6D phase space.</footnote>
for each of the six composite forward models described above, we generate mock galaxy catalogues on an LH of cosmologies.

For the N-body simulations, we use the publicly available Quijote LH subset (Villaescusa-Navarro et al. 2020). It consists of 2000 simulations varying five cosmology parameters over the prior range:

$$ \Omega_m \sim \mathcal{U}[0.1, 0.5], \quad \sigma_8 \sim \mathcal{U}[0.6, 1.0], \quad \Omega_b \sim \mathcal{U}[0.03, 0.07], \\ n_s \sim \mathcal{U}[0.8, 1.2], \quad h \sim \mathcal{U}[0.5, 0.9] $$

For exact comparison, we generated the PM simulations using the same cosmological parameters and the Gaussian initial density field as of Quijote LH. In both cases, we use 1500 of these simulations for training, 200 for validation, and 300 for testing SBI.

Next, we find haloes in these simulations. For the N-body simulations, we use both Rockstar and FoF. For the PM simulations, we only use FoF for the reasons explained in Section 2.4.

Finally, we populate each of these three cases, we populate the halo catalogues with galaxies using the two HOD models described above. For each halo catalogue, we sample 20 different HOD parameter values, resulting in a total of 40 000 galaxy catalogues per forward model. Seven of these HOD parameters are sampled from the following fixed priors to be consistent with previous SBI analysis for galaxy clustering (Hahn et al. 2023)

$$ \alpha \sim \mathcal{U}[0.4, 1.0], \quad \sigma_{\log M} \sim \mathcal{U}[0.3, 0.5], \\ \eta_{\text{conc}} \sim \mathcal{U}[0.2, 2.0], \quad \eta_c \sim \mathcal{U}[0., 0.7], \quad \eta_s \sim \mathcal{U}[0.2, 2.0], \\ A_c \sim \mathcal{N}(0, 0.2) \text{ over } [-1, 1], \quad A_s \sim \mathcal{N}(0, 0.2) \text{ over } [-1, 1]. $$

For the three mass-based HOD parameters, we define priors that vary with cosmology ($\theta$) as follows:

$$ \log M_{\text{min}} \sim \mathcal{U}[\log M_{\text{min}}^{\theta} \pm 0.15], \\ \log M_0 \sim \mathcal{U}[\log M_0^{\theta} \pm 0.2], \\ \log M_1 \sim \mathcal{U}[\log M_1^{\theta} \pm 0.3]. $$

For each cosmology, $M_{\text{min}}^{\theta}$, $M_0^{\theta}$, and $M_1^{\theta}$ are set to ensure that the number density of generated galaxy catalogues is close to the target number density of $\bar{n} = 4 \times 10^{-4}$. This increases sample efficiency over using the same priors for all the cosmologies, which will need to be quite broad. We estimate $M_{\text{min}}^{\theta}$, $M_0^{\theta}$, and $M_1^{\theta}$ as follows: given the target number density $\bar{n}$ and average satellite fraction of 0.2, we estimate the average number of centrals $\bar{N}_{\text{cen}}$. For every cosmology, we use this to determine the halo mass $M_h$ above which the number of haloes is the same as $\bar{N}_{\text{cen}}$ and set $\log M_{\text{min}}^{\theta} = \log M_0^{\theta} = M_h$. With this, we then set $\log M_1^{\theta}$ to match the average number of satellites assuming a fiducial value of $\alpha = 0.7$.

## 3.2 Summary statistics

In this work, we restrict ourselves to analysing only the power spectrum multipoles $P_{\ell}(k)$ for ($\ell = 0, 2, 4$) and bispectrum monopole $B_0(k_1, k_2, k_3)$. The power spectrum multipoles are measured with fast Fourier transforms using Nbodykit (Hand et al. 2018) on a 512³ mesh. These multipoles are measured in the range $k \in [0.007, 0.5] h \text{ Mpc}^{-1}$, in bins of width $\Delta k = 2\pi/1000 h \text{ Mpc}^{-1}$. This leads to a data vector of $79 \times 3$ power spectrum coefficients. During training and testing, we also add to the power spectrum monopole a randomly sampled shot-noise contribution beyond the Poisson shot noise $S_n \sim \mathcal{U}[10^3, 10^4]$, and marginalize over it during inference. This is done to be consistent with previous $P_{\ell}(k)$ analyses (Beutler et al. 2017; Ivanov et al. 2020; Kobayashi et al. 2021; Hahn et al. 2023). The shot noise range is chosen to be comparable, but not dominant, over the expected shot noise for our number density of

tracers. However, we found that our conclusions remain the same without it.

Bispectrum is measured on a 360³ mesh using the PYSPECTRUMPYTHON package,<fnref n="3" /> which implements the Scoccimarro (2015) redshift-space bispectrum estimator. We measure bispectrum in triangle configurations defined by $k_1, k_2, k_3$ bins of width $\Delta k = 3k_f$, where $k_f = 2\pi/(1000 h^{-1} \text{ Mpc})$ is the fundamental mode. We impose the same scale cut of $k_{\text{max}} = 0.5 h \text{ Mpc}^{-1}$ as power spectrum, and this leaves us with 1980 triangle configurations.

We compare the summary statistics of our galaxy catalogues for different forward models in Fig. 1. In each column, we vary one component of the simulation at a time and show the ratio of the three summary statistics – monopole, quadrupole, and bispectrum (rows) – for the two different models considered for each component. For consistency, all the lines of the same colour have same HOD parameters (except Zheng07 model does not include the five assembly bias parameters of the extended model). The largest difference is caused by varying the HOD model between the 5- and 10-parameter models. However even with the same HOD model and parameters, changing gravity models and halo finder can lead to 10–20 per cent difference in quadrupole and bispectrum.

## 4 SIMULATION-BASED INFERENCE

Next, we outline the details of our SBI pipeline using the Latin-hypercubes generated in the previous section as the training data sets.

**Methodology:** We have generated a training data set of ($\theta, \mathbf{x}$) pairs where $\theta$ denotes the cosmology and HOD parameters, and $\mathbf{x}$ denotes the corresponding observations, i.e. the power spectrum multipoles and bispectrum. To infer the posterior $p(\theta|\mathbf{x})$, we train a conditional neural density estimator $q_{\phi}(\theta|\mathbf{x})$ with parameters $\phi$ that are fitted by maximizing the log-probability of the model parameters conditioned on the data over this training data set.

**Implementation:** We use the SNPE-C algorithm implemented in sbi<fnref n="4" /> package to train masked autoregressive flows (MAF; Papamakarios, Pavlakou & Murray 2017) as conditional neural density estimators and learn the posterior $q_{\phi}(\theta|\mathbf{x}) \sim p(\theta|\mathbf{x})$. For robustness, we train 400 networks for each data statistic by varying hyperparameters corresponding to the width and the number of layers in a single MAF block, number of MAF blocks, learning rate, and the batch size. We use Weights-and-Biases<fnref n="5" /> package for this hyperparameter exploration. After training, we collect 10 neural density estimators with best validation loss and use them as an ensemble, i.e. we construct a mixture distribution with uniform weighting to approximate the posterior. For posterior inference over a test observation $\mathbf{x}'$, we query the trained ensemble estimator $q_{\phi^*}(\theta)$ to generate samples from the posterior, i.e. $\theta \sim q_{\phi^*}(\theta|\mathbf{x}')$.

**Validation:** To validate that our posteriors are well specified, we use our trained ensemble to predict the cosmology parameters over the held-out test-data set from the same forward model as was used for training the ensemble. We use these to do calibration tests as described in Talts et al. (2020) and Hahn et al. (2023). Specifically, for each test simulation, we generate 200 samples from the posterior and find the rank of the true cosmology parameters of that particular simulation amongst these posterior samples. For correct posteriors,

<footnote n="3">³https://github.com/changhoonhahn/pySpectrum</footnote>

<footnote n="4">⁴https://github.com/mackelab/sbi</footnote>

<footnote n="5">⁵https://wandb.ai/site</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** *Comparison of summary statistics for different forward models:* We show the ratio of summary statistics for galaxy catalogues generated by varying one stage of the forward model, as indicated by the title of columns, while keeping the other two stages fixed. The three rows show the ratios for power spectrum monopole (top), quadrapole (middle) and bispectrum (showing equilateral configuration only for clarity, bottom), respectively. The three colours show three different HOD realizations (different parameter values) for the same cosmology. HOD parameters are kept consistent across the columns. The first column shows the ratio for FastPM and *N*-body simulations (with FoF halo finder and 10-parameter Zheng07-ex HOD model), the second column for simulations with FoF and Rockstar (with *N*-body gravity and Zheng07-ex HOD model), and the third column varies HOD model between 5-parameter Zheng07 and 10-parameter Zheng07-ex model (for *N*-body simulation with Rockstar halo finder).</figcaption></figure>

the distribution of these ranks should follow a uniform distribution (Talts et al. 2020). Thus, for normalized ranks (normalized between 0 and 1), the empirical cumulative distribution (CDF) corresponds to diagonal line with unit slope. The region below this line corresponds to the posteriors being overconfident and region above the line corresponds to posteriors being underconfident. We will show these rank statistic plots in the next section. This calibration is similar to checking the coverage of a credible interval under the assumed model (Hermans et al. 2021; Lemos et al. 2023). Note that this is a necessary but not a sufficient test to ensure that the posteriors are well calibrated. Furthermore since we use the same forward model for training and testing the SBI procedure in this validation, note that this does not test for model-misspecification.

# 5 RESULTS

We now perform the sensitivity analysis of SBI by looking at the impact of using different component models in training and testing the SBI procedure.

**Set-up:** We have generated mock data from six different forward models. We will use these to vary one of the three components (gravity model, halo finder, and HOD model) at a time between the two choices that are described in Section 2, while keeping the other two components fixed. In each case, we will consider inference in the two scenarios – when the test data are generated from the same forward model as the training data set, and when the test data are generated from another forward model that varies one of the three components. The first scenario validate that our
SBI procedure has been trained properly and our posteriors are well calibrated, while the second scenario gauges the impact of model misspecification.

In all cases, we infer the five cosmology and all HOD parameters using power spectrum multipoles and bispectrum. However for the sake of clarity, we present the results only for $\Omega_m$ and $\sigma_8$ that are the two parameters best constrained by these statistics. We present our results in the form of residuals, i.e. the difference between the true and the inferred mean estimate of the parameters over the held out test-data set, as well one standard deviation of uncertainties in the posterior. Additionally, we also show the rank statistic plot to verify if the posteriors are well calibrated, when relevant. In all the figures, we will use blue (and orange) colour to show the results for the case when SBI is trained and tested on the same (and different) forward model.

## 5.1 Gravity models

We begin by investigating the impact of varying gravity model between the *N*-body and PM simulations. The halo finder is fixed to FoF since, as discussed earlier, Rockstar halo finder is incompatible with PM simulations. The HOD model is fixed to 10-parameter Zheng07-ex model.

In Fig. 2(a), we show the residuals for SBI trained on both the gravity models when the true data are generated from the *N*-body simulations. For both the summary statistics (rows) and parameters (columns), the residuals are consistent, *indicating that we are not sensitive to model misspecification in this case*. This suggests that marginalizing over the HOD parameters due to the uncertainty in galaxy models indeed outweighs the refinements that happen at small scales with using more accurate gravity models. We note that there is a slight negative slope in the $\sigma_8$ residuals with power spectrum. This effect is consistent with the bounded prior on $\sigma_8$, and would likely go away with a broader prior (relative to the constraint level). However since the same trends exist in both the FastPM and Quijote posteriors, ensuring that the predictive posteriors are consistent, our conclusions regarding model misspecification still hold.

In Fig. 2(b), we show the rank statistic plots indicating that all the posteriors are also well-calibrated and do not underestimate or overestimate the posterior widths. Though not shown here, we have checked for consistency that same conclusions hold when the test observations are generated from PM simulations instead of *N*-body simulations, other components kept the same. Overall, these results are promising as they indicate that at least for this particular experimental setting, one could generate cheaper training data from PM simulation to infer parameters for the mock data generated from the expensive *N*-body simulations.

## 5.2 Halo finders

Next, we vary the halo finder in the simulations between FoF and Rockstar. The gravity model is fixed to *N*-body simulation and the HOD model is fixed to 10-parameter Zheng07-extended model.

Figs 3(a) and (b) show the residuals and rank statistic plots for SBI trained on the two halo finders and applied to test-data generated from the Rockstar halo finder. In all cases considered, the posterior for $\Omega_m$ seems to be well-calibrated and unbiased. For $\sigma_8$, the posteriors are unbiased when the summary statistic is power spectrum. However when we use bispectrum, SBI trained on Rockstar haloes infers well-calibrated posteriors for Rockstar data, but the SBI trained on FoF haloes consistently underpredicts $\sigma_8$. While not shown here, we observe similar results when the test-data are generated from FoF

catalogues: all $\Omega_m$ posteriors and $\sigma_8$ posteriors inferred from power spectrum are well calibrated, but $\sigma_8$ inferred with bispectrum from SBI trained on Rockstar catalogues is consistently biased high.

Together, these results clearly indicate that bispectrum statistic is sensitive to differences in halo finder when inferring $\sigma_8$ and SBI suffers from model misspecification. We note that similar analyses were conducted in the robustness tests of SIMBIG (Hahn et al. 2023). The test sets Test I and Test II of SIMBIG were designed to assess the sensitivity of an SBI model trained with Rockstar to the choice of the halo finder (FoF and CompaSO). However, a direct comparison is not possible since other components of the forward models were varied simultaneously (for Test I, the HOD model was also changed to the five-parameter Zheng07 HOD model, while for Test II the gravity model was changed to Abacus). These tests were also done only on a single cosmology. Despite this, similar robustness issues were observed for wavelet scattering statistics in Régaldo-Saint Blancard et al. (2023) which forced them to use aggressive scale-cuts to mitigate model misspecification.

## 5.3 Galaxy models

Finally, we change the galaxy occupation model for the simulations between the five-parameter Zheng07 and 10-parameter Zheng07-extended HOD models. The gravity model is fixed to *N*-body and we use Rockstar halo finder.

We begin by considering the test-data generated from five-parameter HOD model in Figs 4(a) and (b). SBI trained on both the HOD models gives consistent inference for both the parameters and using either of the summary statistics. This is not completely surprising given that the five-parameter HOD model is a subset of the 10-parameter HOD model, it can simply be generated by setting the assembly bias, concentration and velocity bias parameters to zero.

We turn to the more interesting case in Figs 5(a) and (b) where the test-data are generated from 10-parameter HOD model. In this case, SBI trained on the correct forward model results in well-calibrated posteriors for both the parameters from both summary statistics. However for SBI trained on the five-parameter HOD, both power spectrum and bispectrum suffer from model misspecification albeit to different degree. While posterior inferred by power spectrum is still sometimes consistent with the truth, bispectrum almost always leads to incorrect posteriors for both the parameters. This suggests that when trained on a simplistic galaxy occupation model, SBI struggles in doing inference with more complex galaxy models and this is aggravated as the summary statistics used become more informative.

Based on the results of this and the previous section, it is clear that access to accurate galaxy models will likely be the limiting factor in moving forward with all the methods that try to construct models for small scales using cosmological simulations (for e.g. SBI, machine-learning and emulator-based approaches; Yuan et al. 2022).

# 6 DISCUSSION AND OUTLOOK

We have taken the first steps towards a sensitivity analysis of SBI for galaxy clustering to answer the question – how sensitive are we to different components of our simulations? Studies like this are necessary to scale SBI approaches for the future cosmological surveys, especially as these surveys increase in volume and require higher resolution simulations to model observables. It is becoming increasingly urgent to consider the trade-offs between accuracy and
<figure><img src="image_4.png" /><figcaption>(a) Residuals and $1-\sigma$ posterior width in inferring $\Omega_m$ and $\sigma_8$ (columns) for 100 different simulations with the power spectrum (top row) and bispectrum (bottom row) statistic using SBI. The labels indicate the forward model used for training. In blue we show results for SBI trained on the correct forward model and in orange with the alternate forward model. The test-data is generated from Quijote simulations.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Rank statistics corresponding to the posteriors for which residuals are shown in Fig. 2a. We use the same color scheme, with labels indicating the forward model used for training. The test-data is generated from Quijote simulations. Power spectrum and bispectrum results are in solid and dashed lines respectively. Two columns show the two parameters. Diagonal lines following $x=y$ indicate well-calibrated posteriors.</figcaption></figure>

**Figure 2.** *Gravity models:* Varying gravity models for training SBI between the correct forward model (N-body Quijote, in blue) and the alternate forward model (FastPM simulations, in orange). In all cases, the test data are generated from the Quijote N-body simulations. The halo finder is fixed to FoF and the galaxy model is 10-parameter Zheng07-extended HOD.

the number of simulations that can be run to generate training datasets.

In this work, we have considered the problem of constraining $\sigma_8$ and $\Omega_m$ from galaxy catalogue using power spectrum and bispectrum statistics. We have varied three components of the forward simula-

tions – gravity evolution, halo finders, and galaxy occupation, and investigated their impact on inference. We find that inference in the current set-up is not sensitive to changing the gravity model between N-body and particle mesh (PM) simulations. However surprisingly, changing the halo finder between FoF and Rockstar leads to biased
<figure><img src="image_6.png" /><figcaption>(a) Same as Fig.2a, but for varying halo-finders with the test-data generated from Rockstar halos.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) Same as Fig.2b, but for varying halo-finders with the test-data generated from Rockstar halos.</figcaption></figure>

<caption>**Figure 3.** *Halo finders:* Varying halo finders for training SBI between Rockstar (blue) and FoF (orange). In all cases, the test-data are generated from the Rockstar halo finder. The gravity model is fixed to *N*-body and the galaxy model is 10-parameter Zheng07-extended HOD.</caption>

estimate of $\sigma_8$ with bispectrum. For varying galaxy models, SBI results in consistent inference when trained on a 10-parameter HOD model and tested on five-parameter HOD model, but not the other way round. When trained on 5-parameter HOD and tested on the 10-parameter model, both power spectrum and bispectrum can lead to biased results but the degree of bias for bispectrum is much larger than power spectrum.

We summarize these findings below with discussions on a more general outlook for SBI in large-scale structures.

(i) We have demonstrated that with 2000 cosmology simulations, carefully sampling HOD parameters to maximize sampling

efficiency, and properly combining neural density estimators into ensembles, we are able to obtain well-calibrated posteriors for galaxy clustering analysis with SBI. Hence for most intents and purposes, we are not limited anymore by methodological challenges in using SBI for cosmological parameter inference, at least for the realistic configurations discussed in this work. Moving forward, the primary factor in driving the quality of SBI will be the forward models used for the simulations.

(ii) We find that using PM simulations instead of *N*-body simulations does not lead to any biases in inference with power spectrum and bispectrum when combined with FoF haloes and HOD with assembly bias. Taken on its own, this has the potential of making the
<figure><img src="image_8.png" /><figcaption>(a) Same as Fig.2a, but for varying halo-galaxy occupation models with the test-data generated from 5-parameter HOD.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) Same as Fig.2b, but for varying halo-galaxy occupation models with the test-data generated from 5-parameter HOD.</figcaption></figure>

**Figure 4.** *Galaxy Model I:* Varying the halo-galaxy occupation model between the 5-parameter (blue) and 10-parameter HOD model (orange). In all cases, the test-data are generated from the five-parameter HOD. The gravity model is fixed to *N*-body and we use Rockstar haloes.

computational cost of SBI comparable to traditional analyses where similar number of simulations are required to estimate the covariance matrix (Beutler et al. 2017).

(iii) However as we move towards more powerful statistics like bispectrum, wavelet coefficients, learnt neural summary statistics, etc. to extract more information in cosmology, we become increasingly more sensitive to model misspecification in our simulators. For instance as shown in examples in Section 5, bispectrum can lead to biased results under model misspecification when power spectrum does not. Hence, we argue that robust-

ness of our inference is increasingly becoming a more challenging problem than developing summary statistics for optimal inference.

(iv) While we have focused on SBI as a specific tool for inference, the challenge of robustness is faced by all methods that use simulations for building a data model (i.e. most machine-learning or emulator-based frameworks; Yuan et al. 2022; Moran et al. 2023; Nygaard et al. 2023; Conceição et al. 2024) on small scales where simulations can be unreliable. Since SBI learns the full likelihood (or the posterior) distribution of the data, it is
<figure><img src="image_10.png" /><figcaption>(a) Same as Fig.2a, but for varying halo-galaxy occupation models with the test-data generated from 10-parameter HOD.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) Same as Fig.2b, but for varying halo-galaxy occupation models with the test-data generated from 10-parameter HOD.</figcaption></figure>

<caption>**Figure 5.** *Galaxy Model II:* Varying the galaxy model between the 10-parameter (blue) and 5-parameter HOD model (orange). In all cases, the test-data are generated from 10-parameter HOD. The gravity model is fixed to N-body and we use Rockstar haloes.</caption>

simply more suited to highlight these issues than the approaches that learn only the mean prediction and assume a Gaussian likelihood.

(v) To do sensitivity analysis of SBI, it is important to consider the end-to-end simulations rather than separately gauging the accuracy of every component. For instance, changing the halo finder from FoF to Rockstar can cause up to ~20 per cent bias in both, quadrapole and bispectrum statistics. However, marginalizing over the HOD parameters results in consistent posteriors for the former, while the HOD parametrization is not flexible enough to do the same for the latter.

(vi) Testing pipelines end to end can also lead to surprising results, for instance finding that our inference is not sensitive to the gravity model but is sensitive to the halo finder.<fnref n="6" />

(vii) This also serves to guide the new methodologies being developed to accelerate forward simulations (Dai & Seljak 2021; Jamieson et al. 2022; Lanzieri, Lanusse & Starck 2022), i.e. while it is important to report the accuracy of the simulated summary statistics, it is non-trivial to translate these to the expected results of doing inference using these accelerated simulations.

<footnote n="6">⁶Which was surprising at least for the authors.</footnote>
(viii) SBI for galaxy clustering is the most sensitive to the galaxy models. Hence for robustness, while one can train SBI on the most flexible HOD parametrization (Hahn et al. 2023), we still lack validation data for making sure that our inference is not susceptible to model misspecification. This cannot be done on models with less complex HOD parametrization. To build confidence, we need access to simulations that can accommodate different halo- galaxy occupation models such as complex HOD parametriza- tions, subhalo-abundance matching and semi-analytical models (Wechsler & Tinker 2018; Contreras, Angulo & Zennaro 2021; Yuan et al. 2022; Modi & Philcox 2023; Nguyen et al. 2023), both for training and validating our inference on the scales of future surveys.

(ix) Finally, we note that we have performed a sensitivity analysis only for two summary statistics (power spectrum multipoles and bispectrum) in inferring only two cosmology parameters (Ω_m and σ_8). The results here cannot be directly translated for other statistics and parameters (for instance, there are configurations when Ω_m is well-constrained and unbiased even when σ_8 is not). However, we argue that such a sensitivity analysis should be performed for any SBI-based data analysis to ensure that our inference is reliable. In the same vein, our work also motivates further research to develop filters beyond simple scale-cuts to make SBI analyses with higher order statistics more robust to model misspecification. Our approach provides a straightforward template to study this.

ACKNOWLEDGEMENTS

FastPM simulations were run on the KNL nodes on the Cori
supercomputer at National Energy Research Scientific Computing
Center (NERSC). MH and SP were supported by the Simons
Collaboration on Learning the Universe. The Center for Compu-
tational Astrophysics at the Flatiron Institute is supported by the
Simons Foundation. We would like to thank the Implicit Likelihood
working group of Learning the Universe collaboration for useful
discussions.

DATA AVAILABILITY

Quijote data are publicly available athttps://qui jote-simulations.
readthedocs.io/en/latest/. We plan on making FastPM simulations
online through the same channel. Access to summary statistics used
in this work can be requested by reaching out to any of the authors.
The code used for analysis is available at https://github.com/modic
hirag/contrastive cosmology/tree/main.

REFERENCES

Abareshi B. et al., 2022, preprint (arXiv:2205.10939)

Alsing J., Wandelt B., Feeney S., 2018, preprint (arXiv:1801.01497 [astro-ph])

Alsing J., Charnock T., Feeney S., Wandelt B., 2019, MNRAS, 488, 4440

Banerjee A., Abel T., 2021, MNRAS, 500, 5479

Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013, ApJ, 762, 109

Berlind A. A., Weinberg D. H., 2002, ApJ, 575, 587

Berti A. M., Dawson K. S., Dominguez W., 2023, preprint (arXiv:2303.16096)

Beutler F. et al., 2017, MNRAS, 466, 2242

Cannon P., Ward D., Schmon S. M., 2022, preprint (arXiv:2209.01845)

Chen S.-F., Vlah Z., White M., 2022, J. Cosmol. Astropart. Phys., 2022, 008

Conceição M., Krone-Martins A., da Silva A., Moliné Á., 2024, A&A, 681, A123

Contreras S., Angulo R. E., Zennaro M., 2021, MNRAS, 508, 175

Cranmer K., Brehmer J., Louppe G., 2020, Proc. Natl. Acad. Sci., 117, 30055

D'Amico G., Donath Y., Lewandowski M., Senatore L., Zhang P., 2022

Dai B., Seljak U., 2021, Proc. Natl. Acad. Sci., 118, e2020324118

Dawson K. S. et al., 2013, AJ, 145, 10

DESI Collaboration, 2016a, preprint (arXiv:1611.00036 [astro-ph])

DESI Collaboration, 2016b, preprint (arXiv:1611.00037 [astro-ph])

Dutton A. A., Macciò A. V., 2014, MNRAS, 441, 3359

Eickenberg M. et al., 2022, Wavelet Moments for Cosmological Parameter Estimation

Feng Y., Chu M.-Y., Seljak U., McDonald P., 2016, MNRAS, 463, 2273

Garrison L. H., Eisenstein D. J., Ferrer D., Maksimova N. A., Pinto P. A., 2021, MNRAS, 508, 575

Hahn C., Scoccimarro R., Blanton M. R., Tinker J. L., Rodríguez-Torres S. A., 2017, MNRAS, 467, 1940

Hahn C., Beutler F., Sinha M., Berlind A., Ho S., Hogg D. W., 2019, MNRAS, 485, 2956

Hahn C. et al., 2022, preprint (arXiv:2211.00723)

Hahn C. et al., 2023, J. Cosmol. Astropart. Phys., 2023, 010

Hand N., Feng Y., Beutler F., Li Y., Modi C., Seljak U., Slepian Z., 2018, AJ, 156, 160

Hearin A. P., Zentner A. R., van den Bosch F. C., Campbell D., Tollerud E., 2016, MNRAS, 460, 2552

Hermans J., Delaunoy A., Rozet F., Wehenkel A., Begy V., Louppe G., 2021, preprint (arXiv:2110.06581)

Ivanov M. M., Simonović M., Zaldarriaga M., 2020, J. Cosmol. Astropart. Phys., 2020, 042

Jamieson D., Li Y., Alves de Oliveira R., Villaescusa-Navarro F., Ho S., Spergel D. N., 2022, Field Level Neural Network Emulator for Cosmological N-body Simulations

Jeffrey N., Alsing J., Lanusse F., 2021, MNRAS, 501, 954

Knebe A. et al., 2011, MNRAS, 415, 2293

Kobayashi Y., Nishimichi T., Takada M., Miyatake H., 2021, preprint (arXiv:2110.06969 [astro-ph])

Lanzieri D., Lanusse F., Starck J.-L., 2022, in Machine Learning for Astrophysics. p. 60

Laureijs R. et al., 2011, preprint (arXiv:1110.3193)

Lemos P., Coogan A., Hezaveh Y., Perreault-Levasseur L., 2023, 40th International Conference on Machine Learning, Vol. 202, p. 19256

Li Y., Modi C., Jamieson D., Zhang Y., Lu L., Feng Y., Lanusse F., Greengard L., 2022, preprint (arXiv:2211.09815)

Modi C., Philcox O. H. E., 2023, preprint (arXiv:2309.10270)

Modi C., Lanusse F., Seljak U., 2021, Astron. Comput., 37, 100505

Moran K. R. et al., 2023, MNRAS, 520, 3443

Naidoo K., Massara E., Lahav O., 2022, MNRAS, 513, 3596

Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493

Nguyen T., Modi C., Yung L. Y. A., Somerville R. S., 2023, preprint (arXiv:2308.05145)

Nygaard A., Holm E. B., Hannestad S., Tram T., 2023, J. Cosmol. Astropart. Phys., 2023, 025

Papamakarios G., Pavlakou T., Murray I., 2017, Advances in Neural Infor-
mation Processing Systems, Vol. 30

Philcox O. H. E., Ivanov M. M., 2022, Phys. Rev. D, 105, 043517

Régaldo-Saint Blancard B. et al., 2023, preprint (arXiv:2310.15250)

Scoccimarro R., 2015, Phys. Rev. D, 92, 083532

Somerville R. S., Davé R., 2015, ARA&A, 53, 51

Spergel D. et al., 2015, Wide-Field InfrarRed Survey Telescope-Astrophysics Focused Telescope Assets, WFIRST-AFTA 2015 Report

Springel V., 2005, MNRAS, 364, 1105

Takada M. et al., 2014, PASJ, 66, R1

Talts S., Betancourt M., Simpson D., Vehtari A., Gelman A., 2020, preprint (arXiv:1804.06788 [stat])

Tamura N. et al., 2016, in Evans C. J., Simard L., Takami H., eds, Proc. SPIE Conf. Ser. Vol. 9908, Ground-Based and Airborne Instrumentation for Astronomy VI. SPIE, Bellingham, p. 99081M

Tassev S., Zaldarriaga M., Eisenstein D. J., 2013, J. Cosmol. Astropart. Phys., 6, 036

Valogiannis G., Dvorkin C., 2022, Phys. Rev. D, 105, 103534
Villaescusa-Navarro F. et al., 2020, ApJS, 250, 2

Wang Y. et al., 2022, ApJ, 928, 1

Wechsler R. H., Tinker J. L., 2018, ARA&A, 56, 435

Yuan S., Garrison L. H., Eisenstein D. J., Wechsler R. H., 2022, MNRAS, 515, 871

Yuan S. et al., 2023, preprint (arXiv:2306.06314)

Zheng Z., Coil A. L., Zehavi I., 2007, ApJ, 667, 760

Zhou R. et al., 2023, AJ, 165, 58

This paper has been typeset from a TeX/LATEX file prepared by the author.