# Simulation-based Inference of Reionization Parameters from 3D Tomographic 21 cm Light-cone Images. II. Application of Solid Harmonic Wavelet Scattering Transform

Xiaosheng Zhao, Yi Mao, Shifan Zuo, Benjamin D. Wandelt

► To cite this version:

Xiaosheng Zhao, Yi Mao, Shifan Zuo, Benjamin D. Wandelt. Simulation-based Inference of Reionization Parameters from 3D Tomographic 21 cm Light-cone Images. II. Application of Solid Harmonic Wavelet Scattering Transform. The Astrophysical Journal, 2024, 973, 10.3847/1538-4357/ad5ff0. insu-04851530

HAL Id: insu-04851530

https://insu.hal.science/insu-04851530v1

Submitted on 20 Dec 2024

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>

# Simulation-based Inference of Reionization Parameters from 3D Tomographic 21 cm Light-cone Images. II. Application of Solid Harmonic Wavelet Scattering Transform

Xiaosheng Zhao<sup>1,2</sup>①, Yi Mao<sup>1</sup>①, Shifan Zuo<sup>3,4</sup>①, and Benjamin D. Wandelt<sup>2,5,6</sup>①

Received 2023 October 26; revised 2024 June 14; accepted 2024 July 5; published 2024 September 16

## Abstract

The information regarding how the intergalactic medium is reionized by astrophysical sources is contained in the tomographic three-dimensional 21 cm images from the epoch of reionization. In Zhao et al. ("Paper I"), we demonstrated for the first time that density estimation likelihood-free inference (DELFI) can be applied efficiently to perform a Bayesian inference of the reionization parameters from the 21 cm images. Nevertheless, the 3D image data needs to be compressed into informative summaries as the input of DELFI by, e.g., a trained 3D convolutional neural network (CNN) as in Paper I (DELFI-3D CNN). Here in this paper, we introduce an alternative data compressor, the solid harmonic wavelet scattering transform (WST), which has a similar, yet fixed (i.e., no training), architecture to CNN, but we show that this approach (i.e., solid harmonic WST with DELFI) outperforms earlier analyses based on 3D 21 cm images using DELFI-3D CNN in terms of credible regions of parameters. Realistic effects, including thermal noise and residual foreground after removal, are also applied to the mock observations from the Square Kilometre Array. We show that under the same inference strategy using DELFI, the 21 cm image analysis with solid harmonic WST outperforms the 21 cm power spectrum analysis. This research serves as a proof of concept, demonstrating the potential to harness the strengths of WST and simulation-based inference to derive insights from future 21 cm light-cone image data.

Unified Astronomy Thesaurus concepts: Reionization (1383); H I line emission (690); Astrostatistics (1882); Bayesian statistics (1900); Wavelet analysis (1918)

## 1. Introduction

The intensity mapping of the 21 cm line associated with the spin-flip transition of H I atoms is a promising probe of the epoch of reionization (EoR; Furlanetto et al. 2006). It contains information regarding when and how the intergalactic medium (IGM) was heated and reionized by the first luminous objects. Upper limits of the 21 cm power spectrum from the EoR (Paciga et al. 2013; Parsons et al. 2014; Jacobs et al. 2015; Mertens et al. 2020; Trott et al. 2020; Yoshiura et al. 2021; Abdurashidova et al. 2022) have been placed by ongoing interferometric array experiments, including the Precision Array for Probing the Epoch of Reionization (Parsons et al. 2010), the Murchison Widefield Array (Tingay et al. 2013), the Low Frequency Array (van Haarlem et al. 2013), and the Giant Meterwave Radio Telescope (Intema et al. 2017). In the foreseeable future, the first measurements of the 21 cm power spectrum from the EoR will be very likely achieved by upcoming array experiments, including the Hydrogen Epoch of Reionization Array (DeBoer et al. 2017) and the Square Kilometre Array (SKA; Mellema et al. 2013). Furthermore, the SKA will have the exciting promise of mapping the three-dimensional (3D) tomographic light-cone images of the 21 cm brightness temperature from the EoR with high sensitivity.

The 21 cm signal is non-Gaussian due to reionization patchiness. Therefore, the 3D light-cone images of the 21 cm signal contain more information than the power spectrum statistic, which is unlike the traditionally well-studied case of Gaussian signals in the cosmic microwave background (CMB). As such, it is of key importance for 21 cm observers to understand how to optimally extract the information in the 3D 21 cm images. The 3D image data needs to be compressed into informative summaries inevitably because it is technically very challenging to process the high-dimensional image data directly. For this purpose, several works have proposed to apply convolutional neural networks (CNNs) to compress the 2D 21 cm image slices (Gillet et al. 2019) or 3D 21 cm light-cone images (Zhao et al. 2022a; hereafter referred to as "Paper I"; Prelogović et al. 2022; Neutsch et al. 2022) into data summaries. However, the practical applications of CNN-based methods are generically computationally expensive both in generating a large volume of simulation data that are needed for training the networks and in the process of training and optimization itself. Even with so much "engineering," the fine-tuned networks may still be suboptimal (see, e.g., Paper I; Prelogović et al. 2022).

To mitigate these problems, it has been proposed to inject the inductive bias into CNNs with scattering transform (Mallat 2012; Allys et al. 2019; Cheng et al. 2020; Pedersen et al. 2022) and utilize the scattering transform to construct scattering or wavelet networks (Gauthier et al. 2021; Pedersen et al. 2022). The scattering transform employs the filters that have

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

<footnote n="1"><sup>1</sup>Department of Astronomy, Tsinghua University, Beijing 100084, People's Republic of China; zhaoxs18@mails.tsinghua.edu.cn, ymao@tsinghua.edu.cn</footnote>

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

<footnote n="3"><sup>3</sup>National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, People's Republic of China</footnote>

<footnote n="4"><sup>4</sup>Key Laboratory of Radio Astronomy and Technology, Chinese Academy of Sciences, A20 Datun Road, Chaoyang District, Beijing, 100101, People's Republic of China</footnote>

<footnote n="5"><sup>5</sup>Sorbonne Université, Institut Lagrange de Paris (ILP), 98 bis bd Arago, 75014 Paris, France</footnote>

<footnote n="6"><sup>6</sup>Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA</footnote>
well-behaved mathematical structures, e.g., the Morlet filters (Mallat 2012; Trott 2016), and under its unique definition, exploits the modulus nonlinearities and hierarchical structures —similar to the multilayers in the CNN—which allows it to extract information across multiscales. Compared with CNNs, the scattering transform has fixed filter parameters, and therefore do not need to be trained, which is a significant advantage over the CNNs. Recently, the scattering transform has been extended to 3D applications. For example, the harmonic-related wavelets are introduced to infer the molecular properties (Eickenberg et al. 2017, 2018), cosmological parameters (Saydjari et al. 2021; Chung 2022; Eickenberg et al. 2022; Valogiannis & Dvorkin 2022), and the CMB B-mode (Jeffrey et al. 2022). Specifically, Eickenberg et al. (2022) employ the first-order wavelet-based features and show the advantage of harmonic wavelets against the isotropic and oriented ones. In the context of extracting the information from the 3D 21 cm light cone (i.e., “light cuboids”), in this paper, we apply the solid harmonic wavelet scattering transform (WST; Eickenberg et al. 2017, 2018) to compress the 3D image data. The solid harmonic WST injects the inductive bias into 3D CNNs with both 3D solid harmonic wavelets and the scattering transform, which outputs multiple-order wavelet-based features.

Traditionally, the posterior distributions for the parameters of the reionization model (hereafter referred to as “reionization parameters”) can be inferred from the measurements of statistical observables of the 21 cm signal—i.e., informative summaries from the data scientific point of view—with the Markov Chain Monte Carlo (MCMC) analysis that uses an explicit (e.g., Gaussian) likelihood approximation (see the 21CMMC code; Greig & Mesinger 2015, 2017, 2018). However, the inference using the explicit likelihood approximation in 21 cm analysis may find itself either biased if the assumption in likelihood function is not exact or the covariance matrices in the likelihood neglect the off-diagonal elements between different wavenumbers and different redshifts, or computationally too expensive if all elements of covariance matrices are properly accounted for (Zhao et al. 2022b; hereafter referred to as Z22b; Prelogović & Mesinger 2023). To address this issue of intractable likelihood, the so-called “implicit likelihood inference” (ILI; Alsing et al. 2018, 2019; Papamakarios 2019; Cranmer et al. 2020; Tejero-Cantero et al. 2020), aka “simulation-based inference” (SBI) or “likelihood-free inference,” has been recently proposed to “learn” the density of the likelihood or posterior directly from data, using advanced methods in deep learning, e.g., conditional masked autoregressive flows (CMAFs; Papamakarios et al. 2017), which is a variant of normalizing flows (Papamakarios et al. 2021).

In Paper I, we introduced the density estimation likelihood-free inference (DELFI) to the 21 cm analysis for the first time and performed the posterior inference of reionization parameters from the 3D tomographic 21 cm light-cone images. The 3D CNNs were adopted therein to compress the 3D 21 cm images into informative summaries. In this paper, we improve the data compressor of 3D 21 cm images and replace the 3D CNNs in Paper I with the solid harmonic WST,<fnref n="7" /> but otherwise still perform the Bayesian inference of reionization parameters in the framework of DELFI. This new approach (i.e., solid harmonic WST with DELFI) is dubbed “3D ScatterNet”

herein. We will compare the inference results using 3D ScatterNet and those using DELFI-3D CNN, in order to demonstrate the improvement of the data compressor. In addition, we will compare the inference results from the 21 cm images using 3D ScatterNet and those from the 21 cm power spectrum analysis using 21cmDELFI-PS (Z22b), both in the framework of DELFI, in order to demonstrate whether solid harmonic WST can extract more information from the 21 cm images than the power spectrum statistic. Similar to the analysis in Z22b, realistic effects, including thermal noise and residual foreground after removal with the singular value decomposition (SVD; Stewart 1993; Masui et al. 2013; Wolz et al. 2015), are applied to the mock observations from the SKA.

The rest of this paper is organized as follows. The 3D ScatterNet methodology is introduced in Section 2. We describe the simulations and the application of realistic effects in Section 3, present the inference results in Section 4, and make concluding remarks in Section 5. Some technical details are left to Appendix A (on the effect of angular frequency), Appendix B (on the effect of smoothing scales), Appendix C (on the dependence of light-cone volumes), Appendix D (on the network setting and sample size), and Appendix E (on the performance of DELFI-3D CNN with SKA noise). Some of our results were previously summarized by us in a conference paper (Zhao et al. 2023b). Our code is publicly available on GitHub.

## 2. The 3D ScatterNet Methodology

### 2.1. Solid Harmonic WST

We briefly summarize the solid harmonic WST in this subsection, following Eickenberg et al. (2017, 2018).

Solid harmonic WST convolves the original fields with a cascade of solid harmonic wavelets, performs nonlinear moduli on these convolved fields, and integrates over the coordinate space. The solid harmonic wavelet is defined as the solid harmonic multiplying a Gaussian

$$ \psi_{\ell}^{m}(\mathbf{u}) = \frac{1}{(\sqrt{2\pi})^3} e^{-|u|^2/2} |u|^{\ell} Y_{\ell}^{m}\left(\frac{\mathbf{u}}{|u|}\right), \quad (1) $$

where $|u|^{\ell} Y_{\ell}^{m}\left(\frac{\mathbf{u}}{|u|}\right)$ is the solid harmonic, evaluated at the coordinate $\mathbf{u}$. The Gaussian part serves as localizing the wavelet support around zero. Here, we omit the further normalization factors for simplicity. In order to capture features of the field at multiple scales, the mother solid harmonic wavelet in Equation (1) is dilated at the scale<fnref n="8" /> $2^j$, i.e.,

$$ \psi_{j,\ell}^{m}(\mathbf{u}) = 2^{-3j} \psi_{\ell}^{m}(2^{-j}\mathbf{u}). \quad (2) $$

The Euclidean norm as the modulus operator, also dubbed “first-order modulus coefficient” herein, is defined by

$$ U[j, \ell] \mathbf{d}(\mathbf{u}) = \left( \sum_{m=-\ell}^{\ell} |\mathbf{d}^* \psi_{j,\ell}^{m}(\mathbf{u})|^2 \right)^{1/2}, \quad (3) $$

where the field $\mathbf{d}$ is convolved (denoted by “*” in Equation (3)) by the dilated solid harmonic wavelets at the scale $j$ with the angular frequency band $\ell$. The additional rotational phase

<footnote n="7">⁷ The solid harmonic WST is implemented with the Kymatio package (Andreux et al. 2018; https://www.kymatio.io).</footnote>

<footnote n="8">⁸ Without loss of clarity, we refer to the expression of “at the scale $2^j$” as “at the scale $j$” hereafter for simplicity.</footnote>
subspace information represented by *m* is aggregated to produce coefficients that are covariant to rotation. The translational covariance is also guaranteed by the convolution operation.

The first-order solid harmonic wavelet scattering coefficient is given by

$$S_1[d; j, \ell, q] = \int_{\mathbb{R}^3} |U[j, \ell]d(\mathbf{u})|^q d^3\mathbf{u}, \quad (4)$$

where the modulus is raised by the power *q*, which results in the coefficients that are sensitive to the amplitude of the field, in that the small (large) value of *q* gives more weight to the small (large) nonzero values of the integral. Note that the point-wise transformation by the power *q* does not change the covariance property. By integrating these coefficients over the position **u**, we get the coefficients that are invariant to both translation and rotation.

In order to capture the information across multiple scales, the first-order modulus coefficient $U[j, \ell]d(\mathbf{u})$ is convolved with another wavelet at a different scale $j'$ with $j' > j$ but with the same angular frequency band $\ell$, i.e., the second-order modulus coefficient is defined as

$$U[j, j', \ell]d(\mathbf{u}) = \left( \sum_{m=-\ell}^{\ell} |U[j, \ell]d^*\psi_{j,\ell}^m(\mathbf{u})|^2 \right)^{1/2}, \quad j < j'. \quad (5)$$

The covariance property is also maintained.

The second-order solid harmonic wavelet scattering coefficient is defined by integration over the coordinate space, similar to Equation (4),

$$S_2[d; j, j', \ell, q] = \int_{\mathbb{R}^3} |U[j, j', \ell]d(\mathbf{u})|^q d^3\mathbf{u}, \quad (6)$$

which is also invariant to both translation and rotation.

In principle, successive similar operations can be applied to define the higher-order solid harmonic wavelet scattering coefficients.<fnref n="9" /> In this paper, we neglect the information encoded in the higher-order scattering coefficients and only consider the zeroth-, first-, and second-order scattering coefficients, following Allys et al. (2019) and Cheng et al. (2020). In the same spirit as in Equations (4) and (6), the zeroth-order scattering coefficient is defined as the sum of all pixel values raised by the power *q*, i.e.,

$$S_0[d; q] = \int_{\mathbb{R}^3} d(\mathbf{u})^q d^3\mathbf{u}. \quad (7)$$

Basically, while the zeroth-order scattering coefficients highlight the amplitude of the field **d**, the first-order scattering coefficients separate the information of the field into different scales by gathering the specific geometric feature in the original field at the scale *j*. Also, the second-order scattering coefficients represent the information of nonlinear mode mixing between different scales *j* and *j'*. These nonlinear operations obviously extract a richer set of information than the power spectrum analysis that only characterizes the features at separate Fourier modes. In addition, not only does the WST capture the local spatial information with the localized wavelets, but it also saves some large-scale information because of the long tail of the

wavelets, similar to the Morlet transform (Trott 2016). Furthermore, the scattering coefficients are naturally invariant to translation and rotation. Also, they are Lipschitz continuous to deformations (Mallat 2012a), meaning that they are approximately proportional to small deformations of the original field.

Solid harmonic WST is very analogous to 3D CNN in three aspects as follows. First, the modulus operators work similarly to the nonlinear functions in the latter. Second, the integration of the modulus coefficients over the coordinate space is essentially the pooling operation. Third, the orders in scattering coefficients are similar to the layers in the CNN. However, solid harmonic WST has fixed kernels in wavelets, so unlike the CNN, it does not need to be trained in order to output the scattering coefficients. Also, large kernels are usually avoided in the CNN because otherwise, this would require a huge set of training parameters, but large kernels can be applied to solid harmonic WST efficiently.

Figure 1 presents the visualization of selected modulus coefficients. Each coefficient reflects some features of the original field. For example, the first-order modulus coefficients with ($j=0, \ell=1$) highlight the boundaries of H II bubbles, while the second-order modulus coefficients with ($j=0, j'=1, \ell=1$) show more complex structures. Also, as $\ell$ increases (from the left to right in Figure 1), the modulus coefficients encode the information of small structures of the original field. Note that wavelets with $\ell=0$ are Gaussian, while wavelets with nonzero $\ell$ can encode rich structures such as filaments.

To compute the scattering coefficients, we choose the maximum scale $j_{max}=5$ and maximum angular frequency $\ell_{max}=6$, and select the value of the modulus power $q=0.5$, 1, or 2 for the first and second-order scattering coefficients. The half-width parameter (i.e., the standard deviation in the Gaussian factor) of the mother solid harmonic wavelets is set to be unity, in which case the maximum scale $j_{max}=5$, which satisfies the criterion $2 \times (\text{half-width parameter}) \times 2^{j_{max}} \le$ the number of simulation cells on each side (66 in this work, as shown in Section 3.1). For a given $\ell$ (with $0 \le \ell \le \ell_{max}$) and $q$, the number of the first-order coefficients is $j_{max}+1=6$ since $j$ takes the value of $0 \le j \le j_{max}$ and the number of the second-order coefficients is $j_{max}((j_{max}+1)/2 = 15$ since $j$ and $j'$ take the value of $j=0, 1 \le j' \le j_{max}, j=1, 2 \le j' \le j_{max}$, etc. To keep the dimension of coefficients reasonably low, in this paper, we simply average the information over different $\ell$ for a given $q$. (We will discuss the effect of angular frequencies $\ell$ on the parameter inference in Appendix A.) Therefore, the first and second-order coefficients have a total of 63 components for three values of power $q=0.5, 1$, and 2.

The zeroth-order coefficient is ill-defined if $q=0.5$ because the field value at a single pixel can be negative. Also, the zeroth-order coefficient is zero if $q=1$ because it is simply the global sum of the field in this case, which is zero for the 21 cm signal by definition (see Section 3.1). Therefore, instead of $q=0.5$ or 1, we choose (arbitrarily) three higher modulus powers $q=2, 3$, and 4 in this paper for the zeroth-order coefficients.

Altogether, the final concatenated coefficient vector has a dimension of 66, which can be adjusted if necessary. Furthermore, following Allys et al. (2019), we take the logarithms of those coefficients with base 2 as our new coefficients. (If a component is negative, we take the logarithm

<footnote n="9">⁹ Without loss of clarity, we refer to the phrase “solid harmonic wavelet scattering coefficients” as scattering coefficients hereafter for simplicity.</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Visualization of solid harmonic WST. Shown is one slice of the original field of simulated cosmological 21 cm signals along the angular direction (top row), the first-order modulus coefficients $U[j, \ell]d(u)$ (middle row) and the second-order modulus coefficients $U[j, j', \ell]d(u)$ (bottom row), at some representative wavelet scales *j* and/or *j'* and at various angular frequencies *ℓ*. The color map of all fields is the same as shown for the original field.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** The first-order ("$S_1$") and second-order ("$S_2$") scattering coefficients of the simulated cosmological 21 cm signals and their dependence on the reionization parameters. The scattering coefficients are averaged over $0 \le \ell \le \ell_{max}$ and evaluated at $q = 1$. We vary two reionization parameters separately, namely, $\log_{10}(T_{vir})$ (top) and $\log_{10}\zeta$ (bottom)—see their definitions in Section 3.1—with the fiducial model in green lines.</figcaption></figure>

of the absolute value but keep the sign of the component.) The advantage of adopting these new coefficients is that they behave linearly with *j* and *j'*. Hereafter the phrase “scattering coefficients” refers to their logarithms.

We show the representative values of scattering coefficients in Figure 2. We vary the reionization parameters (see their definitions in Section 3.1), and find that the coefficients are generically suppressed (enhanced) as $\log_{10} T_{vir}$ ($\log_{10} \zeta$) increases. This implies that the scattering coefficients are sensitive to the reionization parameters.

## 2.2. SBI with CMAFs

Figure 3 shows the workflow of 3D ScatterNet. The scattering coefficients extracted by the solid harmonic WST serve as the input summaries of the CMAFs,<fnref n="10" /> which is a variant of the neural density estimators (NDEs) that perform the implicit likelihood inference. We refer interested readers to Paper I for an in-depth description of the DELFI and

<footnote n="10">¹⁰ The CMAFs are implemented with `pydelfi` (Alsing et al. 2019; https://github.com/justinalsing/pydelfi).</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** The workflow of 3D ScatterNet. (Left, fixed part) The data compression with solid harmonic WST—the 3D light-cone images *d*, which are simulated with the parameters *θ* are compressed by a cascade of scattering transforms (each containing the convolution with solid harmonic wavelets, harmonic modulus, and integration operation) to form the summaries *t*, which include the zeroth, first, and second-order scattering coefficients {S₀, S₁, S₂}, respectively. (Right, trainable part) The CMAFs are used to learn the summary density conditional on the parameters (i.e., likelihood). With the likelihood, the posterior can be inferred using Bayes' theorem at the data summary *t₀*.</figcaption></figure>

the CMAFs. The detailed settings of NDEs are given in Appendix D. In this paper, we set two neural layers of a single transform and 50 neurons per layer for CMAFs. We also use the ensembles of CMAFs to improve the performance. The number of transforms and the configuration of ensembles are chosen based on the performance of posterior validation.

## 2.3. Validation of the Posterior

As the final step of inference, we validate both marginalized and joint posteriors. In the posterior validation, hypothesis tests are made to check if the posteriors from CMAFs are self-consistent statistically. Note that the posterior validation is sometimes referred to as “posterior calibration.” However, this step is not to calibrate the trained networks, but simply an indication of whether the network complexity and training data are adequate to learn the conditional density accurately in a statistical manner. We follow Appendix A of Paper I and recap the validation statistics in this subsection.

For the marginalized posterior, the probability integral transform (PIT. Gneiting et al. 2007; Mucesh et al. 2021) is defined as

$$ \text{PIT}(\tilde{\theta}) = \int_{-\infty}^{\tilde{\theta}} f(\theta) d\theta, \qquad (8) $$

i.e., the cumulative distribution function (CDF) of the inferred marginal distribution $f(\theta)$ at the true value $\tilde{\theta}$. If the inferred posteriors are accurate, then the PIT should be uniformly distributed.

For the joint posteriors, there are two statistics. The first one is the copula probability integral transformation (copPIT; Ziegel & Gneiting 2014; Mucesh et al. 2021), which is defined as

$$ \text{copPIT} \equiv \text{Pr}\{H(\theta) \le H(\tilde{\theta})\}. \qquad (9) $$

Here, “Pr” represents the probability and $H(\theta)$ is the CDF of the inferred joint distribution. The copPIT is the multivariate extension of PIT.

The second statistic is the highest probability density (HPD. Harrison et al. 2015), which is defined as

$$ \text{HPD}(\tilde{\theta}) = \int_{f(\theta) \ge f(\tilde{\theta})} f(\theta) d^n\theta. \qquad (10) $$

The HPD describes the plausibility of $\tilde{\theta}$ under the distribution $f(\theta)$. A small value of HPD indicates high plausibility. Similar to the PIT, the copPIT and HPD should also be uniformly distributed if the posteriors are accurate.

In order to check the uniformity, we adopt two metrics, the Kolmogorov–Smirnov (KS; Kolmogorov 1992) test and Cramér–von Mises (CvM; Anderson 1962) test, which focus on different aspects of distribution. While the KS test is sensitive to the median, the CvM test captures the tails of a distribution. If the *p*-value from a test is larger than the preset value of significance level, typically 0.01 or 0.05 (Ivezić et al. 2014), then the null hypothesis that these statistics follow a uniform distribution is accepted. Throughout this paper, we adopt the significance level of 0.01 and report the results only if the *p*-value is larger than 0.01 unless stated otherwise.

# 3. Data Preparation

## 3.1. Cosmic 21 cm Signal

The 21 cm brightness temperature relative to the CMB temperature at position *x* can be written (Furlanetto et al. 2006) as

$$ T_{21}(x, z) = \tilde{T}_{21}(z)x_{\text{HI}}(x)[1 + \delta(x)]\left(1 - \frac{T_{\text{CMB}}}{T_S}\right), \qquad (11) $$

where $\tilde{T}_{21}(z) = 27\sqrt{[(1+z)/10](0.15/\Omega_m h^2)(\Omega_b h^2/0.023)}$ in units of millikelvin $x_{\text{HI}}(x)$ is the neutral fraction, and $\delta(x)$ is the matter overdensity, at position *x*. We assume the cold
dark matter can be traced by the baryon perturbation on large scales, so $\delta_{\rho_H} = \delta$. In this paper, we focus on the limit where spin temperature $T_S \gg T_{CMB}$, likely valid soon after reionization begins, though this assumption is strongly model dependent. As such, we can neglect the dependence on spin temperature. Also, as a demonstration of concept, we ignore the effect of peculiar velocity; such an effect can be readily incorporated in forward simulations by the algorithm introduced by Mao et al. (2012).

In this paper, we use the publicly available code 21cmFAST<fnref n="11" /> (Mesinger & Furlanetto 2007; Mesinger et al. 2011), which can be used to perform seminumerical simulations of reionization, as the simulator to generate the data sets. Our simulations were performed on a cubic box of 100 comoving Mpc on each side, with 66³ grid cells. Following the interpolation approach in Paper I, the snapshots at nine different redshifts of the same simulation box (i.e., with the same initial condition) are interpolated to construct a light-cone 21 cm data cube within the comoving distance of a simulation box along the line of sight (LOS). We concatenate 10 such light-cone boxes, each simulated with different initial conditions in density fields but with the same reionization parameters, together to form a full light-cone data cube<fnref n="12" /> of the size 100 × 100 × 1000 comoving Mpc³ (or 66 × 66 × 660 grid cells) in the redshift range $7.51 \le z \le 11.67$. To mimic the observations from radio interferometers, we subtract from the light-cone field the mean of the 2D slice for each 2D slice perpendicular to the LOS, because radio interferometers cannot measure the mode with $k_\perp = 0$.

We parameterize our reionization model as follows, and refer interested readers to Paper I for a detailed explanation of their physical meanings.

(1) $\zeta$, the ionizing efficiency, which is a combination of several parameters related to ionizing photons. In our paper, we vary $\zeta$ as $10 \le \zeta \le 250$.

(2) $T_{vir}$, the minimum virial temperature of halos that host ionizing sources. In our paper, we vary this parameter as $4 \le \log_{10}(T_{vir} / K) \le 6$.

Cosmological parameters are fixed in this paper as ($\Omega_\Lambda$, $\Omega_m$, $\Omega_\rho$, $n_s$, $\sigma_8$, $h$) = (0.692, 0.308, 0.0484, 0.968, 0.815, and 0.678) (Planck Collaboration et al. 2016).

## 3.2. Thermal Noise and Residual Foreground

To apply the realistic effects to the 21 cm signal, we first compute the uv coverage map at each of a total of 660 frequency channels. The uv coverage map can be used to suppress the thermal noise and calculate the telescope response to the cosmic 21 cm signal and foreground. We then generate the thermal noise with a 1000 hr observation. For the thermal noise estimation in this paper, we employ the Tools21cm<fnref n="13" /> (Giri et al. 2020) code to simulate the expected thermal noise in the 3D 21 cm light-cone images as observed with the SKA1-Low. We list our basic assumptions of the SKA configuration in Table 1. When synthesizing the total signal, we applied the same telescope response in the calculation of thermal noise to both the cosmic 21 cm signal component and the radio foreground component. The signal of a pixel in the uv space

<table><caption>**Table 1**
Specifications for the SKA</caption><thead><tr><td>Parameters</td><td>Values</td></tr></thead><tbody><tr><td>System temperature</td><td>60($\frac{v}{300 \text{ MH}}$)<sup>-2.55</sup> + 100 K</td></tr><tr><td>Effective collecting area</td><td>962 m<sup>2</sup></td></tr><tr><td>Integration time</td><td>10 s</td></tr><tr><td>Observation hours per day</td><td>6 hr</td></tr><tr><td>Total observation time</td><td>1000 hr</td></tr></tbody></table>

is kept unchanged for nonzero uv coverage but otherwise set to zero if the pixel is outside the uv coverage. In order to suppress the noise, we smooth the map with a relatively small baseline of 1 km,<fnref n="14" /> which roughly corresponds to the size of the core area of SKA1-Low. Specifically, we convolve the map with a Gaussian filter in the angular direction and a top-hat filter in the LOS direction. The FWHM of the Gaussian filter corresponds to the size of the 1 km baseline and the width of the top-hat filter also matches that FWHM.

To mock the foreground contamination, we employ the pygsm<fnref n="15" /> package that is based on the GSM-building model (Zheng et al. 2017), which interpolates the sky maps with 29 sky map observations using an improved principal component analysis method. We then use these foreground maps to interpolate our images on the grid for each frequency slice. In order to prevent overfitting, we assign a random patch of sky (except for the north and south poles) to model each foreground image. Finally, we employ the SVD technique for foreground removal. Specifically, we remove the largest six singular value modes.

For visualization purposes, we show the mock 21 cm observation with realistic effects in Figure 4. When the cosmic 21 cm signal is applied with the thermal noise of SKA and/or the residual foreground, the features at small scales are washed out but those at the very large scales are still retained.

## 4. Results

For each approach of data compression, we first perform the inference for two representative mock observations—the “Faint Galaxies Model” and the “Bright Galaxies Model,” defined in Tables 2 and 3 (see the “True value” therein), respectively, following the convention of Greig & Mesinger (2017; see also, Paper I and Z22b). They are selected as two illustrative examples with extreme parameter values, which are nevertheless fine-tuned for these models to have similar global reionization histories. We further perform the posterior validation of the trained NDEs using another set of 300 independent samples that are randomly drawn from the allowed region in the parameter space in which the mean neutral fraction satisfies $0.08 \le \bar{x}_{HI} \le 0.81$ at $z=7.1$, corresponding to the 2$\sigma$ credible region as constrained by the IGM damping wing of ULASJ1120+0641 (Greig et al. 2017).

### 4.1. 3D ScatterNet versus 3D CNN

We compare the results using 3D ScatterNet and using DELFI-3D CNN in this subsection. For simplicity, here we

<footnote n="11">¹¹ https://github.com/andreimesinger/21cmFAST</footnote>

<footnote n="12">¹² We investigate the contributions of individual wave bands to the parameter inference in Appendix C.</footnote>

<footnote n="13">¹³ https://github.com/sambit-giri/tools21cm</footnote>

<footnote n="14">¹⁴ Giri et al. (2018), Bianco et al. (2021), Giri & Mellema (2021) choose to smooth with a baseline of 2 km, which corresponds to the size of the central area of SKA1-Low. We discuss the effect of smoothing scales in Appendix B.</footnote>

<footnote n="15">¹⁵ https://github.com/telegraphic/pygsm</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Visualization of the mock 21 cm observation in a slice along the LOS. From top to bottom: (1) the simulated cosmic 21 cm images without realistic effect applied except that the $k_{\perp} = 0$ mode is removed from each 2D slice perpendicular to the LOS (i.e., pure signal); (2) same as the pure signal but with total noise from SKA observation (i.e., SKA noise); (3) same as the pure signal but with total noise from SKA observation and the foreground contamination; (4) same as the pure signal but with total noise from SKA observation and residual foreground after the foreground is removed with the SVD technique (i.e., SKA noise + residual foreground). The map of pure signal retains the original grid size in simulations, but the rest of the three maps are made by convolving the original maps with a Gaussian filter in the angular direction and a top-hat filter in the LOS direction, with the widths of the filters corresponding to the size of the 1 km baseline.</figcaption></figure>

only consider the case of cosmological 3D 21 cm images without realistic effect applied except that the $k_{\perp} = 0$ mode is removed from each 2D slice perpendicular to the LOS (i.e., “pure signal”). Figure 5 shows the credible regions ($1\sigma$ and $2\sigma$, or 0.68 and 0.95 levels) for the Faint Galaxies Model and the Bright Galaxies Model, with quantitative comparisons in terms of medians and $1\sigma$ (16th and 84th percentile) errors presented in Tables 2 and 3.

The results from DELFI-3D CNN were taken from Paper I, in which the ILI was performed with data compression made by a trained 3D CNN. Paper I used a set of 9000 samples for training and validation of the 3D CNN and another set of 9000 samples for training and validation of the density estimators. Nevertheless, if the training sample size was doubled, the results for the Faint Galaxies Model and the Bright Galaxies Model were not improved, which implies that the inference accuracy might be limited by some intrinsic properties in experimental choices (including the network architecture, characteristics in training data sets, and other hyperparameter choices) for training the 3D CNN. However, when we replace the data compression from 3D CNN to solid harmonic WST, we find that the inference results are significantly improved in terms of the location and size of credible regions in the posterior distributions for the reionization parameters, as shown in Figure 5 and Tables 2 and 3. Also, the degeneracy of these two parameters is clearly revealed in the credible region for 3D ScatterNet, which indicates that this degeneracy is intrinsic in the theoretical modeling.

Now that the solid harmonic WST represents a better approach to data compression than 3D CNN, we will focus on the 3D ScatterNet in the remainder of this paper.

## 4.2. 3D ScatterNet versus 21 cm DELFI-PS

Z22b shows that for the power spectrum analysis, the DELFI approach (i.e., 21cmDELFI-PS) outperforms the standard

MCMC analysis. In this subsection, therefore, we compare the results using 3D ScatterNet and 21cmDELFI-PS, both under the same inference strategy using DELFI. On this same footing, we can investigate whether data compression using solid harmonic WST extracts more information from the 3D 21 cm light-cone images than the power spectrum statistic.

### 4.2.1. Pure Signal

We first consider the pure signal case for the comparison, with the credible regions for the Faint Galaxies Model and the Bright Galaxies Model shown in Figure 5 and quantitative results listed in Tables 2 and 3. For the Faint Galaxies Model, the systematic shift (i.e., relative errors of the predicted medians with respect to the true values) and the $1\sigma$ statistical errors are $0.04\% \pm 0.13\%$ ($0\% \pm 0.23\%$) for $\log_{10}(T_{\text{vir}})$ with 3D ScatterNet (21cmDELFI-PS), respectively,$^{16}$ and $0.14\% \pm 0.04\%$ ($0.27\% \pm 0.74\%$) for $\log_{10}\zeta$ with 3D ScatterNet (21cmDELFI-PS), respectively. The estimated statistical errors using 3D ScatterNet are about 1.7 times smaller than using 21cmDELFI-PS. The similar results hold generically for the Bright Galaxies Model, too.

Next, we test the trained NDEs on 300 samples. Table 4 shows the coefficient of determination, $R^2 = 1 - \frac{\sum(y_{\text{pred}} - y_{\text{true}})^2}{\sum(y_{\text{true}} - \bar{y}_{\text{true}})^2}$, where $y_{\text{true}}$ and $y_{\text{pred}}$ are the true value and the predicted median in a sample for the parameter $y$ (e.g., $\log_{10}(T_{\text{vir}})$ and $\log_{10}\zeta$ in this paper), respectively, and the summation is over all testing samples. $\bar{y}_{\text{true}}$ is the average of the true value over all testing samples. A score of $R^2$ close to unity indicates an overall good inference performance of this parameter. For the case of pure signal, both

<footnote>$^{16}$ Here for this mock observation, the predicted median with 3D Scatter-Net has slightly more deviation (but still within the $1\sigma$ credible region) from the true value than that with 21cmDELFI-PS. However, testing the medians with 300 samples, we find that the predicted medians using 3D ScatterNet are statistically closer to the true values than those using 21cmDELFI-PS.</footnote>
<table><caption>**Table 2**
Bayesian Inference with 21cmDELFI-PS, DELFI-3D CNN, and 3D ScatterNet for the Faint Galaxies Model</caption><thead><tr><th rowspan="2">Parameter</th><th rowspan="2">True Value</th><th colspan="3">Pure Signal</th><th colspan="2">SKA Noise</th><th colspan="2">SKA Noise + Residual Foreground</th></tr><tr><th>DELFI-3D CNN</th><th>21cmDELFI-PS</th><th>3D ScatterNet</th><th>21cmDELFI-PS</th><th>3D ScatterNet</th><th>21cmDELFI-PS</th><th>3D ScatterNet</th></tr></thead><tbody><tr><td>log<sub>10</sub>(T<sub>vir</sub>/K)</td><td>4.699</td><td>4.697<sup>+0.024</sup><sub>-0.024</sub></td><td>4.699<sup>+0.011</sup><sub>-0.011</sub></td><td>4.701<sup>+0.006</sup><sub>-0.006</sub></td><td>4.716<sup>+0.081</sup><sub>-0.070</sub></td><td>4.713<sup>+0.085</sup><sub>-0.059</sub></td><td>4.789<sup>+0.185</sup><sub>-0.152</sub></td><td>4.774<sup>+0.133</sup><sub>-0.120</sub></td></tr><tr><td>log<sub>10</sub>(ζ)</td><td>1.477</td><td>1.475<sup>+0.023</sup><sub>-0.023</sub></td><td>1.481<sup>+0.011</sup><sub>-0.011</sub></td><td>1.479<sup>+0.007</sup><sub>-0.006</sub></td><td>1.508<sup>+0.051</sup><sub>-0.043</sub></td><td>1.498<sup>+0.059</sup><sub>-0.038</sub></td><td>1.552<sup>+0.123</sup><sub>-0.099</sub></td><td>1.533<sup>+0.097</sup><sub>-0.077</sub></td></tr></tbody></table>

**Note.** Here, pure signal refers to the mock observations of cosmological 3D 21 cm images (i.e., without thermal noise or foreground contamination, but the $k_{\perp} = 0$ mode is removed from each 2D slice perpendicular to the LOS); SKA noise refers to the mock SKA observations of the 3D 21 cm images with total noise (with the contributions from thermal noise and sample variance errors) yet without foreground contamination; SKA noise + residual foreground refers to the mock SKA observations of the 3D 21 cm images with total noise and residual foreground after the foreground is removed with the SVD technique.
<table><caption>Table 3
Same as Table 2 but for the Bright Galaxies Model</caption><thead><tr><th rowspan="2">Parameter</th><th rowspan="2">True Value</th><th colspan="3">Pure Signal</th><th colspan="2">SKA Noise</th><th colspan="2">SKA Noise + Residual Foreground</th></tr><tr><th>DELFI-3D CNN</th><th>21cmDELFI-PS</th><th>3D ScatterNet</th><th>21cmDELFI-PS</th><th>3D ScatterNet</th><th>21cmDELFI-PS</th><th>3D ScatterNet</th></tr></thead><tbody><tr><td>log<sub>10</sub>(T<sub>vir</sub>/K)</td><td>5.477</td><td>5.485<sup>+0.037</sup><sub>-0.036</sub></td><td>5.480<sup>+0.015</sup><sub>-0.016</sub></td><td>5.476<sup>+0.010</sup><sub>-0.010</sub></td><td>5.464<sup>+0.063</sup><sub>-0.088</sub></td><td>5.446<sup>+0.060</sup><sub>-0.068</sub></td><td>5.341<sup>+0.128</sup><sub>-0.139</sub></td><td>5.349<sup>+0.132</sup><sub>-0.137</sub></td></tr><tr><td>log<sub>10</sub>(ζ)</td><td>2.301</td><td>2.307<sup>+0.036</sup><sub>-0.033</sub></td><td>2.306<sup>+0.023</sup><sub>-0.023</sub></td><td>2.308<sup>+0.021</sup><sub>-0.019</sub></td><td>2.279<sup>+0.078</sup><sub>-0.110</sub></td><td>2.258<sup>+0.077</sup><sub>-0.085</sub></td><td>2.117<sup>+0.149</sup><sub>-0.155</sub></td><td>2.122<sup>+0.164</sup><sub>-0.154</sub></td></tr></tbody></table>
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** The posteriors estimated from 3D cosmic 21 cm light-cone images (i.e., without thermal noise or residual foreground) by three different approaches of data compression—power spectrum using 21cmDELFI-PS (green), 3D CNN using DELFI-3D CNN (blue), and solid harmonic WST using 3D ScatterNet (red)—for two mock observations: the Faint Galaxies Model (left) and the Bright Galaxies Model (right). We show the $1\sigma$ (dark) and $2\sigma$ (light) credible regions. The dashed lines indicate the true parameter values.</figcaption></figure>

<table><caption>**Table 4**
Recovery Performance by the 21cmDELFI-PS and 3D ScatterNet</caption><thead><tr><th rowspan="2"></th><th colspan="2">Pure Signal</th><th colspan="2">SKA Noise</th><th colspan="2">SKA Noise + Residual Foreground</th></tr><tr><th>21cmDELFI-PS</th><th>3D ScatterNet</th><th>21cmDELFI-PS</th><th>3D ScatterNet</th><th>21cmDELFI-PS</th><th>3D ScatterNet</th></tr></thead><tbody><tr><td>R<sup>2</sup> <sup>a</sup></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>log<sub>10</sub>(T<sub>vir</sub>)</td><td>0.9989</td><td>0.9997</td><td>0.9336</td><td>0.9647</td><td>0.7981</td><td>0.8348</td></tr><tr><td>log<sub>10</sub>(ζ)</td><td>0.9978</td><td>0.9990</td><td>0.9353</td><td>0.9604</td><td>0.8028</td><td>0.8254</td></tr><tr><td>ε <sup>b</sup></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>T<sub>vir</sub></td><td>(-0.03, 0.03)</td><td>(-0.02, 0.02)</td><td>(-0.19, 0.30)</td><td>(-0.17, 0.16)</td><td>(-0.28, 0.47)</td><td>(-0.28, 0.41)</td></tr><tr><td>ζ</td><td>(-0.04, 0.03)</td><td>(-0.02, 0.02)</td><td>(-0.15, 0.24)</td><td>(-0.12, 0.14)</td><td>(-0.26, 0.36)</td><td>(-0.22, 0.32)</td></tr></tbody></table>

**Notes.**

<sup>a</sup> The coefficient of determination R<sup>2</sup> is computed for the medians (in the logarithmic scale) of the inferred posteriors from 300 testing samples.

<sup>b</sup> The fractional error ϵ refers to the relative error of the deduced parameters in the linear scale. Here, we present the 68% credible interval of the probability density distribution of ϵ.

<table><caption>**Table 5**
The p-values for the Null Hypotheses That These Statistics Are of a Uniform Distribution with 300 Testing Samples</caption><thead><tr><th rowspan="3">Statistics</th><th colspan="4">Pure Signal</th><th colspan="4">SKA Noise</th><th colspan="4">SKA Noise + Residual Foreground</th></tr><tr><th colspan="2">21cmDELFI-PS</th><th colspan="2">3D ScatterNet</th><th colspan="2">21cmDELFI-PS</th><th colspan="2">3D ScatterNet</th><th colspan="2">21cmDELFI-PS</th><th colspan="2">3D ScatterNet</th></tr><tr><th>KS</th><th>CoM</th><th>KS</th><th>CoM</th><th>KS</th><th>CoM</th><th>KS</th><th>CoM</th><th>KS</th><th>CoM</th><th>KS</th><th>CoM</th></tr></thead><tbody><tr><td>PIT (T<sub>vir</sub>)</td><td>0.37</td><td>0.13</td><td>0.80</td><td>0.65</td><td>0.01</td><td>0.02</td><td>0.20</td><td>0.52</td><td>0.19</td><td>0.17</td><td>0.42</td><td>0.37</td></tr><tr><td>PIT (ζ)</td><td>0.14</td><td>0.05</td><td>0.92</td><td>0.75</td><td>0.05</td><td>0.05</td><td>0.78</td><td>0.65</td><td>0.01</td><td>0.02</td><td>0.18</td><td>0.24</td></tr><tr><td>copPIT</td><td>0.54</td><td>0.39</td><td>0.67</td><td>0.92</td><td>0.01</td><td>0.03</td><td>0.64</td><td>0.85</td><td>0.53</td><td>0.31</td><td>0.46</td><td>0.46</td></tr><tr><td>HPD</td><td>0.04</td><td>0.01</td><td>0.63</td><td>0.48</td><td>0.01</td><td>0.01</td><td>0.41</td><td>0.47</td><td>0.74</td><td>0.67</td><td>0.05</td><td>0.03</td></tr></tbody></table>

3D ScatterNet and 21cmDELFI-PS give a very high R<sup>2</sup> score, but the 3D ScatterNet slightly outperforms the 21cmDELFI-PS.

Table 4 also presents the 68% credible interval of the probability density distribution of ϵ, the fractional errors of the deduced parameter x in the linear scale (e.g., T<sub>vir</sub> and ζ in this paper). The fractional error refers to that of the predicted median with respect to the true value, i.e., $ϵ = (x_{pred} - x_{true})/x_{true}$. For

the case of a pure signal, the typical fractional error (represented by the length of the interval) of 3D ScatterNet is about 1.50 times smaller than that of 21cmDELFI-PS for T<sub>vir</sub>, and 1.75 times smaller for ζ.

Lastly, we perform the posterior validation. Table 5 shows the p-value for some hypothesis tests using 300 testing samples. Note that all p-values are larger than 0.01 and most are larger than 0.5, which implies that our results are at least
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 6.** Same as Figure 5, but the estimations are made from mock observations of the 3D 21 cm light-cone images from SKA, which includes the total noise (with the contributions from thermal noise and sample variance errors) yet without foreground contamination. Here, we only compare two different approaches to data compression—power spectrum using 21cmDELFI-PS (green) and solid harmonic WST using 3D ScatterNet (red).</figcaption></figure>

reliable with a significance of 0.01. Also, the *p*-values of 3D ScatterNet are generically larger than those of 21cmDELFI-PS.

### 4.2.2. SKA Noise and Residual Foreground

We now consider the case of “SKA noise,” which refers to the mock SKA observations of the 3D 21 cm images with total noise (with contributions from thermal noise and sample variance errors) yet without foreground contamination. Since the images after smoothing with 1 km baseline lose small-scale information, we discard the components of scattering coefficients with *j* = 0 and hence the final scattering coefficients have the dimension of 48. For the power spectrum, we also discard the largest *k*-modes in a single box and generate the final power spectrum vector with a dimension of 120.

For the case of SKA noise, Figure 6 shows the credible regions for the Faint Galaxies Model and the Bright Galaxies Model, with quantitative comparisons in terms of medians and $1\sigma$ errors presented in Tables 2 and 3. The inference results of the SKA noise case have larger uncertainties than those of pure signal, which is reasonable, given the noise. The statistical errors of 3D ScatterNet are smaller than those of 21cmDELFI-PS. This improvement is also verified by Table 4: the $R^2$ values of 3D ScatterNet are larger than those of 21cmDELFI-PS, and the typical fractional error of the former is about 1.49 (1.50) times smaller for $T_{\text{vir}}$ (ζ) than that of the latter.

Lastly, we consider the case of “SKA noise + residual foreground,” which refers to the mock SKA observations of the 3D 21 cm images with total noise and residual foreground after the foreground is removed with the SVD technique. For the same reason of smoothing as in the SKA noise case, we reduce the dimension of the scattering coefficients to be 48 by discarding the components with *j* = 0. For the power spectrum, in order to train reliable NDEs, we have to further discard more large *k*-modes with the upper limit $k_{\text{max}}$ = 0.45 Mpc$^{-1}$ and the final vector of the power spectrum has the dimension of 70. In

this case, the information for parameter inference is only from large-scale modes because the images are smoothed with a rather coarse resolution.

For the case of SKA noise + residual foreground, Figure 7 shows the inference results for two mock observations, with the medians and $1\sigma$ errors listed in Tables 2 and 3. As expected, given more uncertainties due to residual foreground, the posteriors have larger errors in this case than those of SKA noise. The statistical errors of 3D ScatterNet are on average smaller than those of 21cmDELFI-PS. This improvement is also verified by Table 4: the $R^2$ values of 3D ScatterNet are higher than those of 21cmDELFI-PS, and the typical fractional error of the former is about 1.09 (1.15) times smaller for $T_{\text{vir}}$ (ζ) than that of the latter.

Regarding the posterior validation, all *p*-values are larger than 0.01 both for the SKA noise case and the case of SKA noise + residual foreground, which implies that our results are at least reliable with a significance of 0.01. Also, the *p*-values of 3D ScatterNet are generically larger than those of 21cmDELFI-PS.

In sum, for all cases of different assumptions in noise and residual foreground, 3D ScatterNet outperforms 21cmDELFI-PS. Our results demonstrate that the solid harmonic WST can extract the information from the 3D 21 cm light-cone images in a more effective way than the power spectrum statistic itself. These results are consistent with the previous findings that the *n*-th-order scattering transform captures the information up to the $2^n$-point function (Cheng & Ménard 2021; Mallat 2012). In addition, Sui et al. (2023) introduced the concept of mutual information to evaluate the comparisons between different data summaries by giving them scores and found that the solid harmonic WST not only outperforms the power spectrum but also surpasses the combined information, including both power spectrum and bispectrum, in the context of extracting the information from the 21 cm data cube.
<figure><img src="image_12.png" /><figcaption>**Figure 7.** Same as Figure 6, but with mock observations of the 3D 21 cm light-cone images from SKA, which includes the total noise and residual foreground after the foreground is removed using SVD.</figcaption></figure>

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

## 4.3. Discussion

In this subsection, we attempt to provide some insights regarding why the 3D CNN (trained with the particular set of experimental choices in Paper I) performs less effectively than the solid harmonic WST (and power spectrum) in compressing the 3D 21 cm light-cone images, and the limitations of the 3D CNN configurations in Paper I.

The solid harmonic WST is analogous to a 3D CNN, but unlike the latter, it essentially gives a fixed kernel without the training process. Furthermore, the resulting scattering coefficients are naturally invariant to translations and rotations and are Lipschitz continuous to deformations (Mallat 2012a). These properties are particularly useful for 3D 21 cm light-cone images because the data summaries have the translational invariance in the angular direction and the rotational invariance along the LOS,<fnref n="17" /> and can be stable under slight deformations due to variations in the initial density fields. These intrinsic invariance properties can enhance conditional density learning of the solid harmonic WST.

In principle, thanks to the similarity between solid harmonic WST and 3D CNN, a deep and/or wide 3D CNN should perform at least as well as the WST in extracting information, as long as three conditions are all met—training data is adequate, network is sufficiently expressive, and training process is successful. Our results imply that the 3D CNN configurations in Paper I might have limitations in these three aspects.

The first limitation is insufficient training data. A data set that samples the complete parameter-data distribution is required to train a robust neural network. For training a 2D CNN (Gillet et al. 2019), a data set comprising 9000 training samples might suffice, but the same number of training samples is likely not adequate for 3D CNN.<fnref n="18" /> While we find that the

results for the Faint Galaxies Model and the Bright Galaxies Model were not improved if the training sample size was doubled, a thorough scaling test on the number of training samples is necessary to make a robust diagnosis.

The second limitation is the insufficient complexity of neural network architectures. While the solid harmonic WST exploits its intrinsic invariance properties, the standard CNNs cannot because they are not invariant to rotations due to inherent architecture. Thus the architecture in CNNs needs to be fine-tuned during the training. Alternatively, a likely more effective approach is to integrate additional inductive biases into the architecture (Kauderer-Abrams 2017; Weiler et al. 2018a, 2018b; Semih Kayhan & van Gemert 2020). In these variants of CNNs, the invariance properties can be enhanced, which might optimize the network architecture. Also, in Neutsch et al. (2022), an elongated kernel is designed to detect the patterns evolving in the redshift direction, which reduces the number of trainable parameters.<fnref n="19" />

The third limitation is the suboptimal network training process. Network optimization involves refined practices for neural network initialization and learning rate determination (DeZoort & Hanin 2023). Some generic optimization strategies (Akiba et al. 2019; Neutsch et al. 2022) may be useful as well.

Note that the 3D CNN trained in Paper I performs even less effectively than the power spectrum analysis using 21cmDELFI-PS. This might be explained by the limitations of 3D CNN discussed above. While the 3D CNN configurations in Paper I might be improved in these three aspects, we leave a thorough exploration along these lines to future work.

The comparison involving 3D CNNs in this work has been limited to the case of a pure signal. This is because the comparison when instrumental effects are applied would require additional fine-tuning efforts for the 3D CNN to fully assess its potential, which is not the focus of this work. Appendix E provides some insight into such a comparison by

<footnote n="17">¹⁷ In practice, we obtain approximate invariance that is limited by the boundaries and the discrete sampling of images.</footnote>

<footnote n="18">¹⁸ We note that Neutsch et al. (2022) trained a 3D CNN with a smaller number of data than Paper I, but it performed comparably well. This might benefit from the specific architecture adopted therein.</footnote>

<footnote n="19">¹⁹ Similarly, in this work, we also use wavelet filters with elongated shapes for solid harmonic WST. Prelogović et al. (2022) included some discussions on the neural network architectures used for interpreting the 21 cm light cones.</footnote>
including the SKA noise in the application of DELFI-3D CNN, nevertheless without additional fine-tuning, i.e., using the same hyperparameters for training the 3D CNN as detailed in Paper I.

## 5. Summary

In this paper, we introduce the solid harmonic WST for compressing the 3D image data and generating meaningful low-dimensional summaries. We apply this technique to the data compression of 3D tomographic 21 cm light-cone images and use the resulting scattering coefficients as the input summaries of DELFI. With DELFI, we perform the Bayesian inference of the reionization parameters where the likelihood is implicitly defined by the forward simulations. This new technique, dubbed 3D ScatterNet (i.e., solid harmonic WST with DELFI), recovers accurate posterior distributions for the reionization parameters.

We compare the inference results with two different approaches of data compression, the solid harmonic WST and the 3D CNN, for two representative mock observations in the pure signal case, and demonstrate that 3D ScatterNet outperforms DELFI-3D CNN significantly. Our results imply that the solid harmonic WST extracts the information from the 3D 21 cm light-cone images in a more informative manner than a 3D CNN given reasonable fine-tuning. We also offer insights into enhancing the 3D CNN regarding the data sets, network architecture, and training process. Moreover, the solid harmonic WST has fixed kernels in wavelets, which means that it does not need to be trained in order to output the data summaries. This highlights its robustness and efficiency, which is another advantage of the solid harmonic WST over the 3D CNN.

We then make another comparison of the inference results between data compression using the solid harmonic WST (3D ScatterNet) and using the power spectrum statistic (21cmDELFI-PS), both under the same inference strategy using DELFI. The comparisons were made for three cases with different assumptions on the noise and residual foreground— pure signal, SKA noise, and SKA noise + residual foreground. We find that the 3D ScatterNet outperforms 21cmDELFI-PS in all cases. This implies that the summaries compressed with the solid harmonic WST contain more (i.e., non-Gaussian) information from the 3D 21 cm light-cone images than the power spectrum analysis.

Our results demonstrate that combining a WST (3D solid harmonic WST in this paper) with the SBI will be a promising tool for the scientific interpretation of future 21 cm light-cone image observation data. Based on the findings of this paper, there is room for improvement with regard to the design of summary statistics from solid harmonic WST or new WST. For example, instead of treating the light cone as a whole to build a statistic, a variant approach might be to treat discrete boxes separately and concatenate the scattering coefficients of each box together as the data summaries. In this way, the correlations between different stages of reionization can be exploited by the CMAFs. Nevertheless, this will increase the dimension of data summaries, and therefore, incur higher computational costs, so the trade-off between accuracy and efficiency in the solid harmonic WST is yet to be further explored. Also, further cross correlation between scales resulting from rescaling the phase information is likely to be exploited in a 3D extension of the wavelet phase harmonics

(Mallat et al. 2019; Allys et al. 2020) for the 2D case. On the other hand, when scaling the similar analysis to a high-dimensional parameter and/or data space, the diffusion model (Luo 2022; Zhao et al. 2023a; Legin et al. 2024) is a promising alternative to the CMAFs, or other normalizing flows, to perform the SBI because of its simplified training objective. We leave the exploration of these directions to future works.

*Note*.—In the final stage of preparing this manuscript, the papers by Eickenberg et al. (2022) and Greig et al. (2022) were posted on arXiv. These papers applied the scattering transform (or wavelet moments, varied band-limited first-order scattering coefficients defined in the Fourier space, in Eickenberg et al. 2022) for parameter inference, but both used the Fisher matrix formalism. In comparison, we apply the scattering transform for parameter inference using SBI in this paper. In addition, Greig et al. (2022) focus on the reionization parameter estimation using the 2D 21 cm images; Eickenberg et al. (2022) focus on the cosmological parameter estimation using the 3D cosmological density fields. In comparison, our paper has a distinct focus, i.e., on the reionization parameter estimation using the 3D 21 cm light-cone images.

## Acknowledgments

This work is supported by the National SKA Program of China (grant No. 2020SKA0110401), NSFC (grant No. 11821303), and the National Key R&D Program of China (grant No. 2018YFA0404502). B.D.W. acknowledges support from the Simons Foundation. We thank Sihao Cheng, Paulo Montero-Camacho, and Bohua Li for their useful discussions and help. We acknowledge the Tsinghua Astrophysics High-Performance Computing platform at Tsinghua University for providing computational and data storage resources that have contributed to the research results reported within this paper.

*Software:* 21CMMC (Greig & Mesinger 2015, 2017, 2018), 21cmFAST (Mesinger & Furlanetto 2007; Mesinger et al. 2011), pydelfi (Alsing et al. 2019), TensorFlow (Abadi et al. 2016), GetDist (Lewis 2019), NumPy (Harris et al. 2020), Matplotlib (Hunter 2007), SciPy (Virtanen et al. 2020), scikit-learn (Pedregosa et al. 2011), Python2 (Van Rossum & Drake 1995), Python3 (Van Rossum & Drake 2009), Kymatio (Andreux et al. 2018), Tools21cm (Giri et al. 2020), pygsm (de Oliveira-Costa et al. 2008; Zheng et al. 2017), Mayavi (Ramachandran & Varoquaux 2011), galpro (Mucesh et al. 2021), seaborn (Waskom 2021), Astropy (Astropy Collaboration et al. 2013, 2018), healpy (Zonca et al. 2019), HEALPix (Górski et al. 2005).

## Appendix A
### The Effect of Angular Frequency $\ell$

In the main text of this paper, the scattering coefficients are averaged over $0 \le \ell \le \ell_{\text{max}}$ with $\ell_{\text{max}} = 6$. In this section, we explore the effect of the angular frequency $\ell$ on the parameter inference.

Figure 8 shows the scattering coefficients evaluated at given values of $\ell$. While the first-order scattering coefficients show tilting when varying $\ell$, the second-order scattering coefficients decrease in overall amplitude as $\ell$ increases. Also, for $\ell=0$ (i.e., Gaussian wavelets), the scattering coefficients are flat at the scale of $j'$.

We construct two sets of experiments of posterior inference over 300 testing samples and show the coefficient of
<figure><img src="image_14.png" /><figcaption>**Figure 8.** The first-order ($S_1$) and second-order ($S_2$) scattering coefficients of the simulated cosmological 21 cm signals, evaluated at given values of $\ell$ and at $q = 1$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 9.** The coefficient of determination $R^2$ as a function of the angular frequency $\ell$. The phrase $\ell$ single refers to the scenario where the scattering coefficients are evaluated at a given single value of $\ell$, while the phrase $\ell$ average refers to the scenario where the scattering coefficients are averaged over $0 \le \ell' \le \ell$ for a given value of $\ell$. The parameter after this phrase (i.e., $T_{vir}$ or $\zeta$) indicates the parameter of the $R^2$ calculation.</figcaption></figure>

determination $R^2$ of the predicted medians in Figure 9. In the first set of experiments (“$\ell$ single”), the scattering coefficients are evaluated at each single value of $\ell$. We find that using the $\ell = 0$ information has the smallest values of $R^2$. When using the information at larger $\ell$, the $R^2$ value increases and then slightly decreases. This implies that the solid harmonic wavelet has better performance than the Gaussian wavelet (i.e., the case with $\ell = 0$). In the second set of experiments (“$\ell$ average”), the scattering coefficients are averaged over $0 \le \ell' \le \ell$ for a given value of $\ell$. Comparing these two scenarios, we find that the combined information ($\ell$ average) gets a higher value of $R^2$ generically than the single information ($\ell$ single). When increasing the value of $\ell$, in the case of $\ell$ average, the $R^2$ value increases and then slightly decreases, with the peak at $\ell = 6$. For this reason, in the main text of this paper, we choose the $\ell$ average scenario with $\ell_{max} = 6$.

## Appendix B

### The Effect of Smoothing Scales

In the main text of this paper, the images are smoothed with the scale corresponding to the size of the 1 km baseline. In this section, we test the effect of smoothing scales on parameter

<table><caption>**Table 6**
Recovery Performance by 3D ScatterNet for Smoothing with the Baseline of 1 km and 2 km, Respectively</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">SKA Noise</th><th colspan="2">SKA Noise + Residual Foreground</th></tr><tr><th>1 km</th><th>2 km</th><th>1 km</th><th>2 km</th></tr></thead><tbody><tr><td rowspan="2">R<sup>2</sup></td><td>log<sub>10</sub>(T<sub>vir</sub>)</td><td>0.9647</td><td>0.8808</td><td>0.8348</td><td>0.7853</td></tr><tr><td>log<sub>10</sub>(ζ)</td><td>0.9604</td><td>0.9153</td><td>0.8254</td><td>0.8212</td></tr><tr><td rowspan="2">ε</td><td>T<sub>vir</sub></td><td>(-0.17, 0.16)</td><td>(-0.18, 0.27)</td><td>(-0.28, 0.41)</td><td>(-0.27, 0.44)</td></tr><tr><td>ζ</td><td>(-0.12, 0.14)</td><td>(-0.16, 0.20)</td><td>(-0.22, 0.32)</td><td>(-0.25, 0.29)</td></tr></tbody></table>

inference. We compute the values of $R^2$ and typical fractional error $\epsilon$ for the predicted medians of 300 testing samples. Table 6 shows that smoothing with the 1 km baseline has a larger $R^2$ value and smaller fractional error than smoothing with the 2 km baseline. This is because the smoothed images in the former have a higher signal-to-noise ratio, although in coarser resolution. This result is the reason why we chose the size of the 1 km baseline for smoothing in this paper.

## Appendix C

### The Dependence of Light-cone Volumes

In the main text of this paper, the full light-cone images, which concatenate ten light-cone boxes, are exploited for parameter inference. In this section, we investigate whether the information from one of these boxes dominates over the others. (If this were the case, then observations could be optimized targeting at the “sweet-spot” redshift.)

For this purpose, we perform the parameter inference from discrete boxes, both using the power spectrum and using scattering coefficients as summaries. Figure 10 shows the posteriors for two mock observations and Table 7 shows the results of recovery performance for the predicted medians of 300 testing samples. For 3D ScatterNet, the box at the low-redshift results in better inference performance in terms of larger $R^2$, smaller typical fractional error, and smaller credible region for mock observations than that at the high redshift. However, no single box appears to dominate the information, because the full light-cone images always yield the optimal result that is significantly better than any of the single boxes. These results hold similarly for 21cmDELFI-PS. We
<figure><img src="image_16.png" /></figure>

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

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

<figure><img src="image_19.png" /><figcaption>**Figure 10.** The posteriors estimated from 3D cosmic 21 cm light-cone images (i.e., pure signal) for the scenario where only the information within the first (red), third (green), fifth (blue) box, and the full (black) light-cone images is exploited for parameter inference, respectively. Shown are the results using 21cmDELFI-PS (top) and using 3D ScatterNet (bottom) for the Faint Galaxies Model (left) and the Bright Galaxies Model (right), respectively.</figcaption></figure>

<table><caption>**Table 7**
Recovery Performance for Information in Discrete Boxes</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="4">21cmDELFI-PS</th><th colspan="4">3D ScatterNet</th></tr><tr><th>First</th><th>Third</th><th>Fifth</th><th>Full</th><th>First</th><th>Third</th><th>Fifth</th><th>Full</th></tr></thead><tbody><tr><td rowspan="2">R<sup>2</sup></td><td>log<sub>10</sub>(T<sub>vir</sub>)</td><td>*<sup>a</sup></td><td>0.9948</td><td>0.9863</td><td>0.9989</td><td>0.9991</td><td>0.9987</td><td>0.9978</td><td>0.9997</td></tr><tr><td>log<sub>10</sub>(ζ)</td><td>*</td><td>0.9909</td><td>0.9754</td><td>0.9978</td><td>0.9960</td><td>0.9953</td><td>0.9924</td><td>0.9990</td></tr><tr><td rowspan="2">ε</td><td>T<sub>vir</sub></td><td>*</td><td>(-0.06, 0.08)</td><td>(-0.07, 0.08)</td><td>(-0.03, 0.03)</td><td>(-0.03, 0.03)</td><td>(-0.03, 0.04)</td><td>(-0.05, 0.04)</td><td>(-0.02, 0.02)</td></tr><tr><td>ζ</td><td>*</td><td>(-0.05, 0.08)</td><td>(-0.08, 0.07)</td><td>(-0.04, 0.03)</td><td>(-0.04, 0.05)</td><td>(-0.04, 0.05)</td><td>(-0.06, 0.04)</td><td>(-0.02, 0.02)</td></tr></tbody></table>

**Notes.** The parameter inference is made for the case of pure signal. Here, “First,” “Third” and “Fifth” refer to the scenario where only the information within the first, third, and fifth box (counting from the low redshift, corresponding to the central redshift of 7.67, 8.36, and 9.11) is exploited for parameter inference, respectively. “Full” refers to the case of exploiting the full light-cone images (i.e., concatenating all ten boxes), which is the case in the main text of this paper.

<footnote><sup>a</sup> We do not show these results for the first box because the hypothesis test with the HPD values failed, likely due to the lack of training samples.</footnote>
conclude that there is no optimal redshift. As such, the full light-cone images, which gather the complete information, provide the parameter inference with the best performance.

## Appendix D
### The NDE Setting and Sample Size

In this section, we give the details of the networks as follows. We choose the CMAFs as the NDEs throughout this paper. In all architectures, we set two neural layers of a single transform, represented by the masked autoencoders for density estimation (Germain et al. 2015), with 50 neurons per layer. We also use the ensembles of NDEs. The final posterior is the stacked one from individual posteriors with weights according to their training errors. We fine-tune the training sample size and the CMAFs' architecture for every method and experiment based on the outcomes of posterior validation (calibration). The guiding principle behind this fine-tuning is that data (summaries) with greater uncertainties might require more intricate CMAFs, such as increased transformations within a single CMAF and the increased number of CMAFs within the ensemble. Meanwhile, data (summaries) of higher dimensions might necessitate a more substantial training sample due to the vast feature space.

The procedure for scattering coefficients is as follows.

(1) The scattering coefficients for the pure signal case with the $k_{\perp} = 0$ mode removed (see Section 4.1). We use 18,000 light-cone images and an ensemble of eight NDEs (5, 6, 7, 8)\*2, which means that we have two NDE blocks, each having the NDE with the number of transformations 5, 6, 7, and 8, respectively. Hereafter we use the same terminology. Notice that the final performance can be further enhanced with double the size of the light-cone images. To fairly compare with 21cmDELFI-PS, we use 18,000 images for reports in the paper.

(2) The coefficients from pure discrete boxes (see Appendix C). We use 27,000 samples. The ensembles for the first, third, and fifth boxes: (5, 6, 7, 8)\*3, (5, 6, 7, 8)\*3, (5, 6, 7, 8)\*2, respectively.

(3) The coefficients with the information of averaged $\ell$ from pure boxes (see Appendix A). We use 36,000 samples. The ensembles from $\ell = 1$ to $\ell = 8$: (5)\*4, (6, 7, 8, 9)\*3, (6, 7, 8, 9)\*2, (5, 6, 7, 8)\*2, (6, 7, 8, 9)\*3, (5, 6, 7, 8)\*2, (5)\*4, (5, 6, 7, 8)\*2, respectively.

(4) The coefficients with the information of single $\ell$ from pure boxes (see Appendix A). We use 36,000 samples.

The ensembles from $\ell = 0$ to $\ell = 8$: (5)\*4, (5)\*4, (5, 6, 7, 8)\*2, (5)\*4, (5, 6, 7, 8)\*2, (5, 6, 7, 8)\*2, (5, 6, 7, 8)\*2, and (5, 6, 7, 8)\*2, respectively.

(5) The coefficients from the signal with SKA noise. We use 27,000 samples (see Section 4.2.2). The ensembles for the signal with 1 km smoothing: (5, 6, 7, 8)\*2, and 2 km smoothing: (5, 6, 7, 8)\*3, respectively.

(6) The coefficients from the signal with SKA noise and residual foreground (see Section 4.2.2). We use 36,000 samples. The ensembles for the signal with 1 km smoothing: (20)\*2, and 2 km smoothing: (30)\*2, respectively.

The procedure for the power spectrum is as follows.

(1) The power spectrum from a pure signal with the $k_{\perp} = 0$ mode removed (see Section 4.1). We use 18,000 samples and the ensembles (5)\*4.

(2) The power spectrum from pure discrete boxes (see Appendix C). For the first box, we use 32,596 samples and the ensembles (6, 7, 8, 9)\*3; for the third box, we use 34,044 samples and the ensembles (5, 6, 7, 8)\*3; and for the fifth box, we use 27,000 samples and the ensembles (5, 6, 7, 8)\*2.

(3) The power spectrum from the signal with SKA noise (see Section 4.2.2). We use 27,000 samples and the ensembles (5)\*4.

(4) The power spectrum from the signal with SKA noise and residual foreground (see Section 4.2.2). We use 36,000 samples and the ensembles (20)\*3.

## Appendix E
### Performance of DELFI-3D CNN with SKA Noise

In this section, we extend the comparison of 3D CNN using DELFI-3D CNN, power spectrum using 21cmDELFI-PS and solid harmonic WST using 3D ScatterNet to the mock observations of the 3D 21 cm light-cone images from SKA, which includes the total noise yet without foreground contamination. Ideally, this would require additional fine-tuning efforts for the 3D CNN to fully assess its potential. Here, we present the results of DELFI-3D CNN without additional fine-tuning, i.e., using the same hyperparameters for training the 3D CNN and NDEs as detailed in Paper I. Figure 11 shows that in this realistic scenario, a trained 3D CNN still performs suboptimal compared to 21cmDELFI-PS and 3D ScatterNet.
<figure><img src="image_20.png" /><figcaption>Figure 11. Same as Figure 6 but including the results from DELFI-3D CNN.</figcaption></figure>

## ORCID iDs

Xiaosheng Zhao https://orcid.org/0000-0002-8328-1447

Yi Mao https://orcid.org/0000-0002-1301-3893

Shifan Zuo https://orcid.org/0000-0003-3858-6361

Benjamin D. Wandelt https://orcid.org/0000-0002-5854-8269

## References

Abadi, M., Barham, P., Chen, J., et al. 2016, in Proc. of the 12th USENIX Conference on Operating Systems Design and Implementation, OSDI’16 (Berkeley, CA: USENIX Association), 265

Abdurashidova, Z., Aguirre, J. E., Alexander, P., et al. 2022, ApJ, 925, 221

Akiba, T., Sano, S., Yanase, T., Ohta, T., & Koyama, M. 2019, arXiv:1907.10902

Allys, E., Levrier, F., Zhang, S., et al. 2019, A&A, 629, A115

Allys, E., Marchand, T., Cardoso, J. F., et al. 2020, PhRvD, 102, 103506

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

Alsing, J., Wandelt, B., & Feeney, S. 2018, MNRAS, 477, 2874

Anderson, T. W. 1962, Ann. Math. Stat., 33, 1148

Andreux, M., Angles, T., Exarchakis, G., et al. 2018, arXiv:1812.11214

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Bianco, M., Giri, S. K., Iliev, I. T., & Mellema, G. 2021, MNRAS, 505, 3982

Cheng, S., & Ménard, B. 2021, arXiv:2112.01288

Cheng, S., Ting, Y.-S., Ménard, B., & Bruna, J. 2020, MNRAS, 499, 5902

Chung, D. T. 2022, MNRAS, 517, 1625

Cranmer, K., Brehmer, J., & Louppe, G. 2020, PNAS, 117, 30055

de Oliveira-Costa, A., Tegmark, M., Gaensler, B. M., et al. 2008, MNRAS, 388, 247

DeBoer, D. R., Parsons, A. R., Aguirre, J. E., et al. 2017, PASP, 129, 045001

DeZoort, G., & Hanin, B. 2023, arXiv:2306.11668

Eickenberg, M., Allys, E., Moradinezhad Dizgah, A., et al. 2022, arXiv:2204.07646

Eickenberg, M., Exarchakis, G., Hirn, M., & Mallat, S. 2017, in Proc. of the 31st Int. Conf. on Neural Information Processing Systems, NIPS’17 (Red Hook, NY: Curran Associates Inc.), 6543

Eickenberg, M., Exarchakis, G., Hirn, M., Mallat, S., & Thiry, L. 2018, JChPh, 148, 241732

Furlanetto, S. R., Oh, S. P., & Briggs, F. H. 2006, PhR, 433, 181

Gauthier, S., Thérien, B., Alsène-Racicot, L., et al. 2021, arXiv:2107.09539

Germain, M., Gregor, K., Murray, I., & Larochelle, H. 2015, ICML, 37, 881, https://proceedings.mlr.press/v37/germain15.html

Gillet, N., Mesinger, A., Greig, B., Liu, A., & Ucci, G. 2019, MNRAS, 484, 282

Giri, S. K., & Mellema, G. 2021, MNRAS, 505, 1863

Giri, S. K., Mellema, G., & Ghara, R. 2018, MNRAS, 479, 5596

Giri, S. K., Mellema, G., & Jensen, H. 2020, JOSS, 5, 2363

Gneiting, T., Balabdaoui, F., & Raftery, A. E. 2007, J. R. Stat, 69, 243

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Greig, B., & Mesinger, A. 2015, MNRAS, 449, 4246

Greig, B., & Mesinger, A. 2017, MNRAS, 472, 2651

Greig, B., & Mesinger, A. 2018, MNRAS, 477, 3217

Greig, B., Mesinger, A., Haiman, Z., & Simcoe, R. A. 2017, MNRAS, 466, 4239

Greig, B., Ting, Y.-S., & Kaurov, A. A. 2022, MNRAS, 513, 1719

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Harrison, D., Sutton, D., Carvalho, P., & Hobson, M. 2015, MNRAS, 451, 2610

Hunter, J. D. 2007, CSE, 9, 90

Intema, H. T., Jagannathan, P., Mooley, K. P., & Frail, D. A. 2017, A&A, 598, A78

Ivezić, Ž., Connolly, VanderPlas, A. J., & Gray, A. J. T. 2014, Statistics, Data Mining, and Machine Learning in Astronomy: A Practical Python Guide for the Analysis of Survey Data (Princeton, NJ: Princeton Univ. Press)

Jacobs, D. C., Pober, J. C., Parsons, A. R., et al. 2015, ApJ, 801, 51

Jeffrey, N., Boulanger, F., Wandelt, B. D., et al. 2022, MNRAS, 510, L1

Kauderer-Abrams, E. 2017, arXiv:1801.01450

Kolmogorov, A. 1992, in On the Empirical Determination of a Distribution Function, ed. S. Kotz & N. L. Johnson (New York: Springer), 106

Legin, R., Ho, M., Lemos, P., et al. 2024, MNRAS, 527, L173

Lewis, A. 2019, arXiv:1910.13970

Luo, C. 2022, arXiv:2208.11970

Mallat, S. 2012, Commun. Pure Appl. Math., 65, 1331

Mallat, S., Zhang, S., & Rochette, G. 2019, Inf. Inference, 9, 721

Mao, Y., Shapiro, P. R., Mellema, G., et al. 2012, MNRAS, 422, 926

Masui, K. W., Switzer, E. R., Banavar, N., et al. 2013, ApJL, 763, L20

Mellema, G., Koopmans, L. V. E., Abdalla, F. A., et al. 2013, ExA, 36, 235

Mertens, F. G., Mevius, M., Koopmans, L. V. E., et al. 2020, MNRAS, 493, 1662

Mesinger, A., & Furlanetto, S. 2007, ApJ, 669, 663

Mesinger, A., Furlanetto, S., & Cen, R. 2011, MNRAS, 411, 955

Mucesh, S., Hartley, W. G., Palmese, A., et al. 2021, MNRAS, 502, 2770

Neutsch, S., Heneka, C., & Brüggen, M. 2022, MNRAS, 511, 3446

Paciga, G., Albert, J. G., Bandura, K., et al. 2013, MNRAS, 433, 639

Papamakarios, G. 2019, arXiv:1910.13233

Papamakarios, G., Nalisnick, E., Rezende, D. J., Mohamed, S., & Lakshminarayanan, B. 2021, JMLR, 22, 1, https://www.jmlr.org/papers/v22/19-1028.html
Papamakarios, G., Pavlakou, T., & Murray, I. 2017, in Proc. of the 31st Int. Conf. on Neural Information Processing Systems, NIPS’17 (Red Hook, NY: Curran Associates Inc.), 2335

Parsons, A. R., Backer, D. C., Foster, G. S., et al. 2010, ApJ, 139, 1468

Parsons, A. R., Liu, A., Aguirre, J. E., et al. 2014, ApJ, 788, 106

Pedersen, C., Ho, S., & Eickenberg, M. 2022, in Machine Learning for Astrophysics (Baltimore, MD: ICML), 40

Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, JMLR, 12, 2825

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016, A&A, 594, A13

Prelogović, D., & Mesinger, A. 2023, MNRAS, 524, 4239

Prelogović, D., Mesinger, A., Murray, S., Fiameni, G., & Gillet, N. 2022, MNRAS, 509, 3852

Ramachandran, P., & Varoquaux, G. 2011, CSE, 13, 40

Saydjari, A. K., Portillo, S. K. N., Slepian, Z., et al. 2021, ApJ, 910, 122

Semih Kayhan, O., & van Gemert, J. C. 2020, in 2020 IEEE/CVF Conf. on Computer Vision and Pattern Recognition (CVPR) (Piscataway, NJ: IEEE), 14262

Stewart, G. W. 1993, SIAM Rev., 35, 551

Sui, C., Zhao, X., Jing, T., & Mao, Y. 2023, arXiv:2307.04994

Tejero-Cantero, A., Boelts, J., Deistler, M., et al. 2020, JOSS, 5, 2505

Tingay, S. J., Goeke, R., Bowman, J. D., et al. 2013, PASA, 30, e007

Trott, C. M. 2016, MNRAS, 461, 126

Trott, C. M., Jordan, C. H., Midgley, S., et al. 2020, MNRAS, 493, 4711

Valogiannis, G., & Dvorkin, C. 2022, PhRvD, 105, 103534

van Haarlem, M. P., Wise, M. W., Gunst, A. W., et al. 2013, A&A, 556, A2

Van Rossum, G., & Drake, F. L. 2009, Python 3 Reference Manual (Scotts Valley, CA: CreateSpace)

Van Rossum, G., & Drake, F. L., Jr 1995, Python Reference Manual (Amsterdam: Centrum voor Wiskunde en Informatica)

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Waskom, M. L. 2021, JOSS, 6, 3021

Weiler, M., Geiger, M., Welling, M., Boomsma, W., & Cohen, T. S. 2018a, in Advances in Neural Information Processing Systems, Vol. 31, ed. S. Bengio et al. (Red Hook, NY: Curran Associates, Inc.)

Weiler, M., Hamprecht, F. A., & Storath, M. 2018b, in 2018 IEEE/CVF Conf. on Computer Vision and Pattern Recognition (New York: IEEE), 849

Wolz, L., Abdalla, F. B., Alonso, D., et al. 2015, in Proc. Advancing Astrophysics with the Square Kilometre Array (AASKA14) (Trieste: SISSA), 35

Yoshiura, S., Pindor, B., Line, J. L. B., et al. 2021, MNRAS, 505, 4775

Zhao, X., Mao, Y., Cheng, C., & Wandelt, B. D. 2022a, ApJ, 926, 151

Zhao, X., Mao, Y., & Wandelt, B. D. 2022b, ApJ, 933, 236

Zhao, X., Ting, Y.-S., Diao, K., & Mao, Y. 2023a, MNRAS, 526, 1699

Zhao, X., Zuo, S., & Mao, Y. 2023b, arXiv:2307.09530

Zheng, H., Tegmark, M., Dillon, J. S., et al. 2017, MNRAS, 464, 3486

Ziegel, J. F., & Gneiting, T. 2014, EJSta, 8, 2619

Zonca, A., Singer, L., Lenz, D., et al. 2019, JOSS, 4, 1298