#  LtU-ILI: An All-in-One Framework for Implicit Inference in Astrophysics and
Cosmology

Matthew Ho1∗<EMAIL_ADDRESS>Deaglan J. Bartlett1 Nicolas Chartier2
Carolina Cuesta-Lazaro3,4,5 Simon Ding1 Axel Lapel1,15 Pablo Lemos6,7,8,9
Christopher C. Lovell10 T. Lucas Makinen11 Chirag Modi9,12 Viraj Pandya13,14
Shivam Pandey13 Lucia A. Perez9 Benjamin Wandelt1,9 Greg L. Bryan13 1 CNRS &
Sorbonne Université, Institut d’Astrophysique de Paris (IAP), UMR 7095, 98 bis
bd Arago, F-75014 Paris, France  2 Center for the Gravitational-Wave Universe,
Astronomy Program Department of Physics and Astronomy, Seoul National
University, Seoul 08826, Korea  3 Center for Astrophysics — Harvard &
Smithsonian, 60 Garden St, Cambridge, MA 02138, USA  4 The NSF AI Institute
for Artificial Intelligence and Fundamental Interactions  5 Department of
Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA  6
Department of Physics, Université de Montréal, Montréal, Canada  7 Mila -
Quebec Artificial Intelligence Institute, Montréal, Canada  8 Ciela - Montreal
Institute for Astrophysical Data Analysis and Machine Learning, Montréal,
Canada  9 Center for Computational Astrophysics, Flatiron Institute, 162 5th
Avenue, New York, NY 10010, USA  10 Institute of Cosmology and Gravitation,
University of Portsmouth, Burnaby Road, Portsmouth, PO1 3FX, United Kingdom
11Imperial Centre for Inference and Cosmology (ICIC) & Astrophysics Group,
Imperial College London, Blackett Laboratory, Prince Consort Road, London SW7
2AZ, United Kingdom  12 Center for Computational Mathematics, Flatiron
Institute, 162 5th Avenue, New York, NY 10010, USA  13 Columbia Astrophysics
Laboratory, Columbia University, 550 West 120th Street, New York, NY 10027,
USA  14 Hubble Fellow  15 Sorbonne Université, Université Paris Diderot,
Sorbonne Paris Cité, CNRS, Laboratoire de Physique Nucléaire et de Hautes
Energies (LPNHE). 4 place Jussieu, F-75252, Paris Cedex 5, France

###### Abstract

This paper presents the Learning the Universe Implicit Likelihood Inference
(LtU-ILI) pipeline, a codebase for rapid, user-friendly, and cutting-edge
machine learning (ML) inference in astrophysics and cosmology. The pipeline
includes software for implementing various neural architectures, training
schema, priors, and density estimators in a manner easily adaptable to any
research workflow. It includes comprehensive validation metrics to assess
posterior estimate coverage, enhancing the reliability of inferred results.
Additionally, the pipeline is easily parallelizable, designed for efficient
exploration of modeling hyperparameters. To demonstrate its capabilities, we
present real applications across a range of astrophysics and cosmology
problems, such as: estimating galaxy cluster masses from X-ray photometry;
inferring cosmology from matter power spectra and halo point clouds;
characterising progenitors in gravitational wave signals; capturing physical
dust parameters from galaxy colors and luminosities; and establishing
properties of semi-analytic models of galaxy formation. We also include
exhaustive benchmarking and comparisons of all implemented methods as well as
discussions about the challenges and pitfalls of ML inference in astronomical
sciences. All code and examples are made publicly available at
https://github.com/maho3/ltu-ili.

###### keywords:

cosmology, astrophysics, statistical methods, machine learning

##  1 Introduction

Statistical inference of unknown quantities is a fundamental problem in
science. In the astronomical sciences, we largely rely on Bayesian inference
to test new physical models (Feigelson and Babu, 2012; Dodelson and Schmidt,
2020; Eadie  _et al._ , 2023), wherein we assume some prior hypothesis of an
observed system and calculate constraints on model parameters which would be
consistent with our observations. For nearly a century, the practice of
building linear or perturbative physical models from first-principles allowed
us to make substantial progress towards understanding the universe. Yet, as
highlighted in the 2020 Decadal Survey (National Academies of Sciences,
Engineering, and Medicine, 2023), exploring the complex, nonlinear regime of
physical phenomena through data-driven inference promises significant gains in
constraining power. The large data volume of next-generation surveys, the
improvements in high-resolution simulations, and the rapid rise of machine
learning (ML) techniques have driven an explosion of inquiry into how one can
automatically and rapidly learn complex physical phenomena (Carleo  _et al._ ,
2019). Despite the growing enthusiasm, a major hurdle to the widespread
adoption of data science methods in astronomy is in their limited
accessibility to the general community, which when compounded by the rapid
pace of developments in the field, makes it challenging for individuals to
stay well-informed on best practices (Huppenkothen  _et al._ , 2023). There is
currently no widespread consensus within the field regarding the approach to
building a framework for implicit inference problems that is both reliable,
robust, and accessible.

Implicit Likelihood Inference (ILI, Cranmer  _et al._ , 2020; Marin  _et al._
, 2012), also known as simulation-based inference (SBI) and likelihood-free
inference (LFI), is an approach for learning the statistical relationship
between parameters and data. ILI is a mathematically rigorous way of framing
ML under the umbrella of Bayesian statistics, a language common to
astrophysicists and cosmologists since the turn of the century (Christensen
_et al._ , 2001). It aims to estimate a full posterior distribution over
parameters, including predictive uncertainties. It is best explained in
contrast to traditional explicit likelihood analyses, wherein one might write
down an analytic likelihood to represent the probability of observations given
a model and its parameters, and then use sampling methods such as Markov Chain
Monte Carlo (MCMC; Brooks  _et al._ , 2011) to sample from the posterior. In
ILI, one seeks to automatically learn the likelihood, i.e. to model the
relationship between model parameters and observations through training on
simulations. As such, ILI can be applied to infer parameters for any
phenomenon that can be simulated, without the need for assumptions about the
analytical form of the likelihood. Furthermore, ILI can accommodate complex
data cuts (provided that consistent filters are applied to both the data and
the forward model) something that can be hard to model in likelihood-based
approaches.

Given a prescribed training set of model parameters and observations or a fast
physical simulator, ILI can produce robust, tightly-constraining inference
without analytic likelihoods. In astronomy and cosmology, ILI has consistently
been shown to accelerate and improve upon traditional analyses such as
cosmological galaxy lensing and clustering (_e.g.,_ Jeffrey  _et al._ , 2021;
Makinen  _et al._ , 2021, 2022a; de Santi  _et al._ , 2023a; Hahn  _et al._ ,
2023), gravitational waves (_e.g.,_ Dax  _et al._ , 2021a; Cheung  _et al._ ,
2022), galaxy cluster mass estimation (_e.g.,_ Ho  _et al._ , 2022; de Andres
_et al._ , 2022), galaxy morphology (_e.g.,_ Walmsley  _et al._ , 2020; Ghosh
_et al._ , 2022), stellar streams (_e.g.,_ Hermans  _et al._ , 2021; Alvey
_et al._ , 2023), and exoplanets (_e.g.,_ Rogers  _et al._ , 2023; Aubin  _et
al._ , 2023).

Despite its popularity, robust applications of ILI have their challenges.
First, fully Bayesian models include capturing epistemic (modeling)
uncertainty, which requires marginalizing over model parameter uncertainty
(_i.e.,_ network weights and biases) and is intractable for large
architectures. However, approximations to this method using dropout
marginalization (Gal and Ghahramani, 2016), model ensembling (Lakshminarayanan
_et al._ , 2017) or stochastic weight averaging (Maddox  _et al._ , 2019;
Lemos  _et al._ , 2023) prove to be empirically sufficient surrogates. Second,
choices of modeling hyperparameters (_e.g.,_ preprocessing, architecture,
learning procedures, etc.) can greatly affect the quality of a trained ILI
model. This is generally addressed through model ensembling or hyperparameter
searches (_e.g.,_ Jin  _et al._ , 2019; White  _et al._ , 2021). Lastly, a
pervasive problem in ILI is model misspecification, wherein the simulator or
training catalog used to learn likelihoods is not representative of reality.
Methods such as SBI++ (Wang  _et al._ , 2023) deal with missing data and
outliers and present applications to galactic photometric redshift inference
(Modi  _et al._ , 2023). For a thorough discussion of caveats, see Section 6.
We argue that addressing these common problems is paramount to responsible
application of ILI, and testing them exhaustively necessitates a unified,
accessible framework.

We introduce the first version of the Learning the Universe Implicit
Likelihood Inference code (LtU-ILI), a pipeline for training ML models for
regressive parameter estimation. Given a labeled training set of observed data
and parameters, LtU-ILI trains neural networks to emulate posterior
probability distributions. The code also provides scientists with a diverse
toolkit for testing and validation. It was built under the following design
principles:

  * •

The code includes state-of-the-art neural architectures, validation tools, and
samplers for enabling ILI.

  * •

The pipeline is modular, easily customizable, and parallelizable for rapid
testing of design choices and hyperparameters.

  * •

The interface is accessible to non-specialists while being practical for high-
level production.

  * •

The methodologies automatically implement standard best practices for machine
learning (See Huppenkothen  _et al._ , 2023, for a review).

These attributes are paramount to make ILI a relevant tool for a wide range of
inference problems, whether in the context of exploratory analysis or
exhaustive empirical testing.

The novel contributions of this work are as follows:

  * •

We unify several leading implicit inference codes–sbi (Tejero-Cantero  _et
al._ , 2020), pydelfi (Alsing  _et al._ , 2019), and lampe (Rozet  _et al._ ,
2021)–under a single common interface, for the first time allowing rigorous
apples-to-apples comparisons of their performance.

  * •

We supply extensions to these codes, enabling them to utilize new types of
embedding networks (_e.g.,_ for image- and graph-like datasets) and providing
state-of-the-art validation metrics for enabling robust, maximally-informative
inference.

  * •

We provide an accessible, comprehensive, dual Jupyter and command-line
interface, for rapid research and development.

  * •

We demonstrate the abilities of ILI in an extensive selection of inference
problems in astrophysics and cosmology, and we provide these as public
benchmarks for future works.

At the time of publishing, LtU-ILI offers the broadest and most complete
toolbox of features among ILI software. Over sbi, we include the capacity for
custom data-loading procedures, additional neural density estimators (NDEs),
and adaptability to exotic embedding architectures (including sequential and
graph encoders). Over pydelfi, we allow the capacity for neural posterior and
ratio estimation. Over lampe, we provide a configuration-based interface to
training procedures, the ability to specify separate proposal and prior
distributions, and automated validation tests. We also include new
multivariate coverage tests such as TARP (Lemos  _et al._ , 2023) and
adaptive, automatic integration of data-loading, inference, and validation. We
also remark that, during the development of this work, nbi (Zhang  _et al._ ,
2023) was released as a neural posterior estimation toolkit focused on spectra
and light-curve analysis. Though we do not natively integrate nbi, we remark
that we encompass the functionality of nbi in our various neural posterior
estimation implementations and extend them with alternatives and additional
features.

In this software release paper, we provide background, details, and example
applications of the LtU-ILI pipeline. The code is publicly available on
Github111https://github.com/maho3/ltu-ili. In Section 2, we review the
theoretical foundations of ILI. In Section 3, we describe the structure of our
pipeline and the various configurations that users can interact with. In
Section 4, we apply LtU-ILI on several synthetic problems where the true
solution is known, to demonstrate and benchmark its performance relative to
traditional methods. In Section 5, we provide numerous diverse examples of
parameter inference in science using LtU-ILI. In 6, we discuss the important
considerations to take into account when designing an ILI pipeline, as well as
possible failure modes and relevant solutions. Lastly, we provide a summary of
our findings in Section 7.

##  2 Implicit Likelihood Inference

Figure 1: Demonstration of the differences between a likelihood (top-left),
posterior (top-right), prior (bottom-right), and joint distribution (bottom-
left) for an arbitrary one-dimensional inference problem. Here, the univariate
data $x$ and parameters $\theta$ are associated by a quadratic relationship
with Gaussian scatter. Given the observed data point $x_{o}$ (in black), our
goal is to recover posterior constraints on $\theta$. In an ILI framework, we
seek to capture knowledge of the joint distribution, $P(x,\theta)$, by
building neural-network emulators for the likelihood (in red) or the posterior
(in blue).

The goal of ILI is the same as any Bayesian inference problem: given a model
and observed data $\bm{x}_{\rm o}$, we wish to know the distribution of model
parameters $\bm{\theta}$ which are consistent with our data, _i.e.,_ the
posterior distribution $\mathcal{P}\left(\bm{\theta}|\bm{x}_{\rm o}\right)$.
The data and parameters can take many forms. In astronomy, we might consider
the observed data $\bm{x}_{\rm o}$ to be, for example, the observed X-ray
photon count of a galaxy cluster, the gravitational wave signal from a
collapsing black hole binary, or the full 2D sky projection of galaxy weak
lensing measurements. From this data, we might wish to infer parameters
$\bm{\theta}$ such as, respectively, the mass of a galaxy cluster, the mass
ratio and ellipticity of the black hole binary, or the cosmological parameters
governing dark matter and dark energy. For Bayesians, the posterior
$\mathcal{P}\left(\bm{\theta}|\bm{x}_{\rm o}\right)$ is the key to
constraining our belief in physical laws after observing new data $\bm{x}_{\rm
o}$.

According to Bayes’ theorem, the posterior is composed of

$\mathcal{P}\left(\bm{\theta}|\bm{x}\right)\propto\mathcal{P}\left(\bm{x}|\bm{\theta}\right)\
\mathcal{P}\left(\bm{\theta}\right),$ (1)

where $\mathcal{P}\left(\bm{x}|\bm{\theta}\right)$ is the likelihood, which
defines the probability of the data given the model and some parameters, and
$\mathcal{P}\left(\bm{\theta}\right)$ is the prior, which specifies our prior
belief on the true distribution of the parameters. A demonstration of the
differences of the prior, posterior, and likelihood in an inference problem is
given in Figure 1.

There are two approaches to modeling the terms in Equation 1: explicit
(theory-driven) analyses or implicit (data-driven) analyses. In a classical
explicit likelihood-based analysis, we assume an analytic form for the
likelihood, $\mathcal{P}\left(\bm{x}|\bm{\theta}\right)$, and substitute it
into Equation 1. Typically, this expression is derived from analytic theory or
simple simulations, coupled with strong assumptions on the statistical
distribution of noise or unknown physics. Sometimes, these statistical
assumptions are relatively simple, _e.g.,_ a Gaussian distribution for
continuous data (_e.g.,_ Tegmark  _et al._ , 2004; Christensen and Meyer,
2022; Porqueres  _et al._ , 2022) or a Poisson distribution for discrete
counts (_e.g.,_ Boese and Doebereiner, 2001; Braun  _et al._ , 2008; Tsaprazi
_et al._ , 2022). Although these may be appropriate distributions for many
situations, for a general inference problem we may not know which distribution
to choose or whether this choice is valid. Also, MCMC-like samplers scale
poorly with dimensionality and have bad mode coverage. This can cause problems
for high-dimensional inference problems, especially if the likelihood is not
differentiable with respect to the parameters. However, if we have access to a
set of examples of the data-parameter relationship (_i.e.,_ samples from the
joint distribution,
$(\bm{x},\bm{\theta})\sim\mathcal{P}(\bm{x},\bm{\theta})$), we can use
implicit likelihood analyses to learn the shape of a likelihood or posterior
automatically with full differentiability. This flexibility makes ILI an
extremely powerful tool for modeling complex or otherwise poorly-understood
physical phenomena.

ILI learns approximate posteriors, $\hat{\mathcal{P}}(\bm{\theta}|\bm{x})$, or
likelihoods from the distribution of example data-parameter pairs in a
training catalog, $\mathcal{D}_{\rm
train}:=\\{(\bm{x}_{i},\bm{\theta}_{i})\\}_{i=1}^{N_{\rm train}}$. These data
can take any form; they can be sensor readings in an experiment (_e.g.,_
Dingeldein  _et al._ , 2023), images of the night sky (_e.g.,_ Ntampaka  _et
al._ , 2019), or even compressed summaries of high-dimensional stochastic
properties (_e.g.,_ Modi  _et al._ , 2023). They can be gathered from a pre-
run simulation (_e.g.,_ Fluri  _et al._ , 2022), be simulated on-the-fly
(_e.g.,_ Alsing  _et al._ , 2019), or curated from manually-labeled
observations (_e.g.,_ Lanusse  _et al._ , 2018), but they must be
representative of the true underlying problem. Given these data-parameter
pairs, we wish to construct a model capable of doing inference on real data,
_i.e.,_ to evaluate the posterior at some observed value $\bm{x}=\bm{x}_{\rm
o}$.

In this section, we describe the theoretical foundations of LtU-ILI. We begin
by providing a simple introductory example for posterior inference with
Approximate Bayesian Computation (ABC; Section 2.1). We then explain how one
can vastly accelerate this process using neural networks to emulate
probability distributions (Section 2.2). We then discuss how to practically
frame machine learning in the context of Bayesian inference (Sections 2.2.1,
2.2.2 & 2.2.3). Lastly, we discuss auxiliary topics on utilizing active
learning to accelerate training (Section 2.3) and validating black-box models
to ensure robustness in real-world applications (Section 2.4).

###  2.1 Approximate Bayesian Computation

Approximate Bayesian Computation (ABC; Rubin, 1984) is the original basis of
ILI. Given a prior $\mathcal{P}(\bm{\theta})$, a simulator
($\bm{\theta}\rightarrow\bm{x}$), and an observation $\bm{x}_{\rm o}$, ABC can
construct a posterior. The general procedure of Rejection ABC is
straightforward: (1) first, draw candidate parameter values from the prior
$\bm{\theta}\sim\mathcal{P}(\bm{\theta})$, (2) forward-simulate these
parameter values to observations $\bm{x}\sim\mathcal{P}(\bm{x}|\bm{\theta})$,
(3) measure a distance metric between the candidate observations and the real
observation, $d(\bm{x},\bm{x}_{\rm o})$, and (4) if the distance is
sufficiently close to the real observation (_i.e.,_ $d(\bm{x},\bm{x}_{\rm
o})<\epsilon$ for small $\epsilon$), accept the candidate parameters
$\bm{\theta}$ as samples from the inferred posterior. ABC has been widely
adopted in statistical inference because it does not require explicit
knowledge of the likelihood and the accepted samples are guaranteed to
converge to the true posterior as $\epsilon\rightarrow 0$. When $\epsilon$ is
non-zero, the inferred posterior is guaranteed to be broader than the true
posterior, enabling conservative inference. Despite successful applications in
astrophysics and cosmology (_e.g.,_ Schafer and Freeman, 2012; Cameron and
Pettitt, 2012; Hahn  _et al._ , 2017), the ultimate problem with ABC methods
is that the acceptance rate vanishes exponentially as the dimensionality of
$\bm{\theta}$ increases, requiring thus significantly more simulations, which
can be alleviated with a larger $\epsilon$ at the cost of targeting a broader
posterior than the true posterior (Alsing  _et al._ , 2018). This limitation
makes ABC highly intractable in regimes where running simulations is
expensive. Occasionally ABC is used as a synonym for LFI/SBI/ILI, however
throughout the paper we exclusively use ABC to refer to the algorithm
discussed above.

###  2.2 Neural Density Estimation

In modern applications, implicit inference is made tractable by using machine
learning models to emulate conditional probability distributions, a procedure
called Neural Density Estimation (Papamakarios, 2019). As an example, consider
a situation wherein we seek to train a mixture density network (MDN; Bishop,
1994) to directly model the conditional distribution
$\mathcal{P}(\bm{v}|\bm{u})$ using our training set of data-parameter pairs
$(\bm{x},\bm{\theta})\sim\mathcal{D}_{\rm train}$. For this subsection only,
we have $(\bm{u},\bm{v})=(\bm{x},\bm{\theta})$ or
$(\bm{u},\bm{v})=(\bm{x},\bm{\theta})$ interchangeably. The goal is to build a
neural architecture $q_{\bm{w}}(\bm{v}|\bm{u})$ with weights $\bm{w}$ that
outputs a probability distribution over $\bm{v}$ which targets the conditional
probability $\mathcal{P}(\bm{v}|\bm{u})$. An MDN outputs probability
distributions by allowing the neural network to specify the statistical
parameters of a mixture density distribution. In the independent Gaussian
case, this means the neural network outputs a mean and a variance for each
component of the mixture density distribution, _i.e.,_

$q_{\bm{w}}(\bm{v}|\bm{u})=\frac{1}{N_{\rm c}\sqrt{2\pi}}\sum_{i=1}^{N_{\rm
c}}\frac{1}{\sigma_{i}(\bm{u};\bm{w})}e^{-\frac{(\bm{v}-\mu_{i}(\bm{u};\bm{w}))^{2}}{2\sigma_{i}^{2}(\bm{u};\bm{w})}},$
(2)

where $N_{c}$ is the chosen number of components of the mixture density
distribution, $\bm{w}$ are the learned weights and biases of the neural
network, and $\mu_{i}(\bm{u};\bm{w})$ and $\sigma_{i}(\bm{u};\bm{w})$ are the
neural network outputs representing the mean and standard deviation (or scale)
of the $i$-th Gaussian component as a function of the input $\bm{u}$ and the
weights.

To train such a network to properly emulate the conditional probability, all
we need to do is maximize the joint likelihood of our training data
$\mathcal{D}_{\rm train}$. This is equivalently and often more conveniently
expressed in the form of minimizing the negative log-probability loss
(typically log-likelihood or log-posterior, see sections 2.2.1 to 2.2.3)

$\mathcal{L}_{\rm MDN}:=-\mathbb{E}_{\mathcal{D}_{\rm train}}\left[\log
q_{\bm{w}}(\bm{v}|\bm{u})\right],$ (3)

where the expectation $\mathbb{E}_{\mathcal{D}_{\rm train}}$ is taken over all
data-parameter pairs in the training set $\sim\mathcal{D}_{\rm train}$. For
the case of the the posterior distribution
($\bm{v}=\bm{\theta},\bm{u}=\bm{x}$), if the training data is sufficiently
diverse and if the MDN is sufficiently flexible, then the minimization of
$\mathcal{L}_{\rm MDN}$ over the network weights $\bm{w}$ will converge
$q_{\bm{w}}(\bm{\theta}|\bm{x})$ to the true posterior
$\mathcal{P}(\bm{\theta}|\bm{x})$ (Papamakarios and Murray, 2016).

The MDN is just one example of an ever-increasing variety of neural density
estimators in the literature. Another very popular class of methods is
normalizing flows (Papamakarios  _et al._ , 2021), which output a flexible
conditional distribution defined via learnable, invertible transformations of
a base Gaussian distribution. Normalizing flows are often better equipped for
handling cases wherein the target distribution is non-Gaussian, but low-
dimensional and uni-modal. For any type of neural density estimator, the
training procedure always involves minimizing the negative log-probability of
training data, as in Equation 3.

Now we have established the tools for neural density estimation, we can apply
them to fitting components of Bayes’ theorem (Equation 1).

####  2.2.1 Neural Posterior Estimation (NPE)

The most straightforward method to do this is to train neural networks to
directly emulate the posterior distribution, also known as Neural Posterior
Estimation (NPE; Papamakarios and Murray, 2016; Greenberg  _et al._ , 2019).
As in the MDN example in the previous section, the loss function for NPE is
simply the negative log-posterior of our learned neural density estimator
$\hat{\mathcal{P}}(\bm{\theta}|\bm{x})$, and can be trained using a loss
function like that of Papamakarios and Murray (2016):

$\displaystyle\mathcal{L}_{\rm NPE}$
$\displaystyle:=-\mathbb{E}_{\mathcal{D}_{\rm
train}}\log\hat{\mathcal{P}}(\bm{\theta}_{i}|\bm{x}_{i})$ (4)
$\displaystyle=-\mathbb{E}_{\mathcal{D}_{\rm
train}}\log\left[\frac{p(\bm{\theta})}{\tilde{p}(\bm{\theta})}q_{\bm{w}}(\bm{\theta}|\bm{x})\right],$

where our neural posterior $\hat{\mathcal{P}}(\bm{\theta}_{i}|\bm{x}_{i})$ is
decomposed into a neural network output $q_{\bm{w}}(\bm{x},\bm{\theta})$ and a
weighting factor taken as the ratio of the assumed prior $p(\bm{\theta})$ to
the proposal prior $\tilde{p}(\bm{\theta})$. The proposal prior is
representative of the distribution of $\bm{\theta}$ in the training dataset
$\mathcal{D}_{\rm train}$ and is often taken to be equal to the assumed prior.
Note, that the normalization factor $\log
p(\bm{\theta})/\tilde{p}(\bm{\theta})$ cancels out when optimizing
$\mathcal{L}_{\rm NPE}$ over $\bm{w}$, but is still required for constructing
the posterior estimator $\hat{\mathcal{P}}(\bm{\theta}|\bm{x})$. This method
has the advantage that it allows for quick evaluation and sampling of the full
posterior at inference time. However, the disadvantage is that it requires
knowledge of the analytic form of the proposal prior $\tilde{p}(\bm{\theta})$
from which training data was sampled, which may not be accessible for
astrophysical training data (_e.g.,_ cluster properties in cosmological
simulations). Recently, Vasist  _et al._ (2023) used an NPE strategy to
constrain exoplanetary atmospheric models. NPE has also become a known tool in
Gravitational-Wave astronomy where analytical Compact Binary Coalescent (CBC)
waveforms play the role of the forward model (_e.g.,_ with DINGO from Dax  _et
al._ 2021a, and the constraints on GW150914 by Crisostomi  _et al._ 2023,
among many others).

####  2.2.2 Neural Likelihood Estimation (NLE)

An alternative method which circumvents the issue of fixed priors is to only
fit for the likelihood, $\mathcal{P}\left(\bm{x}|\bm{\theta}\right)$, via
Neural Likelihood Estimation (NLE; Alsing  _et al._ , 2018, 2019; Papamakarios
_et al._ , 2019). The loss function for NLE notably does not include the
assumed or proposal prior, and simply maximizes the global log-likelihood,

$\mathcal{L}_{\rm NLE}:=-\mathbb{E}_{\mathcal{D}_{\rm train}}\log
q_{\bm{w}}(\bm{x}|\bm{\theta}),$ (5)

wherein we emphasize the swapped arguments of $q$ as compared to Equation 4.
Multiplying a learned likelihood with an assumed prior results in a proxy
which is proportional to the posterior, _i.e.,_
$\hat{\mathcal{P}}(\bm{\theta}|\bm{x})\propto
q_{\bm{w}}(\bm{x}|\bm{\theta})p(\bm{\theta})$. Generating samples from the
posterior is then simply a matter of using this proxy in running Markov chain
Monte Carlo (MCMC; Robert  _et al._ , 1999) or Variational Inference (Blei
_et al._ , 2017). For example, in Metropolis-Hastings MCMC (Robert  _et al._ ,
2004), the learned likelihood and assumed prior can be used to calculate the
acceptance probability of a transition
$\bm{\theta}_{t}\rightarrow\bm{\theta}^{\prime}$ without knowledge of the full
normalized posterior, _i.e.,_

$\alpha:=\min\left(1,\frac{\mathcal{P}(\bm{\theta}^{\prime})\mathcal{P}(\bm{x}|\bm{\theta}^{\prime})q(\bm{\theta}^{\prime}|\bm{\theta}_{t})}{\mathcal{P}(\bm{\theta}_{t})\mathcal{P}(\bm{x}|\bm{\theta}_{t})q(\bm{\theta}_{t}|\bm{\theta}^{\prime})}\right),$
(6)

where $q(\bm{\theta}^{\prime}|\bm{\theta}_{t})$ is the proposal mechanism.
Note, the process of obtaining a posterior is identical to the explicit
likelihood-based approach, but instead of an analytic likelihood function, NLE
has essentially a “black box” function instead. The advantage of this method
is that one only needs to perform training once for a given model, at the cost
of additional sampling once data becomes available. For instance, Jeffrey  _et
al._ (2021) fit a density model for the likelihood of compressed summaries
from DES SV weak lensing maps. The Madminer tool (Brehmer  _et al._ , 2020)
for particle physics notably allows for NLE using ILI tools.

####  2.2.3 Neural Ratio Estimation (NRE)

As in NLE, Neural Ratio Estimation (NRE; Hermans  _et al._ , 2020) targets a
proxy that is proportional to the posterior. Instead of outputting a
conditional density estimate for the likelihood, NRE models instead target the
likelihood ratio,

$r(\bm{x},\bm{\theta}):=\frac{\mathcal{P}(\bm{x}|\bm{\theta})}{\mathcal{P}(\bm{x}|\bm{\theta}_{0})},$
(7)

equal to the ratio of the likelihood at a given point in parameter space to
the likelihood at a reference point, $\bm{\theta}_{0}$. This formulation is
convenient because it can be folded into an acceptance ratio like Equation 6
and also be cast as a classification problem (Cranmer  _et al._ , 2015). The
output of a classification network, $d_{\bm{w}}(\bm{x},\bm{\theta})\in[0,1]$,
can be considered equivalent to a likelihood ratio as

$\hat{r}(\bm{x},\bm{\theta})=\frac{d_{\bm{w}}(\bm{x},\bm{\theta})}{1-d_{\bm{w}}(\bm{x},\bm{\theta})}.$
(8)

The classification problem can be interpreted as quantifying whether observed
data $\bm{x}$ is consistent with $\bm{\theta}$. We can then train this neural
network model using the loss:

$\mathcal{L}_{\rm NRE}:=\mathbb{E}_{\mathcal{D}_{\rm
train}}\left[d_{\bm{w}}(\bm{x},\bm{\theta})+\mathbb{E}_{\bm{\theta}^{\prime}\sim
p(\bm{\theta})}\left[1-d_{\bm{w}}(\bm{x},\bm{\theta}^{\prime})\right]\right],$
(9)

wherein we associate a positive classification with assigning a data $\bm{x}$
to the correct $\bm{\theta}$, and a negative classification with those
assigning data to a $\bm{\theta}^{\prime}$ sampled from the prior. The main
benefit of NRE over NPE and NLE is that one need not specify an approximating
distribution such as MDNs or normalizing flows to emulate a posterior. The
complexity is purely limited in terms of the depth and design of the neural
architecture. Among others, Cole  _et al._ (2022) and Karchev  _et al._ (2023)
applied a variation of NRE (namely Truncated Marginal NRE, Miller  _et al._
2021) to CMB data and Supernova data, respectively. Additionally, Delaunoy
_et al._ (2020) and Bhardwaj  _et al._ (2023) use NRE for Gravitational Wave
parameter inference.

####  2.2.4 Choosing an Estimator

The optimal choice of NPE, NLE, or NRE is highly dependent on the properties
on the underlying problem, the dimensionalities of $\bm{x}$ and $\bm{\theta}$,
and the intended application. First, the ease of learning the shape of a
likelihood or posterior distribution can vary dramatically between different
physical problems. For example, the shape of the unimodal posterior in Figure
1 is considerably simpler than that of the bimodal likelihood, and thus NPE
will be easier to converge than NLE. However, this is the opposite for some
applications presented in Section 4, where the likelihood is simpler and thus
NLE converges much faster than NPE or NRE. In cases where both the posterior
and likelihood exhibit complicated shapes, NRE models are a strong option, as
they do not require an explicit choice of NDE. Secondly, a good rule of thumb
is that neural networks are easier to train when using high-dimensional input
to produce low-dimensional output than vice versa. Thus, if
$\operatorname{dim}(\bm{x})$ is large (_e.g.,_ for image- or graph-like data)
then NPE may work better than NLE models. Once again, NRE models are a nice
middle-ground as they take in both $\bm{x}$ and $\bm{\theta}$. Note, this
heuristic may be changing with the advent of diffusion generative models (Song
_et al._ , 2020), but those are not yet integrated into LtU-ILI. Lastly, it is
imperative to consider the downstream applications of the inference. NPE
models are much better for cases in which the inference must be repeated for
many test $\bm{x}_{\rm o}$ (_e.g.,_ estimating morphology for many JWST
galaxies), as it can be prohibitively expensive to run MCMC chains for NLE and
NRE for many test points. However, if the learned posteriors will simply be
injected into a broader hierarchical likelihood sampling (_e.g.,_ using
cluster masses to constrain cosmology), NLE or NRE would be a more natural
choice.

For a quantitative comparison of NPE, NLE, and NRE on machine learning
benchmarks, see Lueckmann  _et al._ (2021). In any case, the best way to
choose a method for a new problem is to try each and perform quantitative
comparisons, tests which are made considerably easier with the configuration
interface in LtU-ILI.

###  2.3 Sequential Learning

Thus far, we have phrased our problem in such a way that we first have a set
of pre-run simulations, and only later do we obtain the posterior distribution
given our observations. This can be the case when the simulations are
particularly expensive (_e.g.,_ $N$-body simulations in cosmology, or Earth
systems models for climate science _e.g.,_ from which Watson-Parris  _et al._
2021 built emulators). In these cases, we must be careful about any difference
between the proposal prior for the parameters from which the simulations were
run (_e.g.,_ a Latin hypercube) and the chosen prior for the parameters used
in the inference (_e.g.,_ a multivariate Gaussian). It may be the case that
this proposal prior is significantly wider than the posterior distribution
inferred once the real data have been observed, meaning that many of the
simulations used in training for NPE/NLE/NRE have been “wasted” as the
posterior has almost no support in the region of parameter space occupied by a
(potentially significant) fraction of the simulations.

One can instead adopt a multi-round inference approach to improve simulation
efficiency, referred to as Sequential NPE/NLE/NRE or SNLE/SNLE/SNRE. In this
method, we begin by running a fraction of the total desired simulations with
parameters sampled from the broad prior, conducting the first round of
training for the NPE/NLE/NRE model. This weak inference is then applied to
$\bm{x_{\rm o}}$ to derive a first estimate for the posterior. The algorithm
then identifies regions of high posterior density and conducts additional
simulations within these areas to produce a refined posterior estimate. This
can be repeated until a convergence criterion is met or the total simulation
budget is exhausted. This iterative refinement of the posterior around the
region of interest can lead to improved posterior estimates for a given
experiment. However, it is important to note that these results are
deliberately optimized for a singular observed data point. The sampled
simulations may not be appropriate for a different experiment, potentially
resulting in a less accurate estimate of the posterior distribution if no
further simulations are run.

Note that an alternative method for generating parameters of interest using
acquisition functions that minimize the expected uncertainty in the posterior
density approximation exists. This is based on the Bayesian Optimization for
Likelihood-Free Inference (BOLFI) framework (Gutmann  _et al._ , 2016) has
also been suggested and applied to cosmological data (Leclercq, 2018),
although this is not yet part of LtU-ILI.

###  2.4 Validation

The goal of model validation is to assess whether the posteriors
$\hat{\mathcal{P}}(\bm{\theta}|\bm{x}_{\rm o})$ learned in Section 2.2 will be
accurate and reliable when applied to new data. Explicitly, we want to
evaluate: (1) whether the learned posteriors are maximally constraining of
$\bm{\theta}$ given the observed data $\bm{x}_{\rm o}$, and (2) whether the
predictive uncertainties quoted by our learned model are accurately calibrated
to our training data. These two criteria are naturally adversarial, and is
often referred to as the bias-variance tradeoff. For example, a model can
satisfy condition (1) by reporting tight error bars, but it will then fail
condition (2) when its error bars are smaller than the true distribution of
training data. Similarly, a model that produces a posterior equal to the prior
will easily satisfy condition (2) but will ultimately be uninformative for
solving condition (1). Despite the ‘black box’ nature of neural networks,
these criteria control the quality of the learned implicit inference
posteriors, mitigating the chance that the learned model is not representative
of the data used to train it.

We evaluate the constraints and coverage of our learned posteriors on a
labeled test set of data-parameter pairs $\mathcal{D}_{\rm
test}:=\\{(\bm{x}_{i},\bm{\theta}_{i})\\}_{i=1}^{N_{\rm test}}$. This test set
is meant to represent a fully independent sampling of the real data
distribution, unseen by the model prior to validation. This test set may
originate from the same source as our labeled training set $\mathcal{D}_{\rm
train}$ but must be held independent to avoid overfitting and underestimation
of the uncertainty (see Huppenkothen  _et al._ 2023 for a review of best
practices). If the distribution of data-parameter pairs in the test set
matches that of the predictions of our learned posterior, then our model is
calibrated and can be reliably extended to new, unlabeled data.

We build tests by comparing the true parameters from the test set to posterior
samples derived from our models. Sampling-based comparisons provide a unified
framework for all NPE/NLE/NRE learning strategies since they do not require
normalization of the posterior PDF. The technical details for sampling from
neural networks depend on the implemented learning strategy and are discussed
in more detail in Section 3.3. Here, we will assume that, for a given input
$\bm{x}$, we can draw independent and identically distributed (i.i.d.) samples
from the learned posterior distributions
$\hat{\bm{\theta}}\sim\hat{\mathcal{P}}(\bm{\theta}|\bm{x})$. These samples
can then be used to compare the level of constraint on $\bm{\theta}$ and the
consistency with the true values.

To test precision, which reflects the ability of a learned posterior to
constrain parameter values, we examine the distribution and shape of posterior
samples $\hat{\bm{\theta}}$ around the true value $\bm{\theta}$. These
distributions are commonly visualized using multivariate corner plots or
marginal true vs. predicted plots. The contours in these plots help
qualitatively assess constraining power, identify posterior degeneracies, and
detect systematic biases.

A quantitative measure of precision is the cumulative likelihood of the test
dataset, denoted as $\hat{\mathcal{P}}(\mathcal{D}_{\rm
test})=\prod_{i=1}^{N_{\rm
test}}\hat{\mathcal{P}}(\bm{\theta}_{i}|\bm{x}_{i})$. Larger test likelihoods
indicate a model posterior that concentrates more probability mass around the
true value, implying greater constraining power. The test likelihood is simple
to calculate for NPE methods, which directly estimate
$\mathcal{P}(\bm{\theta}|\bm{x})$. However, for NLE and NRE methods,
evaluating the log-likelihood involves an additional step of training a
generative model $\mathcal{G}(\bm{\theta})$ from the inferred posterior
samples $\hat{\bm{\theta}}$. This process is computationally expensive,
especially for large datasets, and not sensitive to overconfidence.

Lastly, given samples from a reference posterior, we can evaluate the
discrepancy between a model’s prediction and its optimum. Reference posteriors
can often be attained via long-run MCMC chains for explicit likelihood
approaches, and then used to compare with ILI methods. Lueckmann  _et al._
(2021) performed an extensive quantitative analysis on various sample-wise
distances for benchmarking ILI methods and found that Classifier 2 Sample
Tests (C2ST; Lopez-Paz and Oquab, 2016) were the most constraining measure of
accuracy. C2ST is defined as the accuracy of a trained classifier (often a
neural network itself) to distinguish between the true and inferred posterior
samples. High C2ST suggests that the true and inferred posterior are very
different, whereas a C2ST of $\sim 0.5$ suggests they are nearly
indistinguishable. C2ST is later used to benchmark LtU-ILI in Section 4.1.

Posterior samples are also useful for quantifying the calibration of our model
uncertainty in parameter space. We can construct a direct comparison of the
predicted percentiles by our inference engine to the true error observed in
our validation dataset. We first define the Probability Integral Transform
(PIT; Cook  _et al._ , 2006) as the cumulative density function (CDF) of our
model posterior given an input $\bm{x}_{\rm o}$,

${\rm PIT}(\bm{\theta};\bm{x}_{\rm o})=\int_{-\infty}^{\bm{\theta}}\
d\bm{\theta}\ \hat{\mathcal{P}}(\bm{\theta}|\bm{x}_{\rm o}).$ (10)

Given i.i.d. samples from the posterior, we can construct an estimator for the
PIT value as

$\hat{\rm PIT}(\bm{\theta};\bm{x}_{\rm
o})=\mathbb{E}_{\hat{\bm{\theta}}\sim\hat{\mathcal{P}}(\bm{\theta}|\bm{x}_{\rm
o})}\left[\Theta\left(\hat{\bm{\theta}}-\bm{\theta}\right)\right],$ (11)

where $\Theta$ is the Heaviside step function. Stated plainly, the PIT value
counts the number of times that our posterior samples $\hat{\bm{\theta}}$ fall
below the true parameter value $\bm{\theta}$. If our model posteriors are
globally consistent with the truth, i.e. we match the true posterior
everywhere in data-parameter space, then the distribution of PIT values
evaluated on data-parameter pairs from our test set must be uniformly
distributed (Zhao  _et al._ , 2021):

${\rm PIT}(\bm{\theta}_{i};\bm{x}_{i})\sim U(0,1),\ \forall\
(\bm{x}_{i},\bm{\theta}_{i})\in\mathcal{D}_{\rm test}\,.$ (12)

As a result of this property, the distribution of PIT values is a common
goodness-of-fit metric for conditional density models (e.g. Cook  _et al._ ,
2006; Bordoloi  _et al._ , 2010; Tanaka  _et al._ , 2018). This is also known
as Simulation-Based Calibration (Talts  _et al._ , 2018).

To better interpret errors in the PIT distribution, we can use a percentile-
percentile (P-P) plot, explicitly comparing the CDF of PIT values to the CDF
of a uniform random variable. Intuitively, the P-P plot measures: ‘What
percentile level is my posterior model assigning to the true value, and with
what frequency does this occur in the test set?’ If the predictive posterior
is well-calibrated, these should be the same, _e.g.,_ we should predict the
true value below the 50-th percentile 50% of the time. If it is not, the shape
of the P-P plot acts as a sensitive probe of global bias or over-/under-
dispersion of predictive uncertainty (see Figure 1 of Zhao  _et al._ 2021 for
interpretations of error with P-P plots). This is an extremely practical tool
for understanding and correcting biases and over/under-dispersion in complex
ILI posteriors.

As the dimensionality of $\bm{\theta}$ increases, proper coverage of the
posterior distribution requires exponentially more samples. As a result, the
PIT estimator (Equation 11) becomes intractably difficult to measure with low
variance in the high-dimensionality regime. In this regime, we use approximate
methods to bound the constraints of posterior coverage tests. Firstly, we can
construct a PIT value from the marginal posterior over each $i$-th component
of $\bm{\theta}$, _i.e.,_

$\hat{\rm
PIT}(\theta_{i};\bm{x})=\mathbb{E}_{\hat{\bm{\theta}}\sim\hat{\mathcal{P}}(\bm{\theta}|\bm{x})}\left[\Theta\left(\hat{\theta}_{i}-\theta_{i}\right)\right],$
(13)

where $\hat{\theta}_{i}$ and $\theta_{i}$ are the $i$-th components of
$\hat{\bm{\theta}}$ and $\bm{\theta}$ respectively. If our model posterior is
globally consistent with the true posterior, then this marginal PIT value has
the same properties as Equation 11, _i.e.,_ when evaluated on the test set
they should be uniformly distributed.

It is also sufficient to use approximations to check multivariate posterior
coverage. The Tests of Accuracy with Random Points (TARP; Lemos  _et al._ ,
2023) method uses samples in high-dimensional parameter space to estimate
expected coverage probabilities. By examining the density of posterior samples
within regions of parameter space near true values, TARP constructs estimates
of posterior coverage which are guaranteed to converge to the true posterior
coverage with sufficient samples. With enough samples, this test is necessary
and sufficient to guarantee an optimal posterior estimate.

##  3 Code Structure

Figure 2: Overview of the LtU-ILI pipeline, displaying the procedural
processes when running each stage. Each of the three stages (Data, Inference,
Validation) is independently configurable, meaning any setup of one stage will
automatically link to the others.

The LtU-ILI software is a pipeline for training implicit inference models to
infer scalar parameters from observational data. A typical inference pipeline
consists of three main stages: Data, Inference, and Validation. The Data stage
involves loading datasets from file or setting up on-the-fly simulators (e.g.
for the methodologies in Section 2.3) in order to provide data-parameter pairs
for implicit inference. The Inference stage then takes these data-parameter
pairs, trains neural networks to connect them with amortized posterior or
likelihood representations (as in Section 2.2), and saves the fitted models to
file. Finally, the Validation stage loads a test dataset as well as the fitted
models from Inference in order to evaluate goodness-of-fit metrics (e.g.
constraint level, posterior coverage, etc.) and produces inference on
unlabeled data. A schematic of this pipeline is shown in Figure 2.

Each stage is independent and adaptive, such that the configuration of one
stage can be changed while keeping the others the same and the pipeline will
still run. For example, if you had a fixed dataset but wanted to try out many
different neural architectures, you would only need to change the Inference
configuration. Similarly, if you wanted to see the gain in constraining power
of using a short, conservative data vector versus a longer, more informative
data vector, you could simply change the Data configuration to load the
different datasets and the Inference and Validation would proceed exactly the
same way. This design principle is extremely practical for doing apples-to-
apples comparisons of the many design choices and hyperparameters associated
with ILI, a key tenet of robust and responsible machine learning (Huppenkothen
_et al._ , 2023).

The pipeline is designed to promote exploratory research while also supporting
efficient production-level testing. We provide a fully functional Jupyter
interface and detailed examples to access each stage of LtU-ILI. The Jupyter
interface is a widely-used tool for examining new problems, building
educational material, and facilitating collaborative efforts (Wang  _et al._ ,
2019). The Jupyter interface of LtU-ILI provides step-by-step guidance for
establishing and testing a new inference pipeline that is familiar to novices
and experts alike. This allows LtU-ILI to act both as a user-friendly toolkit
to explore inference on a new problem as well as carry out thorough testing
and hyperparameter tuning through its yaml-based configuration API. Users can
generate many configurations of the pipeline and run multiple fittings in
parallel on large-scale GPU compute clusters.

In this section, we discuss the design of the Data, Inference, and Validation
stages of the LtU-ILI pipeline, including each of their capabilities and
relevant reference works. For specifics on using the API, we refer the reader
to our full code documentation222https://ltu-ili.readthedocs.io.

###  3.1 Data

The objective of the Data stage is to gather data-parameter pairs and present
them in an amenable format to the Inference stage. For standard inference,
with fixed training and test catalogs, the Data stage simply loads data-
parameter pairs from in-memory or on-disk storage and provides get_all_data()
and get_all_parameters() helper functions to pass along data and parameters in
the correct format. These functions are called during the Inference stage for
later pre-processing within the sbi, pydelfi, and lampe backends. Within LtU-
ILI, we provide several convenience classes to prepare data from either
numpy-, xarray-, or torch-like data formats. Users are also can create their
own data-loading objects, so long as they adhere to the aforementioned
conventions.

For multi-round inference (Section 2.3), the Data stage also provides a
simulate() function, wherein the loader will execute an on-the-fly simulation
for given parameters to generate data for efficient sequential learning. In
the example class provided in LtU-ILI, users can initialize this stage with an
arbitrary simulation function that intakes parameters $\bm{\theta}$ and
outputs data $\bm{x}$, and the sequential Inference stage will handle any
calls to the function. Also, users can load a pre-run initial catalog from
file and write any newly-generated simulations to disk during the training
procedure.

###  3.2 Inference

The LtU-ILI code encompasses a broad diversity of methodologies and
architectures for performing implicit inference. However, the core
configurations for all types of inference are universal. Every Inference stage
requires specifying: a posterior-, likelihood-, or ratio-estimation; one or
multiple neural architectures; a prior distribution; and training
hyperparameters (e.g. learning rate, batch size, etc.).

Once these are provided, the Inference stage knows how to pull data-parameter
pairs from the Data stage, perform data preprocessing, construct each training
algorithm, minimize losses (see Equations 4, 5, and 9), and produce a learned
posterior model, $\hat{\mathcal{P}}(\bm{\theta}|\bm{x})$. Each inference
algorithm automatically implements best practices for training, including data
normalization (Singh and Singh, 2020), adaptive stochastic gradient descent
(Kingma and Ba, 2014), and validation-guided early stopping (Bai  _et al._ ,
2021). The remainder of this section is dedicated to the core considerations
of configuring the Inference stage.

The Inference stage is built on three distinct computational backends: pydelfi
(Alsing  _et al._ , 2019), sbi (Tejero-Cantero  _et al._ , 2020), and lampe
(Rozet  _et al._ , 2021). pydelfi is a Tensorflow package (Abadi  _et al._ ,
2016) for performing NLE using active learning and data compression
techniques. sbi is a PyTorch package (Paszke  _et al._ , 2019) for performing
NPE, NLE, and NRE, as well as their sequential analogs. lampe is also a
PyTorch package focused on NPE methods, with substantial support for diverse
embedding networks and normalizing flow models. The implementations of NLE in
pydelfi and sbi are similar in theory (Papamakarios  _et al._ , 2019) but
differ in applications of embedding networks. pydelfi has the option to train
in a two-step process, wherein first it trains an embedding network to
optimally compress data via Fisher information maximization (Charnock  _et
al._ , 2018; Alsing and Wandelt, 2018) and then attaches it to the head of an
NLE model, whereas sbi exclusively uses NLE embedding networks to compress
parameters. pydelfi’s two-step compression-estimation process has been shown
to significantly improve simulation-efficiency during sequential training for
applications in cosmology (Alsing  _et al._ , 2018). The lampe NPE in LtU-ILI
is a custom implementation of Lueckmann  _et al._ (2017), similar to SNPE-B of
sbi but allowing for greater variety of configuration, embedding networks, and
NDEs. LtU-ILI is a unifying framework for all of these mutually-exclusive
backends, allowing, for the first time, convenient apples-to-apples comparison
of each implementation of ILI.

The Tensorflow and PyTorch backends have a variety of pre-implemented NDEs and
ratio classifiers. The NDEs are either Mixture Density Networks (MDNs; See
Section 2.2) or normalizing flows (See Papamakarios  _et al._ , 2021, for a
review). The pydelfi architectures are implemented from scratch, whereas the
sbi flow models utilize nflows (Durkan  _et al._ , 2020) and lampe models
utilize zuko (Rozet  _et al._ , 2022). Both backends also accept various
standard and customizable prior distributions. The primary requirement of a
custom prior is the ability to evaluate $\log p(\bm{\theta})$ (for NDE models)
or a proportional proxy (for NLE/NRE models). If the assumed prior’s support
is bounded (_e.g.,_ a top-hat prior), the code internally implements an affine
parameter transformation to avoid assigning a nonzero probability to regions
of parameter space beyond prior bounds.

A key component of LtU-ILI is model ensembling, wherein we independently train
multiple neural networks on the same dataset and combine their predictions to
improve robustness to overfitting and model uncertainty. Singular NDEs are
prone to overconfidence (Hermans  _et al._ , 2022), _i.e.,_ they have a
tendency to underestimate the predictive uncertainty. Deep ensembling
(Lakshminarayanan  _et al._ , 2017) is a well-tested schema for correcting
overconfidence by averaging NDE predictions from multiple models to inflate
error bars and correctly capture model uncertainty. It is an alternative to
the computationally-expensive Bayesian neural networks (Blundell  _et al._ ,
2015; Gal and Ghahramani, 2016; Cobb and Jalaian, 2021), in which the full
weight posterior is learned during training and sampled during inference. We
recommend that users utilize deep ensembles to ensure robust inference within
LtU-ILI.

The adaptability of neural networks to various data inputs has greatly
contributed to their growth in astronomy and cosmology. The convenient
implementations of the various NDEs and classifiers in LtU-ILI can be further
augmented with customizable embedding networks. These embedding networks can
use classic Keras-like neural layers (Chollet  _et al._ , 2015) in both
Tensorflow and PyTorch, greatly simplifying implementation. The integration
with embedding networks in the lampe backend allow for exotic inputs, like
graphs and sequences. In the code, we provide several examples of embedding
networks, including connections to convolutional and graph neural networks
(LeCun  _et al._ , 1998; Kipf and Welling, 2016). We also note that all
backends can be run on GPUs, to take advantage of computational speed-ups when
using complex embedding networks.

###  3.3 Validation

The Validation stage provides modular functions to implement the validation
metrics described in Section 2.4 on the learned posteriors
$\hat{\mathcal{P}}(\bm{\theta}|\bm{x})$. The execution of the Validation stage
follows:

  1. 1.

Load the learned neural posteriors from the Inference Stage.

  2. 2.

Sample the posterior
$\hat{\bm{\theta}}\sim\hat{\mathcal{P}}(\bm{\theta}|\bm{x})$ at test data
points, $\mathcal{D}_{\rm test}$, or at an unlabeled observational point,
$\bm{x}_{\rm obs}$.

  3. 3.

Construct the validation metric and return or save it to file.

The most impactful configuration of the Validation stage is the choice of
sampling methods. LtU-ILI provides three classes of samplers: MCMC,
variational inference (VI), and direct sampling. All models (NPE/NLE/NRE) in
both backends are natively integrated with the emcee (Foreman-Mackey  _et al._
, 2013) sampling package. emcee implements the Goodman and Weare (2010) MCMC
sampler and improves efficiency through parallelization and affine
transformations but does not take advantage of probability gradients. The sbi
models also have access to PyTorch’s pyro samplers (Bingham  _et al._ , 2019),
including slice samplers (Neal, 2003), Hamiltonian Monte Carlo (HMC; Neal  _et
al._ , 2011), and the No-U-Turn sampler (NUTS; Hoffman  _et al._ , 2014). All
emcee and pyro samplers have a unified interface within LtU-ILI, wherein users
can simply specify their choice of sampler as well as the number, length, and
thinning of MCMC chains. Also, we provide the functionality for sbi models to
fit and sample from NDEs through neural VI (Graves, 2011). VI is a fast,
approximate alternative to MCMC sampling, particularly useful in the high-
dimensional parameter regime. Lastly, solely for the NPE models in sbi, users
can utilize direct sampling, which is faster than both MCMC and VI. We remark
that the choice of sampling configuration might change according to the
validation metric in question, so LtU-ILI allows for the flexibility to
specify this on a case-by-case basis.

##  4 Synthetic Experiments

In this section, we demonstrate the capabilities of LtU-ILI on a number of
synthetic experiments for which an analytic form for the posterior or
likelihood is known. We use these traditional experiments to measure benchmark
performance the variety of inference engines available in the code against
each other and classical methods (_e.g.,_ ABC, HMC).

###  4.1 Toy Problem

Figure 3: Example SNPE posterior inference in comparison to classical implicit
(Rejection ABC) and explicit (HMC) likelihood inference for the known toy
simulator in Equation 14. The true value for each parameter is shown in red
and contours are shown at the central $68\%$ and $95\%$ confidence intervals.
Note that the prior for each $\theta_{i}$ is a standard normal, as
$p(\theta_{i})=\mathcal{N}(0,1)$.

We first build a toy model to demonstrate LtU-ILI’s ability to learn non-
linear posteriors. We consider a 10-dimensional data vector $\bm{x}$
constructed element-wise from the following stochastic simulator:

$x_{i}=3\sin(k_{i}+\phi_{0})+\phi_{1}k_{i}^{2}+\epsilon_{i},$ (14)

where $k_{i}=(2i/3)-3$ and $i\in[0,9]$. Here, the data vector
$\bm{x}:=\\{x_{i}\\}_{i=0}^{9}\in\mathbb{R}^{10}$ is composed of a sinusoidal
signal with phase $\phi_{0}$, a quadratic signal of amplitude $\phi_{1}$, and
an i.i.d. noise component $\epsilon_{i}\sim\mathcal{N}(0,1)$. We further
decompose ${\bm{\phi}}$ into our target parameters
$\bm{\theta}\in\mathbb{R}^{3}$ via $\phi_{0}=\theta_{0}+\theta_{1}$ and
$\phi_{1}=\theta_{1}-3\theta_{2}^{2}$. We expect our inference engine to
forward constraints on latent variables $\bm{\phi}$ onto target variables
$\bm{\theta}$. We assume a standard normal prior on each of our parameters of
interest, $p(\theta_{i})=\mathcal{N}(0,1)$. This is a challenging inference
task because the likelihood is highly non-linear and the posteriors exhibit
strong degeneracies. However, learning parameters from such complex,
correlated data vectors is a common problem in astronomy and cosmology
(_e.g.,_ Dax  _et al._ , 2021a; Modi  _et al._ , 2023; Bartlett  _et al._ ,
2023).

Given this simulator, we train sbi-SNPE using an ensemble of two normalizing
flows (MAFs) to produce the example $\bm{\theta}$ posterior in Figure 3, shown
relative to equivalent constraints with Rejection ABC and long-run HMC. We
train the SNPE sequentially using 10 rounds of 2000 simulations each,
following the procedure described in Section 2.3 and using Equation 14. For
ABC, we run the same number of simulations (20,000) and take the closest
$0.5\%$ to our data vector, using an L2 distance. For HMC, we run eight chains
each with 10,000 tuning steps and 10,000 sampling steps, for a total of
160,000 simulations. This run has an effective sample size of $~{}3,900-8,000$
and a Gelman-Rubin $\hat{R}$ statistic of $1.001$ for all parameters. We
consider these long-run HMC samples to be equivalent to a reference ‘ground
truth’ posterior.

Despite the same simulation budget as ABC, the SNPE posterior constraints have
converged to the reference posterior and are strongly more constraining than
the ABC constraints. In both the SNPE and HMC posteriors, we clearly observe
the expected linear degeneracy between $\theta_{0}$ and $\theta_{1}$ as well
as the quadratic degeneracy between $\theta_{1}$ and $\theta_{2}$, direct
results from the aforementioned definitions of $\phi_{0}$ and $\phi_{1}$.
These two degeneracies imply another quadratic degeneracy between $\theta_{0}$
and $\theta_{2}$. Despite this non-linearity, we see the multi-dimensional 68%
and 95% confidence intervals are entirely consistent between SNPE and HMC.
However, ABC methods struggle in this regime. Because the data dimensionality
is ten, ABC requires orders-of-magnitude more samples to fully capture the
tails of the posterior distribution. This regime demonstrates the utility of
having a low-cost amortized posterior model like SNPE, especially in further
cases where the simulators are not as simple.

###  4.2 Benchmarking

Figure 4: Prediction error as a function of simulation budget for various LtU-
ILI training methodologies on the benchmark SLCP inference problem
(Papamakarios  _et al._ , 2019). For clarity, plotted points for different
models at the same simulation budget are slightly offset horizontally. The
Rejection-ABC method is shown on all subplots as a baseline. Error is defined
in terms of the Classifier 2-sample Test (C2ST; Lopez-Paz and Oquab, 2016)
metric. A lower C2ST indicates a more accurate inferred posterior, with 0.5
being the optimal score. C2ST values are shown at their median and central 95%
confidence interval, calculated over ten independent runs.

We benchmark all models in LtU-ILI on the standard “Simple Likelihood Complex
Posterior” problem (SLCP; Papamakarios  _et al._ , 2019) common in simulation-
based inference experiments (_e.g.,_ Greenberg  _et al._ , 2019; Lueckmann
_et al._ , 2021; Ramesh  _et al._ , 2022). SLCP has a fairly straightforward,
unimodal likelihood but a very complex, multimodal posterior distribution,
making it a difficult, but attainable benchmark for implicit inference. We
follow the implementation of SLCP in Lueckmann  _et al._ (2021), wherein we
train various NPE/NLE/NRE methods and compare their resultant posteriors to
reference posteriors determined with long-run HMC chains. We implement this
for all available NPE/NLE/NRE methods and their sequential analogs in each of
the pydelfi, sbi, and lampe backends. We also include benchmarking of
Rejection ABC as a baseline. Throughout all backends and engines, we use the
same architectures as in Lueckmann  _et al._ (2021), which were determined
through exhaustive hyperparameter searches.

For each trial, we measure the error of an inference engine with respect to a
reference posterior using C2ST (see Section 2.4). For each backend, engine,
and simulation budget, we perform independent training and inference on ten
separate reference posteriors to aggregate test statistics.

In Figure 4, we show the performance of each LtU-ILI engine on the SLCP
benchmark as a function of the simulation budget used for generating the
training set. As we would expect, the accuracy of inferred posteriors
increases with the number of given simulations across all models. We observe
that all models strongly outperform the Rejection ABC baseline due to their
efficient use of simulation data. For the SLCP benchmark, we observe that NLE
methods perform the best, NPE methods are a close second, and NRE methods have
the worst performance. However, we note that the SLCP by definition has a
simple, easy-to-learn likelihood, encouraging strong NLE performance. This
hierarchy of performance could be entirely different for other problems.
Sequential methods (SNPE/SNLE/SNRE) show very little improvement for small
simulation budgets ($10^{3}$), but show massive gains for large budgets
($10^{5}$), especially for the NLE methods. This is likely because, for small
training datasets, the posteriors are poorly constrained and barely deviate
from the prior. However, for the tight posterior constraints in large budget
runs, sequential methods are able to more easily focus simulation resources on
small regions of parameter space. Next, we note that the lampe-NPE
implementation largely mirrors that of sbi-NPE, and pydelfi-NLE slightly
outperforms sbi-NLE for a medium simulation budget ($10^{4}$). Although
hyperparameters (_i.e.,_ architectures, learning rate, validation schema) are
fixed for these experiments, slight differences in backend implementations
might lead to better or worse performance. Lastly, as a check, the
performances of all sbi engines closely match those found in Lueckmann  _et
al._ (2021), suggesting our implementation of SLCP is correct.

##  5 Science Experiments

In this section, we demonstrate various applications of LtU-ILI for common
astrophysics and cosmology problems. These experiments are designed to display
the diverse functionality of the code, including using CNN embedding networks
(Section 5.1), applying posterior coverage validation (Section 5.2),
incorporating graph neural network layers (Section 5.3), using multiround
inference (Section 5.4), quantifying information in different data sources
(Section 5.5), and integrating information-optimal pretrained embeddings
(Section 5.6). We note that these experiments are intended to illuminate
interesting pathways for the scientific application of LtU-ILI and not
necessarily make significant scientific conclusions on their subject material.
For readability, experiment descriptions are kept short, but the code and data
are made publicly available through the codebase.

###  5.1 X-ray Mass Estimation of Galaxy Clusters

Figure 5: Inference of galaxy cluster $M_{500c}$ mass from eROSITA-like mock
X-ray observations, as in Ho  _et al._ (2023). The left four subplots show
random examples of single-band X-ray observables in terms of their sky-
projected surface brightness and target mass (in $h^{-1}M_{\odot}$). The right
plot is the true vs. predicted log-mass estimates on the test set, where we
show the median and central 68% confidence interval of our neural posteriors.
For readability, test points in the right plot have been randomly subsampled.

A prime functionality of LtU-ILI is its flexibility to adapt to many
modalities of data. Through the use of a customizable embedding network, we
process image data to estimate the mass of galaxy clusters from X-ray
observations. Our inference problem mirrors that of Ho  _et al._ (2023),
wherein we try to estimate the $M_{\rm 500c}$ mass of a galaxy cluster from
$128\times 128$ images of sky-projected X-ray photon counts. To train and test
our model, we use a dataset of $3,285$ mock X-ray observations of clusters
derived from the Magneticum hydrodynamical simulation (Dolag  _et al._ , 2016)
designed to mimic those of the eROSITA telescope (Soltis  _et al._ , 2022).
These mock observations are single-band images of the X-ray photon emission
from the hot intra-cluster medium (ICM) and include realistic simulation of
the contaminating systematics which afflict real X-ray measurements, including
cluster morphology, background emission, telescope response, and AGN sources.
Examples of these observables are shown in Figure 5. A good estimator of
cluster mass needs to be able to disentangle the source X-ray signal from the
noise and peripheral emission, while also understanding the physical
connection between X-ray emission, ICM gas content, and system mass.

To recreate this example, we perform NPE using an embedding network of the
same convolutional architecture as Ho  _et al._ (2023). That is, instead of
the last dense layer mapping to a single point estimate, we forward the final
embedding into an NDE. We construct an ensemble of four models, each with the
same embedding architecture, but two with Masked Autoregressive Flow (MAF
Papamakarios  _et al._ , 2017) NDEs and two with MDNs. We then train the
models from scratch using the NDE losses, instead of the mean-squared-error
(MSE) loss of Ho  _et al._ (2023). As in the reference, we assume a uniform
prior on cluster mass. We split our mock catalog into 90% training data and
10% test data, and present our performance results on the latter. The entire
training and testing procedure takes about fifteen minutes on an Nvidia V100
GPU.

The true vs. predicted mass estimates for this model are shown in Figure 5.
Using the original architecture, we achieve a test scatter on the mean
prediction of $0.0782\ {\rm dex}$. When compared to the $0.0773\ {\rm dex}$
scatter of the single-band X-ray images in Ho  _et al._ (2023), we achieve a
very similar level of information extraction, though our uncertainty is
slightly higher than the original. We attribute this to the fact that model
ensembling slightly inflates predictive uncertainty (Hermans  _et al._ ,
2022). However, we remark that by slightly increasing the number of
convolutional filters in the first two layers, we observed mean scatters as
low as $0.071\ {\rm dex}$, suggesting that the NPEs are capable of state-of-
the-art constraints on physical problems when given the right architecture. We
also note that hyperparameter search through neural architectures is made
efficient and accessible through LtU-ILI’s configuration-based interface.

###  5.2 Quijote Dark Matter Power Spectrum

(a) True vs. Predicted (b) Marginal Percentile Coverage Test

(c) TARP Multivariate Coverage Test

Figure 6:  Inference of three cosmological parameters $\Omega_{m}$, $h$, and
$\sigma_{8}$ from halo power spectrum multipoles in Quijote simulations using
SBI. In (a) we show the comparison between true and predicted value of the
parameters in the test simulations. Marginal error bars are shown at the
central $68\%$ confidence intervals. In (b) and (c) we show the coverage tests
using the P-P plots as well as using TARP method. These show bootstrapping
error bars at $68\%$ and $95\%$, indicative of typical variations in an ideal
estimator for a finite test set.

One of the major applications of SBI in the field of cosmology is to infer the
cosmological parameters from observations. This takes the form of using
simulations as a forward model, measuring summary statistics that capture the
relevant cosmological information, and then comparing them to observations.
Here we demonstrate this application by predicting three cosmological
parameters using the power spectrum multipoles measured from the halo catalog
of Quijote simulation (see Villaescusa-Navarro  _et al._ , 2020, for a
description of the simulations). We use the power spectrum multipoles
corresponding to $\ell=0,1,2$ in 23 linearly spaced wavenumber ($k$) bins
ranging from 0.08$h/{\rm Mpc}$ to 2.8$h/{\rm Mpc}$ as our observations, $x$.
We use the SNLE method to learn the likelihood $p(x|\theta)$, where the
parameters $\theta$ correspond to three cosmological parameters: the total
matter density $\Omega_{m}$, amplitude of the matter fluctuations
$\sigma_{8}$, and the Hubble constant $h$ measuring the rate of the expansion
of the Universe.

We train an ensemble of six networks, each using the MAF architecture with
three transformations and 10 hidden features. We use 1800 simulations to train
the network, retaining the remaining 200 to test its performance as shown in
Figure 6. We then use VI sampling to obtain the posterior of the parameters
from this trained likelihood ensemble, weighted by their validation losses.

We show the true vs predicted value of the three parameters on the top-left
panels (Figure 6a), finding good predictive performance for $\Omega_{m}$ and
$\sigma_{8}$. The lack of predictability for $h$ is expected physically, as
the power spectrum multipoles alone are significantly less sensitive to
expansion rate without some external prior information. We also perform
coverage tests on the posteriors, finding uniformly distributed rank
histograms as given by the coverage test in the bottom-left panels (Figure
6b). Finally, we show the coverage test using the TARP method in the center-
right panel (Figure 6c), finding that the posterior is correctly calibrated.

###  5.3 Field-level Inference with Graph Neural Networks

We next demonstrate LtU-ILI’s ability to handle different data structures by
also including a point cloud dataset example. In particular, we train a neural
posterior estimator to recover the cosmological parameters $\Omega_{m}$ and
$\sigma_{8}$ from the positions of the $10,000$ heaviest dark matter haloes in
the Quijote simulations (Villaescusa-Navarro  _et al._ , 2020). The high
resolution information present in the discrete halo positions has been shown
to constrain cosmological parameters to extremely high precision (_e.g.,_
Makinen  _et al._ , 2022b; Cuesta-Lazaro and Mishra-Sharma, 2023; de Santi
_et al._ , 2023b).

As an embedding network, we use a graph neural network. We construct a graph
where dark matter haloes act as nodes connected to any other dark matter halo
separated by a distance smaller than $20\,{\rm Mpc}/h$.

We train an ensemble of four networks, each one using a MAF architecture with
five transformations and $50$ hidden features. The graph neural network is a
message-passing neural network with node and edge neural networks being two
multi-layer perceptrons with $3$ layers of $128$ hidden units each. The
network takes as inputs 3D distances between nodes, and their corresponding
modulus.

In Figure 7, we show an example point cloud together with the resulting graph
neural network inference. The graph neural network obtains well-calibrated and
informative posteriors, consistent with true values. We note that the
constraints are not as tight as those produced for the Quijote dark matter
power spectrum in Figure 6. This is likely because the graph only targets the
most massive halos, _i.e.,_ a more physically realistic probe, while the
results in Figure 6 access the full, high-resolution dark matter distribution.

Figure 7: Cosmological inference from dark matter halo point clouds. Left:
Example point cloud from the Quijote simulations. The graph neural network
takes as an input the distances between halos separated by less than $20\,{\rm
Mpc}/h$. Right: Inferred posterior over $\Omega_{m}$ and $\sigma_{8}$ on a
test set simulation. True values are shown in gray. Contours are shown at the
central $68\%$ and $95\%$ confidence intervals.

###  5.4 Multi-round Inference for Gravitational Waves

(a) Example of GW signal in source and detector frame. (b) Multi-round
training and validation posterior probability.

(c) Posterior corner plot for SNPE inference.

Figure 8: Example parameter inference of a GW150914-like event as observed by
the gravitational wave Hanford (H1) and Livingstone detectors. (a) Example
simulated gravitational signal. Top: Radiation polarizations $h_{+}$ and
$h_{\times}$ in source frame. Bottom: Real part of the strain and noise
realization on Hanford (H1) detector frame. (b) Evolution of the log-
probability for the multi-round GW inference example using SNPE, _i.e.,_ the
negative of the loss function in equation (4) for $n_{r}=5$ rounds. (c)
Resultant posterior for a single run of the multi-round SNPE inference on GW
strain data. Contours are shown at $68\%$ and $95\%$ confidence intervals, and
dashed lines indicate the fiducial parameter values. The fake observation
$\bm{x}_{0}$ was generated with the same simulator function as the training
forward model.

Over recent years, several studies based on ILI methods have addressed the
problem of reconstructing parameters of merger events from gravitational wave
(GW) signals (Cutler and Flanagan, 1994; Thrane and Talbot, 2019). Initial
advancements were marked by significant contributions to posterior estimation
using normalizing flows, as highlighted in Green  _et al._ (2020); Green and
Gair (2020). Subsequent research (DINGO: Dax  _et al._ , 2021a) introduced
variations of NPE taking advantage of Group Equivariant Neural Posterior
Estimation (GNPE: Dax  _et al._ , 2021b) or Flow Matching Posterior Estimation
(Wildberger  _et al._ , 2023). More recently, Crisostomi  _et al._ (2023)
showcased a similar application of multi-round NPE to estimate ringdown
parameters of GW150914. Delaunoy  _et al._ (2020) and Bhardwaj  _et al._
(2023) also serve as recent illustrations of applying NRE to GW analyses – the
latter incorporating multi-round inference as well.

In this example, we explore the application of multi-round inference to
constrain a reduced set of gravitational wave (GW) parameters. Our forward-
modeled data vector is constructed by merging the simulated frequency-domain
strain signals from each of the Hanford (H1) and Livingstone (L1) detectors.
We employ the IMRPhenomPv2 algorithm (Khan  _et al._ , 2016) with a frequency
resolution of $\Delta f=0.125\textrm{ Hz}$, resulting in an input data vector
of $\textrm{dim}(\bm{x})=2\times 7872$. An illustration of this data vector is
presented in Figure 8a for clarity. We stress that, in contrast to likelihood-
based sampling methods implemented in Bilby (Ashton  _et al._ , 2019; Romero-
Shaw  _et al._ , 2020), which only necessitate the signal component of the
forward model for the likelihood computation, our ILI approach requires the
inclusion of noise realizations in the training data vector.

We illustrate a simplified application of LtU-ILI for NPE on GW signals,
fixing parameters related to the time and geometry of the detection and
focusing on a set of four parameters: chirp mass ($\mathcal{M}$), mass ratio
($q$), effective aligned spin ($\chi_{\textrm{eff}}$) and luminosity distance
($d_{L}$). These parameters, chosen for their reduced correlation and
significance in GW parameter estimation (_e.g.,_ Veitch  _et al._ , 2015;
Vitale  _et al._ , 2017), are derived from the masses of the coalescing
objects – with
$\mathcal{M}\equiv(m_{1}m_{2})^{\frac{3}{5}}/(m_{1}+m_{2})^{\frac{1}{5}}$ and
$q\equiv m_{2}/m_{1}\leq 1$ – and spin characteristics (Racine, 2008).

We established baseline values for our simulated event $\bm{x_{0}}$ at
{$d_{L}=390\ \textrm{Mpc},q=0.79,\mathcal{M}=30.2\
M_{\odot},\chi_{\textrm{eff}}=-0.09$)}, akin to the GW150914 event as detailed
in the first column of Table 1 by Abbott  _et al._ (2016). Following the
procedure described in Karathanasis  _et al._ (2023), the defined priors are
uniform ranging from [$100.0\textrm{ Mpc},0.2,12.0\,M_{\odot},-1.0$] to
[$1000.0\textrm{ Mpc},1.0,45.0\,M_{\odot},1.0$], with narrower margins for
chirp mass and luminosity distance to facilitate a straightforward example.
For density estimation, we utilize the sbi backend with a MAF and a 1D
convolutional embedding network, undergoing five rounds of training with 1000
simulations per round.

Figure 8b illustrates the progression of the cumulative posterior probability,
as defined in Equation 4, as a function of training time. The noticeable
increase in probability mass for both training and validation datasets is
attributed to the model’s focus on learning the posterior density near the
target observation $\mathcal{P}(\bm{\theta}|\bm{x}=\bm{x}_{0})$. However, an
early increase in validation loss during subsequent rounds suggests the
potential for overfitting, indicating the model’s limitations in gleaning
further insights from the simulation inputs. Sampling the posterior at
$\bm{x}=\bm{x_{0}}$ depicted in Figure 8c, reveals that we are able to
constrain each of the four parameters around the fiducial point.

###  5.5 Inferring Dust Parameters

Dust makes up approximately 1% of the interstellar medium, but reprocesses
around 30% of stellar emission (Silva  _et al._ , 1998). Understanding the
properties of dust and its impact on observables is, therefore, a key goal of
galactic and extragalactic astronomy, and a fundamental uncertainty in
downstream cosmological inference analyses. Dust leads to an overall reddening
of the emission, parametrized by the optical depth and attenuation law; these
are analogous to the normalization and slope, particularly where the latter is
parametrized as a power law. Below we show how these parameters can be
inferred from observational properties using a simple dust model (Trayford
_et al._ , 2015) applied to the Simba cosmological hydrodynamic simulation
(Davé  _et al._ , 2019) within an SBI framework.

We first briefly describe our forward model for the stellar emission and dust
attenuation. We generated spectra for all galaxies with stellar mass
$>10^{10}\;\mathrm{M_{\odot}}$ in the Simba simulation using Synthesizer
(Lovell et al. in prep.)333https://flaresimulations.github.io/synthesizer. The
integrated intrinsic stellar emission of each galaxy was obtained by linking
each stellar particle to the BC03 stellar population synthesis models (Bruzual
and Charlot, 2003) based on their age and metallicity. We then modelled the
dust attenuation as a wavelength-dependent power law, with slope $\alpha$,

$\displaystyle
T(\lambda,t)=\mathrm{exp}\left[\tau(t)\left(\frac{\lambda}{\lambda_{\mathrm{V}}}\right)^{-\alpha}\right]\;\;,$
(15)

where $\lambda_{\mathrm{V}}$ is the wavelength of the V band
($550\;\mathrm{nm}$). Here the optical depth $\tau$ is dependent on the age of
the stellar population $t$; star particles younger than 10 Myr are assumed to
still reside within their birth clouds (BC), and therefore experience an extra
source of attenuation,

$\displaystyle\tau(t)=\begin{cases}\tau_{\mathrm{BC}}+\tau_{\mathrm{ISM}}&;\quad
t\leqslant 10{\rm Myr}\\\ \tau_{\mathrm{ISM}}&;\quad t>10{\rm
Myr}\end{cases}\;\;.$ (16)

Assuming the same underlying stellar emission model, we generated rest frame
optical photometry in the $g$ and $r$ bands for all galaxies in Simba for 1000
simulations on a Latin hypercube, with flat priors on the parameters in the
range $\alpha\in[0.5,2.0]$, $\tau_{\mathrm{ISM}}\in[0.01,0.5]$ and
$\tau_{\mathrm{BC}}\in[0.3,1.5]$. We then take these simulations and measure
the $r$-band luminosity function and the $g-r$ color distribution, and use
these as our summary statistics. We perform NPE using an MAF architecture;
Figure 9 shows the posteriors on a test parameter – data pair. We show results
from training on just the luminosity function, just the color, and the two
combined. It is clear that the overall normalisation of the luminosity
function constrains the ISM optical depth, $\tau_{\mathrm{ISM}}$, whereas the
colour constrains the slope of the attenuation law, $\alpha$. When combined,
these parameters are simultaneously tightly constrained. The birth cloud
attenuation has a more subtle effect on the luminosity function and colour
distribution, but the combination leads to tighter constraints.

Figure 9: Dust parameter inference example on the Simba simulation. Left:
example posteriors on the ISM and birth cloud optical depth as well as the
slope of the attenuation law, from the distribution of $g-r$ colors (blue),
the $r$-band magnitude luminosity function (orange), and the combination of
both (green). Right: posterior constraints on each parameter for a range of
test simulations.

###  5.6 The Mass and Energy Loading of Galactic Winds

A major open question in galaxy formation is the role of feedback in
regulating star formation and shaping the stellar-to-halo-mass relation.
Recently, Carr  _et al._ (2023) and Pandya  _et al._ (2023) showed that
supernova-driven galactic winds can heat and stir turbulence in the
circumgalactic medium of dwarf galaxies and Milky Way-mass halos. This can
suppress the cooling and accretion of gas onto galaxies thereby helping to
regulate star formation. The strength of this “preventative feedback” is
controlled by two wind parameters: the mass loading factor ($\eta_{\rm M}$)
and the energy loading factor ($\eta_{\rm E}$). Together, these describe the
specific energy of winds and summarize the fraction of mass and energy
produced by supernovae that leave the galaxy and become available for large-
scale heating.

Here we want to explore the connection between $(\eta_{\rm M},\eta_{\rm E})$
and the stellar-to-halo-mass relation predicted by numerous realizations of
the next-generation semi-analytic model sapphire (Pandya et al in prep.). The
model used here is more similar to Carr  _et al._ (2023) than Pandya  _et al._
(2023) in that it neglects turbulence and assumes purely thermal heating of
the circumgalactic medium (CGM) for simplicity. We assume that $\eta_{\rm M}$
and $\eta_{\rm E}$ both follow power law scalings with halo virial velocity
$V_{\rm vir}$ (a proxy for halo mass):

$\eta_{\rm M}=A_{\rm M}\left(\frac{V_{\rm
vir}}{125\rm{km\;s^{-1}}}\right)^{\alpha_{\rm M}}$ (17)

$\eta_{\rm E}=A_{\rm E}\left(\frac{V_{\rm
vir}}{125\rm{km\;s^{-1}}}\right)^{\alpha_{\rm E}}$ (18)

Here, $A_{\rm M}$, $\alpha_{\rm M}$, $A_{\rm E}$ and $\alpha_{\rm E}$ are
hyperparameters that govern these power law scalings, and we have ignored any
redshift dependence. We generate 5,000 model realizations on a Latin hypercube
by uniformly sampling combinations of $A_{\rm M}\in[0.01,100]$, $\alpha_{\rm
M}\in[-2,2]$, $A_{\rm E}\in[0.01,1]$ and $\alpha_{\rm E}\in[-2,2]$. All other
parameters of the model were fixed to reasonable values following Pandya  _et
al._ (2023). With these parameters, the system of ordinary differential
equations describing the model (Equations 1, 2, and 4-8 of Pandya  _et al._
2023) is evolved using mass accretion histories of 120 halos randomly selected
from five subvolumes of the TNG100 simulation (Pillepich  _et al._ , 2018). We
find that this is a sufficient number of halos to capture variations in the
stellar-to-halo-mass relation with model parameters. We restrict to halos with
$\log M_{\rm vir}/M_{\odot}=10-12.4$ at $z=0$, spanning the dwarf and Milky
Way regime in which supernova feedback is thought to dominate.

Figure 10: Sets of galaxies with properties $\textbf{d}=\\{(M_{\rm
vir},M_{*})\\}_{i=1}^{n_{d}=120}$ are sensitive to sapphire mass and energy
loading parameters $\theta=(\eta_{M},\eta_{E})$. Here we utilize the fishnets
architecture to aggregate sets of galaxy properties before performing NPE.
Contours are shown at the central $68\%$ and $95\%$ confidence intervals for a
single data realization (left). Posterior predictive plots show our NPE
network is sensitive to both energy loading hyper-parameters but cannot learn
the mass loading hyper-parameters with the $M_{*}-M_{\rm vir}$ relation alone.

We use the fishnets architecture by Makinen  _et al._ (2023) to extract what
information the stellar-to-halo-mass relation contains about
$\theta=(\eta_{\rm M},\eta_{\rm E})$. The output from each sapphire
realization is a set of 120 independently-evolved halos $\textbf{d}=\\{(M_{\rm
vir},M_{*},\frac{M_{\rm vir}}{M_{*}})\\}_{i=1}^{n_{d}=120}$ where $M_{\rm
vir}$ and $M_{*}$ are the $z=0$ halo mass and stellar mass, respectively. The
likelihood here is a product of each individual halos’ likelihoods,
$\mathcal{L}(\\{d_{i}\\}|\theta)=\prod^{n_{d}=120}_{i=0}\mathcal{L}(d_{i}|\theta)$,
requiring an aggregation over heterogeneous data distributions. Following the
optimal fishnets formalism presented in Makinen  _et al._ (2023), we embed
data into neural score embeddings and Fisher weights, each parameterized by
fully-connected networks of size $[128,128,128]$ with LeakyRelu activations
before passing the weighted scores to two NPE networks.

We show the resulting posterior for one test simulation in Figure 10. We find
that galaxy energy loading factors are extremely sensitive to halo and stellar
mass properties (as expected), while mass loading factors prove more difficult
to constrain. This is consistent with Carr  _et al._ (2023) who showed that
the stellar-to-halo-mass relation is largely insensitive to variations in
$\eta_{\rm M}$ but very responsive to changes in $\eta_{\rm E}$. The physical
reason for this is that increasing $\eta_{\rm M}$ leads to a higher CGM
density and enhanced cooling of gas back into the galaxy without suppressing
star formation, thus leading to a similar $M_{*}-M_{\rm vir}$ relation.
Instead, increasing $\eta_{\rm E}$ can heat the CGM, prevent gas accretion,
reduce star formation and alter the normalization and shape of the
$M_{*}-M_{\rm vir}$ relation. In the future, it will be compelling to use
fishnets to understand what additional galaxy properties output by sapphire
contain the information needed to better learn both $\eta_{\rm M}$ and
$\eta_{\rm E}$.

##  6 Discussion

ILI methods are guaranteed to recover the true posterior under the assumptions
that:

  1. 1.

There is a sufficient amount of training data to cover data and parameter
space.

  2. 2.

The training data or simulator is reflective of reality.

  3. 3.

The chosen neural architecture is sufficiently flexible to capture the data-
parameter relationship.

Despite the best practices implemented in LtU-ILI, biases may occur when these
assumptions are violated. Before applying these approaches to real
observational data, it is crucial to understand these failure modes, as well
as how to build inference pipelines to be robust to them.

When training data is limited, ILI methods may overfit to a subset of
parameter space and suffer from poor out-of-sample performance. This behavior
has led to recent discussion as to whether ILI methods are systematically
overconfident (Hermans  _et al._ , 2022). PIT tests (§2.4) on independent
validation sets are a useful diagnostic tool for detecting overconfidence. For
a fixed simulation budget, this problem can be mitigated through the use of
ensembling (Hermans  _et al._ , 2022). Ensembling averages the predictions of
multiple, independently trained models together, naturally inflating
predictive uncertainty to account for epistemic bias (Lakshminarayanan  _et
al._ , 2017). Another approach is the frequentist recalibration of error bars
(Masserano  _et al._ , 2022).

Model misspecification occurs when the data-generating function for creating
training data is poorly representative of physical processes in the real world
(Cannon  _et al._ , 2022). This can cause an implicit model to interpret the
wrong data-parameter connection, resulting in biased predictions. In
astrophysics and cosmology, combating model misspecification often requires
running more realistic, high-resolution simulators which can be more
computationally demanding. In turn, this reduces the number of available
simulations, damaging assumption (1). Alternative approaches include
increasing the noise model of a given simulation (_e.g.,_ Modi  _et al._ ,
2023), training NDEs with a conservative loss to inflate uncertainties and
counteract misspecification (_e.g.,_ Kelly  _et al._ , 2023; Huang  _et al._ ,
2023), or using latent functions present within intermediate layers of the
Bayesian hierarchical model as diagnostics (_e.g.,_ Leclercq, 2022). In
addition, Bayesian model comparison (Silva  _et al._ , 1998; Handley  _et al._
, 2015) is an effective tool for distinguishing between various models
conditioned on real data, but also relies on massive simulation budgets.
Recent work such as Evidence Networks (Jeffrey and Wandelt, 2024) has shown
that this process can be greatly accelerated through the application of neural
networks.

Lastly, implicit inference can be very challenging in the high-dimensional
parameter regime (_i.e.,_ when $\operatorname{dim}(\bm{\theta})\gg 100$). As
the dimensionality of $\bm{\theta}$ increases, the number of parameters
necessary to specify a highly-correlated posterior increases exponentially.
This problem is of particular relevance in astrophysics and cosmology, in
which we often would like to infer voxelized 2D or 3D fields (_e.g.,_ Jasche
and Wandelt, 2013; de Andres  _et al._ , 2024). Recent empirical evidence
suggests that Moment Networks (Jeffrey and Wandelt, 2020) and diffusion-based
generative models (Song  _et al._ , 2020; Legin  _et al._ , 2024) are
incredibly effective neural architectures for these tasks. Though LtU-ILI does
not currently include these models, we intend to implement them in future
work.

##  7 Conclusion

In conclusion, we have introduced LtU-ILI, a comprehensive framework designed
for implicit inference in scientific applications. LtU-ILI provides extensive
routines to facilitate the training of neural networks for regression
posterior inference, aimed in particular at cases in which the likelihood is
intractable. The code encapsulates best practices and state-of-the-art models,
ensuring adaptability to a broad spectrum of scientific domains. By
consolidating the training, testing, and robustness assessment of ILI models
into a unified, modular package, LtU-ILI streamlines the process of addressing
novel inference problems while providing for production-level inference and
hyperparameter searches, making it accessible to researchers both experienced
and novice.

In this manuscript, we have explained the theoretical basis of ILI, described
the current features of LtU-ILI in detail, and provided benchmarks and science
examples on a variety of astronomy and cosmology inference problems. Our
results demonstrate that LtU-ILI yields competitive constraints on both
synthetic and real science examples, indicating that these methods are ready
out-of-the-box for application to new problems.

Looking ahead, this framework can be applied to accelerate the analyses of
large surveys such as SDSS, Planck, VRO, Euclid, DESI, and JWST. We also aim
to provide exhaustive benchmarks on the variety of problems described in
Section 5, contributing to the formalization of machine learning analyses in
astronomy and beyond. We intend to integrate LtU-ILI with state-of-the-art
hyperparameter tuning (_e.g.,_ Akiba  _et al._ , 2019) to provide an
automatic, yet principled way for doing model selection. In addition, current
approaches for calibrating posterior coverage rely on tweaking hyperparameters
and retraining. An alternative approach is to implement recalibration methods
(_e.g.,_ Masserano  _et al._ , 2022) which correctly calibrate posteriors
given an independent validation set. We reserve implementation and analysis of
these improvements to future work.

###### Acknowledgements.

We thank Rosa Malandrino, Niall Jeffrey, and Justin Alsing for useful comments
and suggestions on developing the codebase and writing the manuscript. This
project was developed as part of the Simons Collaboration on “Learning the
Universe.”

## References

  * Feigelson and Babu (2012) E. D. Feigelson and G. J. Babu,  _Modern statistical methods for astronomy: with R applications_ (Cambridge University Press, 2012).
  * Dodelson and Schmidt (2020) S. Dodelson and F. Schmidt,  _Modern cosmology_ (Academic press, 2020).
  * Eadie  _et al._ (2023) G. M. Eadie, J. S. Speagle, J. Cisewski-Kehe, D. Foreman-Mackey, D. Huppenkothen, D. E. Jones, A. Springford, and H. Tak, arXiv preprint arXiv:2302.04703 (2023).
  * National Academies of Sciences, Engineering, and Medicine (2023) National Academies of Sciences, Engineering, and Medicine,  _Pathways to Discovery in Astronomy and Astrophysics for the 2020s_ (The National Academies Press, Washington, DC, 2023).
  * Carleo  _et al._ (2019) G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, Reviews of Modern Physics 91, 045002 (2019).
  * Huppenkothen  _et al._ (2023) D. Huppenkothen, M. Ntampaka, M. Ho, M. Fouesneau, B. Nord, J. Peek, M. Walmsley, J. F. Wu, C. Avestruz, T. Buck, _et al._ , arXiv preprint arXiv:2310.12528 (2023).
  * Cranmer  _et al._ (2020) K. Cranmer, J. Brehmer, and G. Louppe, Proceedings of the National Academy of Sciences 117, 30055 (2020).
  * Marin  _et al._ (2012) J.-M. Marin, P. Pudlo, C. P. Robert, and R. J. Ryder, Statistics and computing 22, 1167 (2012).
  * Christensen  _et al._ (2001) N. Christensen, R. Meyer, L. Knox, and B. Luey, Classical and Quantum Gravity 18, 2677 (2001).
  * Brooks  _et al._ (2011) S. Brooks, A. Gelman, G. Jones, and X.-L. Meng,  _Handbook of markov chain monte carlo_ (CRC press, 2011).
  * Jeffrey  _et al._ (2021) N. Jeffrey, J. Alsing, and F. Lanusse, Monthly Notices of the Royal Astronomical Society 501, 954 (2021).
  * Makinen  _et al._ (2021) T. L. Makinen, T. Charnock, J. Alsing, and B. D. Wandelt, Journal of Cosmology and Astroparticle Physics 2021, 049 (2021).
  * Makinen  _et al._ (2022a) T. L. Makinen, T. Charnock, P. Lemos, N. Porqueres, A. F. Heavens, and B. D. Wandelt, The Open Journal of Astrophysics 5 (2022a), 10.21105/astro.2207.05202.
  * de Santi  _et al._ (2023a) N. S. de Santi, H. Shao, F. Villaescusa-Navarro, L. R. Abramo, R. Teyssier, P. Villanueva-Domingo, Y. Ni, D. Anglés-Alcázar, S. Genel, E. Hernandez-Martinez, _et al._ , arXiv preprint arXiv:2302.14101 (2023a).
  * Hahn  _et al._ (2023) C. Hahn, P. Lemos, L. Parker, B. R.-S. Blancard, M. Eickenberg, S. Ho, J. Hou, E. Massara, C. Modi, A. M. Dizgah, _et al._ , arXiv preprint arXiv:2310.15246 (2023).
  * Dax  _et al._ (2021a) M. Dax, S. R. Green, J. Gair, J. H. Macke, A. Buonanno, and B. Schölkopf, Phys. Rev. Lett. 127, 241103 (2021a), arXiv:2106.12594 [gr-qc] .
  * Cheung  _et al._ (2022) D. H. Cheung, K. W. Wong, O. A. Hannuksela, T. G. Li, and S. Ho, Physical Review D 106, 083014 (2022).
  * Ho  _et al._ (2022) M. Ho, M. Ntampaka, M. M. Rau, M. Chen, A. Lansberry, F. Ruehle, and H. Trac, Nature Astronomy 6, 936 (2022).
  * de Andres  _et al._ (2022) D. de Andres, W. Cui, F. Ruppin, M. De Petris, G. Yepes, G. Gianfagna, I. Lahouli, G. Aversano, R. Dupuis, M. Jarraya, _et al._ , Nature Astronomy 6, 1325 (2022).
  * Walmsley  _et al._ (2020) M. Walmsley, L. Smith, C. Lintott, Y. Gal, S. Bamford, H. Dickinson, L. Fortson, S. Kruk, K. Masters, C. Scarlata, _et al._ , Monthly Notices of the Royal Astronomical Society 491, 1554 (2020).
  * Ghosh  _et al._ (2022) A. Ghosh, C. M. Urry, A. Rau, L. Perreault-Levasseur, M. Cranmer, K. Schawinski, D. Stark, C. Tian, R. Ofman, T. T. Ananna, _et al._ , The Astrophysical Journal 935, 138 (2022).
  * Hermans  _et al._ (2021) J. Hermans, N. Banik, C. Weniger, G. Bertone, and G. Louppe, Monthly Notices of the Royal Astronomical Society 507, 1999 (2021).
  * Alvey  _et al._ (2023) J. Alvey, M. Gerdes, and C. Weniger, arXiv preprint arXiv:2304.02032 (2023).
  * Rogers  _et al._ (2023) J. G. Rogers, C. Janó Muñoz, J. E. Owen, and T. L. Makinen, Monthly Notices of the Royal Astronomical Society 519, 6028 (2023).
  * Aubin  _et al._ (2023) M. Aubin, C. Cuesta-Lazaro, E. Tregidga, J. Viaña, C. Garraffo, I. E. Gordon, M. López-Morales, R. J. Hargreaves, V. Y. Makhnev, J. J. Drake, _et al._ , arXiv preprint arXiv:2309.09337 (2023).
  * Gal and Ghahramani (2016) Y. Gal and Z. Ghahramani, in  _international conference on machine learning_ (PMLR, 2016) pp. 1050–1059.
  * Lakshminarayanan  _et al._ (2017) B. Lakshminarayanan, A. Pritzel, and C. Blundell, Advances in neural information processing systems 30 (2017).
  * Maddox  _et al._ (2019) W. J. Maddox, P. Izmailov, T. Garipov, D. P. Vetrov, and A. G. Wilson, Advances in neural information processing systems 32 (2019).
  * Lemos  _et al._ (2023) P. Lemos, M. Cranmer, M. Abidi, C. Hahn, M. Eickenberg, E. Massara, D. Yallup, and S. Ho, Machine Learning: Science and Technology 4, 01LT01 (2023), arXiv:2207.08435 [astro-ph.CO] .
  * Jin  _et al._ (2019) H. Jin, Q. Song, and X. Hu, in  _Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining_ (2019) pp. 1946–1956.
  * White  _et al._ (2021) C. White, W. Neiswanger, and Y. Savani, in  _Proceedings of the AAAI Conference on Artificial Intelligence_ (2021).
  * Wang  _et al._ (2023) B. Wang, J. Leja, V. A. Villar, and J. S. Speagle, ApJ 952, L10 (2023), arXiv:2304.05281 [astro-ph.IM] .
  * Modi  _et al._ (2023) C. Modi, S. Pandey, M. Ho, C. Hahn, B. Blancard, and B. Wandelt, arXiv preprint arXiv:2309.15071 (2023).
  * Tejero-Cantero  _et al._ (2020) A. Tejero-Cantero, J. Boelts, M. Deistler, J.-M. Lueckmann, C. Durkan, P. J. Gonçalves, D. S. Greenberg, and J. H. Macke, Journal of Open Source Software 5, 2505 (2020).
  * Alsing  _et al._ (2019) J. Alsing, T. Charnock, S. Feeney, and B. Wandelt, Monthly Notices of the Royal Astronomical Society 488, 4440 (2019).
  * Rozet  _et al._ (2021) F. Rozet, A. Delaunoy, B. Miller, _et al._ , “LAMPE: Likelihood-free amortized posterior estimation,”  (2021).
  * Lemos  _et al._ (2023) P. Lemos, A. Coogan, Y. Hezaveh, and L. P. Levasseur, in  _International Conference on Machine Learning, ICML 2023, 23-29 July 2023, Honolulu, Hawaii, USA_, Proceedings of Machine Learning Research, Vol. 202, edited by A. Krause, E. Brunskill, K. Cho, B. Engelhardt, S. Sabato, and J. Scarlett (PMLR, 2023) pp. 19256–19273.
  * Zhang  _et al._ (2023) K. Zhang, J. Bloom, and N. Hernitschek, in  _NeurIPS 2023 Workshop on Deep Learning and Inverse Problems_ (2023).
  * Tegmark  _et al._ (2004) M. Tegmark, M. R. Blanton, M. A. Strauss, F. Hoyle, D. Schlegel, R. Scoccimarro, M. S. Vogeley, D. H. Weinberg, I. Zehavi, A. Berlind, T. Budavari, A. Connolly, D. J. Eisenstein, D. Finkbeiner, J. A. Frieman, J. E. Gunn, A. J. S. Hamilton, L. Hui, B. Jain, D. Johnston, S. Kent, H. Lin, R. Nakajima, R. C. Nichol, J. P. Ostriker, A. Pope, R. Scranton, U. Seljak, R. K. Sheth, A. Stebbins, A. S. Szalay, I. Szapudi, L. Verde, Y. Xu, J. Annis, N. A. Bahcall, J. Brinkmann, S. Burles, F. J. Castander, I. Csabai, J. Loveday, M. Doi, M. Fukugita, I. Gott, J. Richard, G. Hennessy, D. W. Hogg, Ž. Ivezić, G. R. Knapp, D. Q. Lamb, B. C. Lee, R. H. Lupton, T. A. McKay, P. Kunszt, J. A. Munn, L. O’Connell, J. Peoples, J. R. Pier, M. Richmond, C. Rockosi, D. P. Schneider, C. Stoughton, D. L. Tucker, D. E. Vanden Berk, B. Yanny, D. G. York, and SDSS Collaboration, ApJ 606, 702 (2004), arXiv:astro-ph/0310725 [astro-ph] .
  * Christensen and Meyer (2022) N. Christensen and R. Meyer, Reviews of Modern Physics 94, 025001 (2022).
  * Porqueres  _et al._ (2022) N. Porqueres, A. Heavens, D. Mortlock, and G. Lavaux, Monthly Notices of the Royal Astronomical Society 509, 3194 (2022).
  * Boese and Doebereiner (2001) F. Boese and S. Doebereiner, Astronomy & Astrophysics 370, 649 (2001).
  * Braun  _et al._ (2008) J. Braun, J. Dumm, F. De Palma, C. Finley, A. Karle, and T. Montaruli, Astroparticle Physics 29, 299 (2008).
  * Tsaprazi  _et al._ (2022) E. Tsaprazi, N.-M. Nguyen, J. Jasche, F. Schmidt, and G. Lavaux, J. Cosmology Astropart. Phys 2022, 003 (2022), arXiv:2112.04484 [astro-ph.CO] .
  * Dingeldein  _et al._ (2023) L. Dingeldein, P. Cossio, and R. Covino, Machine Learning: Science and Technology 4, 025009 (2023).
  * Ntampaka  _et al._ (2019) M. Ntampaka, J. ZuHone, D. Eisenstein, D. Nagai, A. Vikhlinin, L. Hernquist, F. Marinacci, D. Nelson, R. Pakmor, A. Pillepich, _et al._ , The Astrophysical Journal 876, 82 (2019).
  * Fluri  _et al._ (2022) J. Fluri, T. Kacprzak, A. Lucchi, A. Schneider, A. Refregier, and T. Hofmann, Phys. Rev. D 105, 083518 (2022).
  * Lanusse  _et al._ (2018) F. Lanusse, Q. Ma, N. Li, T. E. Collett, C.-L. Li, S. Ravanbakhsh, R. Mandelbaum, and B. Póczos, Monthly Notices of the Royal Astronomical Society 473, 3895 (2018).
  * Rubin (1984) D. B. Rubin, The Annals of Statistics 12, 1151 (1984).
  * Schafer and Freeman (2012) C. M. Schafer and P. E. Freeman, in  _Statistical Challenges in Modern Astronomy V_ (Springer, 2012) pp. 3–19.
  * Cameron and Pettitt (2012) E. Cameron and A. Pettitt, Monthly Notices of the Royal Astronomical Society 425, 44 (2012).
  * Hahn  _et al._ (2017) C. Hahn, M. Vakili, K. Walsh, A. P. Hearin, D. W. Hogg, and D. Campbell, Monthly Notices of the Royal Astronomical Society 469, 2791 (2017).
  * Alsing  _et al._ (2018) J. Alsing, B. Wandelt, and S. Feeney, MNRAS 477, 2874 (2018), arXiv:1801.01497 [astro-ph.CO] .
  * Papamakarios (2019) G. Papamakarios, arXiv preprint arXiv:1910.13233 (2019).
  * Bishop (1994) C. M. Bishop,  _Mixture density networks_ , WorkingPaper (Aston University, 1994).
  * Papamakarios and Murray (2016) G. Papamakarios and I. Murray, Advances in neural information processing systems 29 (2016).
  * Papamakarios  _et al._ (2021) G. Papamakarios, E. Nalisnick, D. J. Rezende, S. Mohamed, and B. Lakshminarayanan, The Journal of Machine Learning Research 22, 2617 (2021).
  * Greenberg  _et al._ (2019) D. Greenberg, M. Nonnenmacher, and J. Macke, in  _International Conference on Machine Learning_ (PMLR, 2019) pp. 2404–2414.
  * Vasist  _et al._ (2023) M. Vasist, F. Rozet, O. Absil, P. Mollière, E. Nasedkin, and G. Louppe, A&A 672, A147 (2023), arXiv:2301.06575 [astro-ph.EP] .
  * Crisostomi  _et al._ (2023) M. Crisostomi, K. Dey, E. Barausse, and R. Trotta, Phys. Rev. D 108, 044029 (2023), arXiv:2305.18528 [gr-qc] .
  * Alsing  _et al._ (2019) J. Alsing, T. Charnock, S. Feeney, and B. Wandelt, MNRAS 488, 4440 (2019), arXiv:1903.00007 [astro-ph.CO] .
  * Papamakarios  _et al._ (2019) G. Papamakarios, D. Sterratt, and I. Murray, in  _The 22nd International Conference on Artificial Intelligence and Statistics_ (PMLR, 2019) pp. 837–848.
  * Robert  _et al._ (1999) C. P. Robert, G. Casella, and G. Casella,  _Monte Carlo statistical methods_ , Vol. 2 (Springer, 1999).
  * Blei  _et al._ (2017) D. M. Blei, A. Kucukelbir, and J. D. McAuliffe, Journal of the American statistical Association 112, 859 (2017).
  * Robert  _et al._ (2004) C. P. Robert, G. Casella, C. P. Robert, and G. Casella, Monte Carlo statistical methods , 267 (2004).
  * Jeffrey  _et al._ (2021) N. Jeffrey, J. Alsing, and F. Lanusse, MNRAS 501, 954 (2021), arXiv:2009.08459 [astro-ph.CO] .
  * Brehmer  _et al._ (2020) J. Brehmer, F. Kling, I. Espejo, and K. Cranmer, Comput. Softw. Big Sci. 4, 3 (2020), arXiv:1907.10621 [hep-ph] .
  * Hermans  _et al._ (2020) J. Hermans, V. Begy, and G. Louppe, in  _International conference on machine learning_ (PMLR, 2020) pp. 4239–4248.
  * Cranmer  _et al._ (2015) K. Cranmer, J. Pavez, and G. Louppe, arXiv preprint arXiv:1506.02169 (2015).
  * Cole  _et al._ (2022) A. Cole, B. K. Miller, S. J. Witte, M. X. Cai, M. W. Grootes, F. Nattino, and C. Weniger, J. Cosmology Astropart. Phys 2022, 004 (2022), arXiv:2111.08030 [astro-ph.CO] .
  * Karchev  _et al._ (2023) K. Karchev, R. Trotta, and C. Weniger, MNRAS 520, 1056 (2023), arXiv:2209.06733 [astro-ph.CO] .
  * Miller  _et al._ (2021) B. Miller, A. Cole, P. Forré, G. Louppe, and C. Weniger, Advances in Neural Information Processing Systems 34, 129 (2021), arXiv:2107.01214 [stat.ML] .
  * Delaunoy  _et al._ (2020) A. Delaunoy, A. Wehenkel, T. Hinderer, S. Nissanke, C. Weniger, A. R. Williamson, and G. Louppe, arXiv e-prints , arXiv:2010.12931 (2020), arXiv:2010.12931 [astro-ph.IM] .
  * Bhardwaj  _et al._ (2023) U. Bhardwaj, J. Alvey, B. K. Miller, S. Nissanke, and C. Weniger, Phys. Rev. D 108, 042004 (2023), arXiv:2304.02035 [gr-qc] .
  * Song  _et al._ (2020) Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole, arXiv preprint arXiv:2011.13456 (2020).
  * Lueckmann  _et al._ (2021) J.-M. Lueckmann, J. Boelts, D. Greenberg, P. Goncalves, and J. Macke, in  _International conference on artificial intelligence and statistics_ (PMLR, 2021) pp. 343–351.
  * Watson-Parris  _et al._ (2021) D. Watson-Parris, A. Williams, L. Deaconu, and P. Stier, Geoscientific Model Development 14, 7659 (2021).
  * Gutmann  _et al._ (2016) M. U. Gutmann, J. Cor, _et al._ , Journal of Machine Learning Research 17, 1 (2016).
  * Leclercq (2018) F. Leclercq, Phys. Rev. D 98, 063511 (2018), arXiv:1805.07152 [astro-ph.CO] .
  * Lopez-Paz and Oquab (2016) D. Lopez-Paz and M. Oquab, arXiv preprint arXiv:1610.06545 (2016).
  * Cook  _et al._ (2006) S. R. Cook, A. Gelman, and D. B. Rubin, Journal of Computational and Graphical Statistics 15, 675 (2006).
  * Zhao  _et al._ (2021) D. Zhao, N. Dalmasso, R. Izbicki, and A. B. Lee, in  _Uncertainty in Artificial Intelligence_ (PMLR, 2021) pp. 1830–1840.
  * Bordoloi  _et al._ (2010) R. Bordoloi, S. J. Lilly, and A. Amara, Monthly Notices of the Royal Astronomical Society 406, 881 (2010).
  * Tanaka  _et al._ (2018) M. Tanaka, J. Coupon, B.-C. Hsieh, S. Mineo, A. J. Nishizawa, J. Speagle, H. Furusawa, S. Miyazaki, and H. Murayama, Publications of the Astronomical Society of Japan 70, S9 (2018).
  * Talts  _et al._ (2018) S. Talts, M. Betancourt, D. Simpson, A. Vehtari, and A. Gelman, arXiv e-prints , arXiv:1804.06788 (2018), arXiv:1804.06788 [stat.ME] .
  * Wang  _et al._ (2019) A. Y. Wang, A. Mittal, C. Brooks, and S. Oney, Proceedings of the ACM on Human-Computer Interaction 3, 1 (2019).
  * Singh and Singh (2020) D. Singh and B. Singh, Applied Soft Computing 97, 105524 (2020).
  * Kingma and Ba (2014) D. P. Kingma and J. Ba, arXiv preprint arXiv:1412.6980 (2014).
  * Bai  _et al._ (2021) Y. Bai, E. Yang, B. Han, Y. Yang, J. Li, Y. Mao, G. Niu, and T. Liu, Advances in Neural Information Processing Systems 34, 24392 (2021).
  * Abadi  _et al._ (2016) M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, _et al._ , in  _12th USENIX symposium on operating systems design and implementation (OSDI 16)_ (2016) pp. 265–283.
  * Paszke  _et al._ (2019) A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, _et al._ , Advances in neural information processing systems 32 (2019).
  * Charnock  _et al._ (2018) T. Charnock, G. Lavaux, and B. D. Wandelt, Physical Review D 97, 083004 (2018).
  * Alsing and Wandelt (2018) J. Alsing and B. Wandelt, Monthly Notices of the Royal Astronomical Society: Letters 476, L60 (2018).
  * Alsing  _et al._ (2018) J. Alsing, B. Wandelt, and S. Feeney, Monthly Notices of the Royal Astronomical Society 477, 2874 (2018).
  * Lueckmann  _et al._ (2017) J.-M. Lueckmann, P. J. Goncalves, G. Bassetto, K. Öcal, M. Nonnenmacher, and J. H. Macke, Advances in neural information processing systems 30 (2017).
  * Durkan  _et al._ (2020) C. Durkan, A. Bekasov, I. Murray, and G. Papamakarios, “nflows: normalizing flows in PyTorch,”  (2020).
  * Rozet  _et al._ (2022) F. Rozet _et al._ ,  _Zuko: Normalizing flows in PyTorch_ (2022).
  * Hermans  _et al._ (2022) J. Hermans, A. Delaunoy, F. Rozet, A. Wehenkel, V. Begy, and G. Louppe, stat. 1050 (2022).
  * Blundell  _et al._ (2015) C. Blundell, J. Cornebise, K. Kavukcuoglu, and D. Wierstra, in  _International conference on machine learning_ (PMLR, 2015) pp. 1613–1622.
  * Cobb and Jalaian (2021) A. D. Cobb and B. Jalaian, in  _Uncertainty in Artificial Intelligence_ (PMLR, 2021) pp. 675–685.
  * Chollet  _et al._ (2015) F. Chollet _et al._ , “Keras,” https://github.com/fchollet/keras (2015).
  * LeCun  _et al._ (1998) Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner, Proceedings of the IEEE 86, 2278 (1998).
  * Kipf and Welling (2016) T. N. Kipf and M. Welling, arXiv preprint arXiv:1609.02907 (2016).
  * Foreman-Mackey  _et al._ (2013) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, Publications of the Astronomical Society of the Pacific 125, 306 (2013).
  * Goodman and Weare (2010) J. Goodman and J. Weare, Communications in applied mathematics and computational science 5, 65 (2010).
  * Bingham  _et al._ (2019) E. Bingham, J. P. Chen, M. Jankowiak, F. Obermeyer, N. Pradhan, T. Karaletsos, R. Singh, P. A. Szerlip, P. Horsfall, and N. D. Goodman, J. Mach. Learn. Res. 20, 28:1 (2019).
  * Neal (2003) R. M. Neal, The annals of statistics 31, 705 (2003).
  * Neal  _et al._ (2011) R. M. Neal _et al._ , Handbook of markov chain monte carlo 2, 2 (2011).
  * Hoffman  _et al._ (2014) M. D. Hoffman, A. Gelman, _et al._ , J. Mach. Learn. Res. 15, 1593 (2014).
  * Graves (2011) A. Graves, Advances in neural information processing systems 24 (2011).
  * Bartlett  _et al._ (2023) D. J. Bartlett, L. Kammerer, G. Kronberger, H. Desmond, P. G. Ferreira, B. D. Wandelt, B. Burlacu, D. Alonso, and M. Zennaro, arXiv preprint arXiv:2311.15865 (2023).
  * Greenberg  _et al._ (2019) D. S. Greenberg, M. Nonnenmacher, and J. H. Macke, arXiv e-prints , arXiv:1905.07488 (2019), arXiv:1905.07488 [cs.LG] .
  * Ramesh  _et al._ (2022) P. Ramesh, J.-M. Lueckmann, J. Boelts, Á. Tejero-Cantero, D. S. Greenberg, P. J. Gonçalves, and J. H. Macke, arXiv preprint arXiv:2203.06481 (2022).
  * Ho  _et al._ (2023) M. Ho, J. Soltis, A. Farahi, D. Nagai, A. Evrard, and M. Ntampaka, arXiv preprint arXiv:2303.00005 (2023).
  * Dolag  _et al._ (2016) K. Dolag, E. Komatsu, and R. Sunyaev, Monthly Notices of the Royal Astronomical Society 463, 1797 (2016).
  * Soltis  _et al._ (2022) J. Soltis, M. Ntampaka, J. F. Wu, J. ZuHone, A. Evrard, A. Farahi, M. Ho, and D. Nagai, The Astrophysical Journal 940, 60 (2022).
  * Papamakarios  _et al._ (2017) G. Papamakarios, T. Pavlakou, and I. Murray, Advances in neural information processing systems 30 (2017).
  * Villaescusa-Navarro  _et al._ (2020) F. Villaescusa-Navarro, C. Hahn, E. Massara, A. Banerjee, A. M. Delgado, D. K. Ramanah, T. Charnock, E. Giusarma, Y. Li, E. Allys, A. Brochard, C. Uhlemann, C.-T. Chiang, S. He, A. Pisani, A. Obuljen, Y. Feng, E. Castorina, G. Contardo, C. D. Kreisch, A. Nicola, J. Alsing, R. Scoccimarro, L. Verde, M. Viel, S. Ho, S. Mallat, B. Wandelt, and D. N. Spergel, ApJS 250, 2 (2020), arXiv:1909.05273 [astro-ph.CO] .
  * Makinen  _et al._ (2022b) T. L. Makinen, T. Charnock, P. Lemos, N. Porqueres, A. Heavens, and B. D. Wandelt, The Open Journal of Astrophysics 5 (2022b), 10.21105/astro.2207.05202.
  * Cuesta-Lazaro and Mishra-Sharma (2023) C. Cuesta-Lazaro and S. Mishra-Sharma, arXiv preprint arXiv:2311.17141 (2023).
  * de Santi  _et al._ (2023b) N. S. M. de Santi, H. Shao, F. Villaescusa-Navarro, L. R. Abramo, R. Teyssier, P. Villanueva-Domingo, Y. Ni, D. Anglés-Alcázar, S. Genel, E. Hernández-Martínez, U. P. Steinwandel, C. C. Lovell, K. Dolag, T. Castro, and M. Vogelsberger, The Astrophysical Journal 952, 69 (2023b).
  * Cutler and Flanagan (1994) C. Cutler and É. E. Flanagan, Phys. Rev. D 49, 2658 (1994), arXiv:gr-qc/9402014 [gr-qc] .
  * Thrane and Talbot (2019) E. Thrane and C. Talbot, PASA 36, e010 (2019), arXiv:1809.02293 [astro-ph.IM] .
  * Green  _et al._ (2020) S. R. Green, C. Simpson, and J. Gair, Phys. Rev. D 102, 104057 (2020), arXiv:2002.07656 [astro-ph.IM] .
  * Green and Gair (2020) S. R. Green and J. Gair, arXiv e-prints , arXiv:2008.03312 (2020), arXiv:2008.03312 [astro-ph.IM] .
  * Dax  _et al._ (2021b) M. Dax, S. R. Green, J. Gair, M. Deistler, B. Schölkopf, and J. H. Macke, arXiv e-prints , arXiv:2111.13139 (2021b), arXiv:2111.13139 [cs.LG] .
  * Wildberger  _et al._ (2023) J. B. Wildberger, M. Dax, S. Buchholz, S. R. Green, J. Macke, and B. Schölkopf, in  _Machine Learning for Astrophysics_ (2023) p. 34, arXiv:2305.17161 [cs.LG] .
  * Khan  _et al._ (2016) S. Khan, S. Husa, M. Hannam, F. Ohme, M. Pürrer, X. J. Forteza, and A. Bohé, Phys. Rev. D 93, 044007 (2016), arXiv:1508.07253 [gr-qc] .
  * Ashton  _et al._ (2019) G. Ashton, M. Hübner, P. D. Lasky, C. Talbot, K. Ackley, S. Biscoveanu, Q. Chu, A. Divakarla, P. J. Easter, B. Goncharov, F. Hernandez Vivanco, J. Harms, M. E. Lower, G. D. Meadors, D. Melchor, E. Payne, M. D. Pitkin, J. Powell, N. Sarin, R. J. E. Smith, and E. Thrane, ApJS 241, 27 (2019), arXiv:1811.02042 [astro-ph.IM] .
  * Romero-Shaw  _et al._ (2020) I. M. Romero-Shaw, C. Talbot, S. Biscoveanu, V. D’Emilio, G. Ashton, C. P. L. Berry, S. Coughlin, S. Galaudage, C. Hoy, M. Hübner, K. S. Phukon, M. Pitkin, M. Rizzo, N. Sarin, R. Smith, S. Stevenson, A. Vajpeyi, M. Arène, K. Athar, S. Banagiri, N. Bose, M. Carney, K. Chatziioannou, J. A. Clark, M. Colleoni, R. Cotesta, B. Edelman, H. Estellés, C. García-Quirós, A. Ghosh, R. Green, C. J. Haster, S. Husa, D. Keitel, A. X. Kim, F. Hernandez-Vivanco, I. Magaña Hernandez, C. Karathanasis, P. D. Lasky, N. De Lillo, M. E. Lower, D. Macleod, M. Mateu-Lucena, A. Miller, M. Millhouse, S. Morisaki, S. H. Oh, S. Ossokine, E. Payne, J. Powell, G. Pratten, M. Pürrer, A. Ramos-Buades, V. Raymond, E. Thrane, J. Veitch, D. Williams, M. J. Williams, and L. Xiao, MNRAS 499, 3295 (2020), arXiv:2006.00714 [astro-ph.IM] .
  * Veitch  _et al._ (2015) J. Veitch, V. Raymond, B. Farr, W. Farr, P. Graff, S. Vitale, B. Aylott, K. Blackburn, N. Christensen, M. Coughlin, W. Del Pozzo, F. Feroz, J. Gair, C. J. Haster, V. Kalogera, T. Littenberg, I. Mandel, R. O’Shaughnessy, M. Pitkin, C. Rodriguez, C. Röver, T. Sidery, R. Smith, M. Van Der Sluys, A. Vecchio, W. Vousden, and L. Wade, Phys. Rev. D 91, 042003 (2015), arXiv:1409.7215 [gr-qc] .
  * Vitale  _et al._ (2017) S. Vitale, R. Lynch, V. Raymond, R. Sturani, J. Veitch, and P. Graff, Phys. Rev. D 95, 064053 (2017), arXiv:1611.01122 [gr-qc] .
  * Racine (2008) É. Racine, Phys. Rev. D 78, 044021 (2008), arXiv:0803.1820 [gr-qc] .
  * Abbott  _et al._ (2016) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 241102 (2016), arXiv:1602.03840 [gr-qc] .
  * Karathanasis  _et al._ (2023) C. Karathanasis, B. Revenu, S. Mukherjee, and F. Stachurski, Astron. Astrophys. 677, A124 (2023), arXiv:2210.05724 [astro-ph.CO] .
  * Silva  _et al._ (1998) L. Silva, G. L. Granato, A. Bressan, and L. Danese, ApJ 509, 103 (1998).
  * Trayford  _et al._ (2015) J. W. Trayford, T. Theuns, R. G. Bower, J. Schaye, M. Furlong, M. Schaller, C. S. Frenk, R. A. Crain, C. D. Vecchia, and I. G. McCarthy, Monthly Notices of the Royal Astronomical Society 452, 2879 (2015), https://academic.oup.com/mnras/article-pdf/452/3/2879/4922335/stv1461.pdf .
  * Davé  _et al._ (2019) R. Davé, D. Anglés-Alcázar, D. Narayanan, Q. Li, M. H. Rafieferantsoa, and S. Appleby, MNRAS 486, 2827 (2019), arXiv:1901.10203 [astro-ph.GA] .
  * Bruzual and Charlot (2003) G. Bruzual and S. Charlot, MNRAS 344, 1000 (2003), arXiv:astro-ph/0309134 [astro-ph] .
  * Carr  _et al._ (2023) C. Carr, G. L. Bryan, D. B. Fielding, V. Pandya, and R. S. Somerville, ApJ 949, 21 (2023), arXiv:2211.05115 [astro-ph.GA] .
  * Pandya  _et al._ (2023) V. Pandya, D. B. Fielding, G. L. Bryan, C. Carr, R. S. Somerville, J. Stern, C.-A. Faucher-Giguère, Z. Hafen, D. Anglés-Alcázar, and J. C. Forbes, ApJ 956, 118 (2023), arXiv:2211.09755 [astro-ph.GA] .
  * Pillepich  _et al._ (2018) A. Pillepich, V. Springel, D. Nelson, S. Genel, J. Naiman, R. Pakmor, L. Hernquist, P. Torrey, M. Vogelsberger, R. Weinberger, and F. Marinacci, MNRAS 473, 4077 (2018), arXiv:1703.02970 [astro-ph.GA] .
  * Makinen  _et al._ (2023) T. L. Makinen, J. Alsing, and B. D. Wandelt, “Fishnets: Information-optimal, scalable aggregation for sets and graphs,” (2023), arXiv:2310.03812 [cs.LG] .
  * Masserano  _et al._ (2022) L. Masserano, T. Dorigo, R. Izbicki, M. Kuusela, and A. B. Lee, arXiv preprint arXiv:2205.15680 (2022).
  * Cannon  _et al._ (2022) P. Cannon, D. Ward, and S. M. Schmon, arXiv preprint arXiv:2209.01845 (2022).
  * Kelly  _et al._ (2023) R. P. Kelly, D. J. Nott, D. T. Frazier, D. J. Warne, and C. Drovandi, arXiv preprint arXiv:2301.13368 (2023).
  * Huang  _et al._ (2023) D. Huang, A. Bharti, A. Souza, L. Acerbi, and S. Kaski, arXiv preprint arXiv:2305.15871 (2023).
  * Leclercq (2022) F. Leclercq, arXiv e-prints , arXiv:2209.11057 (2022), arXiv:2209.11057 [stat.ME] .
  * Handley  _et al._ (2015) W. Handley, M. Hobson, and A. Lasenby, Monthly Notices of the Royal Astronomical Society: Letters 450, L61 (2015).
  * Jeffrey and Wandelt (2024) N. Jeffrey and B. D. Wandelt, Machine Learning: Science and Technology 5, 015008 (2024).
  * Jasche and Wandelt (2013) J. Jasche and B. D. Wandelt, Monthly Notices of the Royal Astronomical Society 432, 894 (2013).
  * de Andres  _et al._ (2024) D. de Andres, W. Cui, G. Yepes, M. De Petris, A. Ferragamo, F. De Luca, G. Aversano, and D. Rennehan, Monthly Notices of the Royal Astronomical Society , stae071 (2024).
  * Jeffrey and Wandelt (2020) N. Jeffrey and B. D. Wandelt, arXiv preprint arXiv:2011.05991 (2020).
  * Legin  _et al._ (2024) R. Legin, M. Ho, P. Lemos, L. Perreault-Levasseur, S. Ho, Y. Hezaveh, and B. Wandelt, Monthly Notices of the Royal Astronomical Society: Letters 527, L173 (2024).
  * Akiba  _et al._ (2019) T. Akiba, S. Sano, T. Yanase, T. Ohta, and M. Koyama, in  _Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining_ (2019) pp. 2623–2631.