<table><caption>**Algorithm 2.4:** Multiple-Try Metropolis (MTM) algorithm</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>for</td><td>t = 1,..., T<sub>MC</sub> do</td></tr><tr><td>2</td><td></td><td>// Propose K candidates, select one</td></tr><tr><td>3</td><td></td><td>Θ<sub>c</sub><sup>(k)</sup> ~ q(Θ|Θ<sup>(t-1)</sup>) for k = 1,..., K</td></tr><tr><td>4</td><td></td><td>w(Θ<sub>c</sub><sup>(k)</sup>) for k = 1,..., K</td></tr><tr><td>5</td><td></td><td>w<sub>k</sub> for k = 1,..., K</td></tr><tr><td>6</td><td></td><td>i</td></tr><tr><td>7</td><td></td><td>// select one candidate following probabilities w<sub>k</sub></td></tr><tr><td>8</td><td></td><td>// Accept or reject</td></tr><tr><td>6</td><td></td><td>ρ&#771;<sup>(t)</sup></td></tr><tr><td>7</td><td></td><td>Draw ζ ~ Unif(0,1)</td></tr><tr><td>8</td><td></td><td>Θ<sup>(t)</sup> = Θ<sub>c</sub><sup>(i)</sup> if ζ ≤ ρ&#771;<sup>(t)</sup> else Θ<sup>(t-1)</sup></td></tr><tr><td colspan="3">Output: Chain of samples {Θ<sup>(t)</sup>}<sub>t=T<sub>BI</sub>+1</sub><sup>T<sub>MC</sub></sup></td></tr></tbody></table>

**Advantages of MTM compared to MH** – Using the MTM algorithm instead of MH might not seem advantageous at first, as it requires more computations per iteration. It may also seem that one step of MTM would be equivalent to *K* steps of MH, but with loss of the intermediate iterates. However, MTM is roughly to MH what chromatic Gibbs is to Gibbs sampling. The generation of the *K* candidates and the evaluation of their posterior pdf can be vectorized or performed in parallel. One step of MTM is therefore much faster than *K* iterations of MH, and about as fast as 1.

Besides, resorting to an MTM sampling kernel shortens the burn-in phase. In Chapter 5 (Section 5.2), we propose an MCMC algorithm that combines two kernels. Having at least one of the kernels based on the MTM algorithm permits to reach the posterior high probability region with few iterations. The other proposed kernel is based on RMSProp preconditioning (Eq. 2.21), which keeps in memory the history of the log-posterior gradient. Reaching the high probability region faster enables this second kernel to focus on gradient properties within the high probability region.

## 2.3 Comparing and checking observation models

As stated in Gelman et al. (2015, chapter 6): “Once we have accomplished the first two steps of a Bayesian analysis — constructing a probability model and computing the posterior distribution of all [physical parameters] — we should not ignore the relatively easy step of assessing the fit of the model to the data and to our substantive knowledge.”

Solving an inverse problem relies on the assumption that “the forward model accurately simulates the observed phenomena and that the uncertainty model accurately describes the uncertainty sources that affect observations”. For the inverse problem considered in this thesis, one might question the relevance of the Meudon PDR code to model a given region, the validity of the values of the secondary parameters listed in Table 1.3, or the chosen uncertainty model. These choices for the observation model have a great impact on the estimates. Unrealistic choices can greatly impact the estimation relevance for an end user. One thus needs to check whether the aforementioned hypothesis is satisfied or not, or compare models to select the one that is most compatible with the observations.

In this section, after briefly mentioning model selection and evaluation, we describe a well established Bayesian model assessment approach. Applications of both model evaluation and assessment in the ISM are described in Chapter 3 (Section 3.3). In this thesis, as considering more than one model is expensive, we focus on model assessment. The associated Bayesian hypothesis testing method is extended to be more robust in Chapter 5 (Section 5.3).
### 2.3.1 Model selection and evaluation

As already stated in Section 2.1, we group the forward model f and the noise model A in a single observation model M. Model selection consists in comparing the results obtained with different models M<sub>i</sub> and choosing the best one with respect to a quantitative criterion. There are two dominant approaches in Bayesian statistics. The first relies on the Bayesian evidence π(Y|M<sub>i</sub>) (Eq. 2.1), the second on the expected log-predictive density (elpd) and information criteria. We succinctly describe the two approaches.

**Bayesian evidence approach** – Having access to the evidence, also called marginal likelihood, for each considered model M<sub>i</sub> permits to perform Bayesian model selection. Indeed, in the case of two competing models M<sub>1</sub> and M<sub>2</sub>, the Bayes theorem yields the relative model posterior probability

$$ \frac{\pi(M_1|\mathbf{Y})}{\pi(M_2|\mathbf{Y})} = \frac{\pi(\mathbf{Y}|M_1) \pi(M_1)}{\pi(\mathbf{Y}|M_2) \pi(M_2)}, \qquad (2.45) $$

where π(M<sub>1</sub>) and π(M<sub>2</sub>) encode prior model preferences and sum to 1.

Computing the evidence is generally hard. Some methods exploit samples from the posterior – see e.g., McEwen et al. (2022), or Chib and Jeliazkov (2001) for MH samplers. Nested sampling (Skilling, 2004; Skilling, 2006) focuses on computing the evidence. Sequential MC (Del Moral et al., 2006), computes it as a byproduct. Nested sampling and SMC are described in Appendix 2.A. For a review on evidence estimation, see e.g., Friel and Wyse (2012).

**Information criteria approach** – A second option is to measure its so-called predictive accuracy, which relies on predicted observations ŶY. The predicted observations, also called reproduced observations, are distributed according to the posterior predictive distribution, given by

$$ \pi(\tilde{\mathbf{Y}}|\mathbf{Y}, \mathcal{M}) = \int \pi(\tilde{\mathbf{Y}}|\mathbf{\Theta}, \mathcal{M}) \pi(\mathbf{\Theta}|\mathbf{Y}, \mathcal{M}) d\mathbf{\Theta}. \qquad (2.46) $$

In the context of atomic or molecular emission line observations, these predicted observations ŶY can correspond to lines not used during the inversion. They can also simply correspond to synthetic reproduced observations of the same environment, i.e., with the same underlying physical conditions Θ, but with different noise realizations. In the following, we consider the latter case.

The predictive accuracy is usually measured with the expected log-predictive density (elpd) (Gelman et al., 2015, chapter 7),

$$ \text{elpd} = \sum_{n=1}^{N} \sum_{\ell=1}^{L} \int \ln \pi(\tilde{y}_{n\ell} | \mathbf{Y}, \mathcal{M}) \pi_t(\tilde{y}_{n\ell}) d\tilde{y}_{n\ell} \qquad (2.47) $$

$$ = \sum_{n=1}^{N} \sum_{\ell=1}^{L} \int \ln \left[ \int \pi(\tilde{y}_{n\ell} | \mathbf{\Theta}, \mathcal{M}) \pi(\mathbf{\Theta} | \mathbf{Y}, \mathcal{M}) d\mathbf{\Theta} \right] \pi_t(\tilde{y}_{n\ell}) d\tilde{y}_{n\ell}, \qquad (2.48) $$

with π<sub>t</sub>(ỹ<sub>nℓ</sub>) the pdf of the true generating process of new observations ŷ<sub>nℓ</sub>, which is unknown. The elpd therefore cannot be evaluated directly, but can be approximated with cross-validation or information criteria. Cross-validation is computationally demanding, as it requires solving multiple inverse problems instead of one, but is accurate. Information criteria do not require solving additional inverse problems and are quite cheap to evaluate. They rely on a simplified estimator of the elpd, the log-predictive density (lpd)

$$ \text{lpd} = \sum_{n=1}^{N} \sum_{\ell=1}^{L} \ln \pi(y_{n\ell} | \mathbf{Y}, \mathcal{M}) \qquad (2.49) $$

$$ = \sum_{n=1}^{N} \sum_{\ell=1}^{L} \ln \left[ \int \pi(y_{n\ell} | \mathbf{\Theta}, \mathcal{M}) \pi(\mathbf{\Theta} | \mathbf{Y}, \mathcal{M}) d\mathbf{\Theta} \right], \qquad (2.50) $$
which exploits observations **Y** used in the inversion to define the posterior distribution. Using $T_{MC}$ samples $\Theta^{(t)}$ of the posterior distribution, e.g., obtained with an MCMC algorithm, the lpd can be estimated with

$$ \widehat{\text{lpd}} = \sum_{n=1}^{N} \sum_{\ell=1}^{L} \ln \left[ \frac{1}{T_{\text{MC}}} \sum_{t=1}^{T_{\text{MC}}} \pi (y_{n\ell} | \Theta^{(t)}, \mathcal{M}) \right]. \quad (2.51) $$

The lpd uses the observations **Y** in both the posterior distribution $\pi(\Theta|\mathbf{Y})$ and in the predictive distribution $\pi(y_{n\ell}|\Theta^{(t)}, \mathcal{M})$. Therefore, the lpd introduces an optimistic bias in the accuracy measure. In particular, the lpd might be largely overestimated in case of overfitting.

In the information criterion approach, this optimistic bias in the lpd is compensated with a positive correction term that is subtracted from the lpd. For any information criterion IC, the elpd thus is estimated with

$$ \widehat{\text{elpd}}_{\text{IC}} = \widehat{\text{lpd}} - \hat{p}_{\text{IC}} \quad (2.52) $$

where $\hat{p}_{\text{IC}}$ is called an *effective number of parameters*. Its definition depends on the choice of the information criterion. For instance, the deviance information criterion (DIC) (Spiegelhalter et al., 2002), the widely applicable information criterion (WAIC) (Watanabe, 2010), or the WBIC (Watanabe, 2012) rely on different definitions of this effective number of parameters. We now present two examples of common information criteria and the associated effective number of parameters.

**Example 1: Akaike information criterion (AIC)** – Although we favor probabilistic approaches in this thesis, the information criterion approach can be used with a point estimate $\hat{\Theta}$, such as the maximum likelihood estimator or the maximum a posteriori. The posterior $\pi(\Theta|\mathbf{Y}, \mathcal{M})$ is then replaced by a Dirac mass located at the estimate $\hat{\Theta}$. The posterior predictive distribution then simplifies to the likelihood function, and the lpd to $\ln \pi(\mathbf{Y}|\hat{\Theta}, \mathcal{M})$. Besides, directly using the number of inferred parameters – in this thesis the dimension of $\Theta$, $ND$ – is the simplest possible definition for the effective number of parameters. By combining a point estimate and this definition for the effective number of parameters, one obtains the classical Akaike information criterion (AIC) (Gelman et al., 2015, chapter 7)

$$ \widehat{\text{elpd}}_{\text{AIC}} = \widehat{\text{lpd}} - \hat{p}_{\text{AIC}} = \ln \pi (\mathbf{Y}|\hat{\Theta}, \mathcal{M}) - ND \quad (2.53) $$

which is usually defined with a factor -2, such that

$$ \text{AIC} = -2 \ln \pi (\mathbf{Y}|\hat{\Theta}, \mathcal{M}) + 2ND \quad (2.54) $$

**Example 2: the Bayesian information criterion (BIC)** – Limited to the maximum likelihood estimator $\hat{\Theta}_{\text{MLE}}$, the BIC (Schwarz, 1978) uses $\hat{p}_{\text{BIC}} = ND \ln(NL)/2$. The BIC is not exactly a Bayesian information criterion in the sense that it considers a point estimate $\hat{\Theta}_{\text{MLE}}$ and not a posterior distribution $\pi(\Theta|\mathbf{Y}, \mathcal{M})$. Its name “Bayesian” comes from the fact that it is an approximation of the log-evidence based on Laplace’s method (Konishi and Kitagawa, 2008, chapter 9).

For more information on elpd, cross validation and information criteria based on the posterior predictive distribution, see e.g., Gelman et al. (2015, chapter 7) or Vehtari et al. (2017). For a general review on information criteria, see, e.g., Konishi and Kitagawa (2008).

In general, an isolated value of Bayesian evidence lacks interpretability. This quantity is thus mostly used relatively, to compare multiple models. This is also true for information criteria. In particular, there exists no simple and general rule to decide whether a model $\mathcal{M}$ accurately reproduces observations.
### 2.3.2 Model checking using Bayesian hypothesis testing

When only one model $\mathcal{M}$ is considered, one might want to check whether it accurately reproduces the observations $\mathbf{Y}$, using hypothesis testing. One key advantage of this hypothesis testing compared to model selection or evaluation is that its metric, the $p$-value, is easily interpretable. In Bayesian statistics, such a hypothesis test is called *posterior predictive assessment*, or *checking*.

This subsection introduces posterior predictive assessment key notions, such as the underlying null hypothesis, the associated $p$-value, its Monte Carlo estimator, and the discrepancy measure $T$. Chapter 3 (Section 3.3.2) details three applications in interstellar medium studies. In Chapter 5 (Section 5.3), we extend the test to make it more robust to wrong decisions due to the Monte Carlo estimator error.

Posterior predictive checking was introduced in Guttman (1967) and Rubin (1984). Gelman et al. (1996) extended the definition to discrepancy measures $T$, more general than test statistics. It is based on the following null hypothesis: “The observation model $\mathcal{M}$ can reproduce the observations $\mathbf{Y}$”. Like model evaluation methods based on information criteria, this test relies on reproduced observations $\tilde{\mathbf{Y}}$ and on the posterior predictive distribution $\pi(\tilde{\mathbf{Y}}|\mathbf{Y})$. The goal is to compare predicted observations $\tilde{\mathbf{Y}}$ and the true observations $\mathbf{Y}$ with respect to a discrepancy measure $T: (\mathbf{Y}, \mathbf{\Theta}) \mapsto T(\mathbf{Y}, \mathbf{\Theta}) \in \mathbb{R}$. If the null hypothesis is true, then predictions $\tilde{\mathbf{Y}}$ should be as likely as the true observations $\mathbf{Y}$ with respect to the measure $T$. Common discrepancy measures include $T(\tilde{\mathbf{y}}, \mathbf{\Theta}) = \tilde{\mathbf{y}}$ for scalar observations (Gelman et al., 2015), the negative log-likelihood, and the $L_2$-norm

$$T(\tilde{\mathbf{Y}}, \mathbf{\Theta}) = \sum_{n=1}^{N} \sum_{\ell=1}^{L} \frac{(\tilde{y}_{n\ell} - f_{\ell}(\theta_n))^2}{\sigma_{n\ell^2}}, \quad (2.55)$$

that is also called $\chi^2$ loss in astrophysics articles.

The Bayesian $p$-value associated with this test corresponds to the probability of obtaining values $T(\tilde{\mathbf{Y}}, \mathbf{\Theta})$ at least as unlikely as $T(\mathbf{Y}, \mathbf{\Theta})$ under the null hypothesis.

$$p = \mathbb{P}_{(\tilde{\mathbf{Y}}, \mathbf{\Theta})} \left[ T(\tilde{\mathbf{Y}}, \mathbf{\Theta}) \ge T(\mathbf{Y}, \mathbf{\Theta}) \mid \mathbf{Y}, \mathcal{M} \right]. \quad (2.56)$$

Equivalently, it is the measure of the set $I = \{(\tilde{\mathbf{Y}}, \mathbf{\Theta}) \mid T(\tilde{\mathbf{Y}}, \mathbf{\Theta}) \ge T(\mathbf{Y}, \mathbf{\Theta})\}$ when using the model $\mathcal{M}$

$$p = \int \mathbb{1}_I (\tilde{\mathbf{Y}}, \mathbf{\Theta}) \pi (\tilde{\mathbf{Y}} | \mathbf{\Theta}, \mathcal{M}) \pi (\mathbf{\Theta} | \mathbf{Y}, \mathcal{M}) d\mathbf{\Theta} d\tilde{\mathbf{Y}}. \quad (2.57)$$

For classical choices of $T$, if the $p$-value is below a threshold $\alpha$ chosen prior to the analysis (typically 0.05 or 0.01), then the null hypothesis can be rejected with confidence $1-\alpha$. For some other choices of $T$ (as in Gelman et al. (2015, chapter 6)), the null hypothesis can be rejected when the $p$-value is below $\alpha/2$ or above $(1-\alpha)/2$.

Analytical computation of the Bayesian $p$-value can be performed in some simple cases (Meng, 1994). For instance, for a point estimate $\hat{\mathbf{\Theta}}$, i.e., a posterior distribution reduced to a Dirac on $\hat{\mathbf{\Theta}}$, the statistic in Eq. 2.55 follows a $\chi^2_{NL}$ distribution. The associated $p$-value,

$$p = \mathbb{P}_{T(\tilde{\mathbf{Y}}, \hat{\mathbf{\Theta}}) \sim \chi^2_{NL}} \left[ T(\tilde{\mathbf{Y}}, \hat{\mathbf{\Theta}}) \ge T(\mathbf{Y}, \hat{\mathbf{\Theta}}) \mid \mathbf{Y}, \mathcal{M} \right], \quad (2.58)$$

can thus be easily computed with the cdf of the $\chi^2$ distribution.

**Evaluation of the $p$-value** – Figure 2.3 illustrates the presented $p$-value on two examples with a scalar observation $\mathbf{Y} = y \in \mathbb{R}$. In both cases, the statistic $T$ is set to the negative log-likelihood. A point estimate $\hat{\mathbf{\Theta}}$ is considered for simplicity and visualization. The first row shows a Gaussian additive noise model. The test statistic thus boils down to the $\chi^2$ statistic (Eq. 2.55).
The second row corresponds to a more complicated fictitious observation model, defined as a Gaussian mixture model (GMM) with 3 components. The left column shows for both cases the likelihood distribution pdf on y with the parameter Θ fixed to the estimator ̂Θ. The black dashed line shows the attained likelihood pdf value for y, π(y|Θ). The blue area corresponds to the area to be integrated for the p-value computation (Eq. 2.57), i.e., the area for reproduced observations ̃y such that T(̃y, Θ) ≥ T(y, Θ). The right column shows the pdf on the test statistic T. In Figure 2.3b, as T ~ χ², the cdf can be evaluated efficiently. In Figure 2.3d, the considered test statistic T does not follow a simple distribution. Evaluating the cdf on the test statistic T exactly requires integrating with respect to ̃y. When the posterior distribution Θ is not reduced to a dirac, on and , which is unrealistic in high dimensions.

<figure><img src="image_19.png" /><figcaption>Figure 2.3 – Application of the Bayesian p-value (Eq. 2.57) for a point estimate ̂Θ.</figcaption></figure>

For these general cases, the p-value is approximated with a Monte Carlo estimator (Gelman et al., 1996). As in model selection methods based on information criteria, this Monte Carlo estimator exploits TMC samples Θ⁽ᵗ⁾ from the posterior, typically obtained with an MCMC algorithm. In addition, it requires sampling observation reproductions ̃Y⁽ᵗ⁾ from the observation model M and the sample Θ⁽ᵗ⁾. The p-value counts the number of iterations t where T(̃Y⁽ᵗ⁾, Θ⁽ᵗ⁾) ≥ T(Y, Θ⁽ᵗ⁾), i.e., it is an empirical frequency

$$ \hat{p} = \frac{1}{T_{\text{MC}}} \sum_{t=1}^{T_{\text{MC}}} \mathbb{1}_{I} (\tilde{\mathbf{Y}}^{(t)}, \mathbf{\Theta}^{(t)}), \quad \mathbf{\Theta}^{(t)} \sim \pi(\mathbf{\Theta}|\mathbf{Y}, \mathcal{M}), \quad \tilde{\mathbf{Y}}^{(t)} \sim \pi(\tilde{\mathbf{Y}}|\mathbf{\Theta}^{(t)}, \mathcal{M}), \quad (2.59) $$

where Θ⁽ᵗ⁾ ~ π(Θ|Y) indicates that Θ⁽ᵗ⁾ is drawn from the distribution of pdf π(Θ|Y). In an MCMC algorithm, the iterates from burn-in phase should be disregarded for the p-value evaluation, as they are not considered to be samples from the posterior distribution. However, their influence on the estimated p-value decreases as the size of the Markov chain increases. Note that unlike model selection methods based on information criteria or the three-case χ² rule widespread
in astrophysics, this hypothesis testing framework does not penalize the number of free parame-
ters. Overfitting cases thus cannot be detected, regardless of the test statistic $T(\mathbf{Y}, \mathbf{\Theta})$. Indeed,
in such a case, by definition of overfitting, the null hypothesis is verified.

**Marginal predictive checks** – As mentioned in Gelman et al. (2015, chapter 6), a more informative hypothesis testing approach can be applied on individual observations $\tilde{\mathbf{y}}_n$ instead of the full observation set $\tilde{\mathbf{Y}} = (\tilde{\mathbf{y}}_n)_{n=1}^N$. Such tests are called marginal predictive checks as they focus on the marginal predictive distribution $\pi(\tilde{\mathbf{y}}_n|\mathbf{Y}, \mathcal{M})$ instead of the joint predictive distribution $\pi(\tilde{\mathbf{Y}}|\mathbf{Y}, \mathcal{M})$. These tests can be used e.g., to find outliers among the observations $y_n$ or to simplify the identification of issues in a problematic model $\mathcal{M}$.

**Interpretability discussion** – One interpretability limitation of this Bayesian *p*-value is that its null distribution, i.e., the distribution of *p* (Eq. 2.57) when the null hypothesis is true, is not uniform on [0, 1] Bayarri and Berger (2000) and makes rejection rarer. This statistical behavior is illustrated in two examples in Gelman (2013) along with a discussion on the impact for the practical interest of the described *p*-value. Some alternative predictive posterior checks were calibrated to admit the uniform distribution asymptotically (Robins et al., 2000). These alternative checks are generally expensive to compute, while computing the presented *p*-value is very cheap. In Gelman et al. (2015, chapter 6), the authors emphasize that the described *p*-value is by construction a valid probability and should be used as such. In the following, we still use the presented *p*-value while acknowledging this limitation for interpretation.

In Chapter 3 (Section 3.3.2), we show applications of this model assessment approach in interstellar medium studies. Though the estimator $\hat{p}$ (Eq. 2.59) converges to the theoretical *p*-value (Eq. 2.57) with the law of large numbers, using an approximation instead of the theoretical value introduces an error. In Chapter 5 (Section 5.3), we extend the test by taking this error into account. The test is therefore made more robust to wrong decisions.

**2.4 Conclusion**

In this chapter, we covered the three main steps involved in solving an inverse problem: sta-
tistical modeling, statistical inference and model assessment.

Bayesian statistical modeling consists in setting the inverse problem to solve. The distribu-
tion of interest is called the posterior distribution. It is defined with the Bayes' theorem and
proportional to the product of the likelihood function and the pdf of the prior distribution. The
normalization constant, called the Bayesian evidence, is not considered in this thesis. The like-
lihood function encodes the observation model, which includes the forward model and the noise
model. The prior distribution encodes prior knowledge on the physical parameters to infer.

In Chapter 3 (Section 3.1), we review statistical modeling in interstellar medium studies. In Chapter 5 (Sections 5.1.1, 5.1.2, 5.1.3), we set the likelihood, prior and posterior of the general inverse problem addressed in this thesis, respectively.

For statistical inference, we described the two main approaches: optimization-based and
sampling-based. In both cases, the most common methods were recalled. The optimization-based
approach only provides point estimates. Conversely, the sampling-based approach, e.g., MCMC
algorithms, naturally yields uncertainty quantification on the physical parameters Θ. Due to the
absence of ground truth in ISM studies, we favor sampling methods to solve the inverse problem
considered in this thesis. Some advanced MCMC methods such as preconditioned Metropolis ad-
justed Langevin algorithm (PMALA), Gibbs sampling, and multiple-try Metropolis (MTM) were
covered.

Chapter 3 (Section 3.2) reviews applications of statistical inference in ISM studies. The sam-
pler proposed in Chapter 5, combines a PMALA kernel that performs efficient local explorations
and a second kernel that combines chromatic Gibbs and MTM. By exploiting the natural structure
of the parameter space, this second kernel can escape from local modes.

For model assessment, we covered Bayesian model selection and Bayesian model checking. The two main approaches of model selection, based on Bayesian evidence and information criteria, were explained. For model checking, we presented the Bayesian *p*-value from Gelman et al. (1996). Applications of both approaches in ISM studies are reviewed in Chapter 3 (Section 3.3). In ISM studies, often only one model is available, as generating forward models is expensive. In Chapter 5 (Section 5.3), we extend the model checking strategy by including uncertainties on the associated *p*-value.
# Appendix 2.A Samplers widespread in astrophysics dedicated to multimodal distributions

In astrophysics, the considered forward models are almost exclusively non-linear, even in very simple cases – see, e.g., applications in Chapter 3. This non-linearity causes the negative log-posterior to be non-convex with possibly multiple modes. Multimodal distributions are known to be difficult to sample, especially when modes are isolated, i.e., when they are separated by large low-probability regions. In practice, samplers such as MH, MALA, or HMC fail to explore the full distribution. Indeed, these samplers tend to get stuck in one local minimum, and would require an unrealistically large amount of samples to visit all modes.

This appendix presents some existing sampling algorithms dedicated to multimodal distributions. Some of these algorithms are already popular in the ISM community, as we outline in Chapter 3 (Section 3.2.2). We do not use the algorithms listed in this appendix to solve the considered inverse problem. In Chapter 5 (Section 5.2.2), we propose a sampling kernel that addresses multimodality with a combination of chromatic Gibbs and MTM. However, we will compare them with the proposed sampler in Chapter 5 (Section 5.4).

## 2.A.1 Adapting meta-heuristics to MCMC

In optimization, meta-heuristics are common methods to address multimodality. After introducing meta-heuristics and particularly tempering and population methods, we present some samplers adapted to the MCMC context. In particular, sequential MC (SMC) and the so-called affine invariant sampler, two popular methods in ISM studies – see Chapter 3 –, are presented.

**Meta-heuristics algorithms** were initially designed for discrete non-convex optimization problems such as the traveling salesman problem. They are stochastic in two ways: at each step *t*, they randomly generate candidates $\Theta_c$, and then keep them with a certain probability that depends on the loss value $\mathcal{L}(\Theta_c)$. The most popular meta-heuristics are simulated annealing (van Laarhoven and Aarts, 1987), genetic algorithms (Chelouah and Siarry, 2000), particle swarm (Eberhart and Kennedy, 1995; Kennedy and Eberhart, 1995), and some combinations of these algorithms (Kao and Zahara, 2008).

Genetic algorithm and particle swarm are both *population algorithms*. They first initialize a population of points, and then make the population evolve towards local minima with low loss values by combining current points or applying random mutations.

Simulated annealing resorts to *tempering*. Tempering relies on an inverse temperature parameter $\phi^{(t)} \in [0, 1]$. It first flattens the loss function by considering a tempered loss function $\Theta \mapsto \phi^{(t)}\mathcal{L}(\Theta)$, with $\phi^{(t)}$ close to 0. This allows simple transitions between modes and permits to consider points with high loss to escape from local minima. During the optimization procedure, the inverse temperature $\phi^{(t)}$ progressively increases to 1 and transitions become more conservative. Simulated annealing enjoys some strong theoretical properties, though only valid for extremely slow temperature evolutions. For a complete introduction to meta-heuristics for both discrete and continuous optimization problems, see e.g., Dreo et al. (2006).

**Tempering-based samplers without population** – The equi-energy sampler (Kou et al., 2006) and the adaptive parallel tempering algorithm (Miasojedow et al., 2013) exploit tempering, as simulated annealing. These two algorithms run parallel interacting Markov chains at different temperatures. High temperature chains can navigate between modes and only one chain at low temperature is actually used for estimations.

**Population-based samplers without tempering** – Some other samplers adapt evolutionary methods such as the genetic algorithms to the MCMC context. For instance, Evolutionary MC
(EMC) (Liang and Wong, 2000), distributed genetic MC (DGMC) (Hu and Tsui, 2010) do not include any tempering, but run parallel interacting chains to generate relevant candidates. Besides, the so-called affine-invariant MCMC sampler (Goodman and Weare, 2010) is a special case of EMC that only applies the snooker crossover (Liang and Wong, 2001). To the best of our knowledge, this link between EMC and the affine-invariant sampler was never identified before.

The affine-invariant sampler runs $K > 1$ interacting Markov chains in parallel. At each step $t$, it updates an iterate $\Theta^{(k,t-1)}$ by making it interact with another iterate $\Theta^{(j,t-1)}$ with $j \neq k$. First, it draws a scaling from the distribution of pdf

$$ \pi(z) \propto \frac{1}{\sqrt{z}} \mathbb{I}_{[\frac{1}{2}, 2]}(z). \qquad (2.60) $$

Second, it proposes a candidate $\Theta_c^{(k,t)}$ using an affine combination of $\Theta^{(k,t-1)}$ and $\Theta^{(j,t-1)}$

$$ \Theta_c^{(k,t)} = \Theta^{(j,t-1)} + z (\Theta^{(k,t-1)} - \Theta^{(j,t-1)}). \qquad (2.61) $$

The candidate is then accepted with probability

$$ \rho_{\text{emcee}}^{(t)} = \min \left( 1, z^{ND-1} \frac{\pi(\Theta_c^{(k,t)} | \mathbf{Y})}{\pi(\Theta^{(k,t-1)} | \mathbf{Y})} \right). \qquad (2.62) $$

Algorithm 2.5 summarizes the affine-invariant sampler.

<caption>**Algorithm 2.5: Affine-invariant sampler**</caption>

**Input:** Number of parallel chains $K$, starting points $(\Theta^{(k,0)})_{k=1}^K$, number of iterations $T_{BI}$ et $T_{MC}$
**Output:** Set of samples $\{\Theta^{(k,t)}, t = T_{BI} + 1, \dots, T_{MC}, k = 1, \dots, K\}$

1 **for** $t = 1, \dots, T_{MC}$ **do**
2 &nbsp;&nbsp;&nbsp;**for** $k = 1, \dots, K$ **do**
3 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**// generate candidate**
4 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$j \sim \text{Unif}(\{1, \dots, k-1, k+1, \dots, K\})$
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$z \sim \pi(z)$
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**// using Eq. 2.60**
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\Theta_c^{(k,t)} = \Theta^{(j,t-1)} + z (\Theta^{(k,t-1)} - \Theta^{(j,t-1)})$
8 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**// snooker crossover (Eq. 2.61)**
9 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**// Accept-reject step**
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\rho_{\text{emcee}}^{(t)}$
11 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\zeta \sim \text{Unif}(0, 1)$
12 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\Theta_c^{(k,t)} = \Theta_c^{(k,t)}$ **if** $\zeta \le \rho_{\text{emcee}}^{(t)}$ **else** $\Theta_c^{(k,t-1)}$

This affine-invariant MCMC sampler – and the associated EMCEE package (Foreman-Mackey et al., 2013) – is a common sampler in astrophysics, as we state in Chapter 3 (Section 3.2.2). In principle, the snooker crossover enables to jump between modes. In Foreman-Mackey et al. (2013), the authors state the affine-invariant sampler is in practice not suited to multimodal cases. They recommend initializing the chains in a small ball centered around a point expected to be close to the MAP to avoid getting trapped in a low probability local mode.

**Sequential MC (SMC): combining tempering and population methods – SMC (Del Moral et al., 2006)** combines tempering with a population approach. It was originally defined as a generalization of Kalman filters to non-linear or non-Gaussian state spaces. In this context, SMC, also known as Particle filter (Del Moral, 2004), draws samples from a succession of distributions. It was adapted to multimodal sampling by defining a succession of $K$ tempered posterior
distributions

$$ \pi_k(\Theta|\mathbf{Y}) \propto [\pi(\mathbf{Y}|\Theta)]^{k/K} \pi(\Theta). \qquad (2.63) $$

It first generates $T_{MC}$ samples from the prior distribution. These samples are then modified as the likelihood is progressively included in the tempered posterior by increasing an inverse temperature parameter $k/K \in [0, 1]^5$. Other sequences of inverse temperature can also be used.

At each step $k$, in a first correction step, the weights $w^{(k,t)}$ of the $T_{MC}$ samples $\Theta^{(k-1,t)}$ are evaluated. Then, in a selection step, the population is resampled according to the weights $\{w^{(k,t)}\}_{t=1}^{T_{MC}}$. Finally, in a mutation step, each member of the resampled population goes through a few steps of the MH algorithm with a kernel $K_k$ that admits the tempered posterior $\pi_k(\cdot|\mathbf{Y})$ as invariant distribution. This last step can be performed in parallel. Algorithm 2.6 summarizes the SMC algorithm. As a by-product, SMC can also yield the Bayesian evidence, which is useful for model selection – see Section 2.3.

<caption>**Algorithm 2.6: Sequential Monte Carlo (SMC)**</caption>

**Input:** number of intermediate steps $K$, number of samples $T_{MC}$, number of MH steps $N_{MH}$, transition kernels $K_k$
**Initialization:** Starting points $\Theta^{(0,t)} \sim \pi(\Theta)$ for $t=1, \dots, T_{MC}$
**2** for $k=1, \dots, K$ do
    // Correction: compute weights
    3 & $\tilde{w}^{(k,t)} = \frac{\pi_k(\mathbf{Y}|\Theta^{(k-1,t)})}{\pi_{k-1}(\mathbf{Y}|\Theta^{(k-1,t)})}$ for $t=1, \dots, T_{MC}$
    4 & $w^{(k,t)} = \frac{\tilde{w}^{(k,t)}}{\sum_{t=1}^{T_{MC}} \tilde{w}^{(k,t)}}$ for $t=1, \dots, T_{MC}$ // normalize the weights
    5
    // Selection: resample the population
    6 & sample $\tilde{\Theta}^{(k,t)}$ from $\{\Theta^{(k-1,t)}\}_{t=1}^{T_{MC}}$ with selection probabilities $\{w^{(k,t)}\}_{t=1}^{T_{MC}}$
    7 & // Mutation: apply $N_{MH}$ steps of MH
    & $\Theta^{(k,t)} \sim K_k^{N_{MH}}(\tilde{\Theta}^{(k-1,t)}, \cdot)$ for $t=1, \dots, T_{MC}$
    **Output:** Set of samples $\{\Theta^{(K,t)}\}_{t=1}^{T_{MC}}$

## 2.A.2 Prior or parallel identification of the modes

Some methods resort to an augmented distribution with a latent mode index. The posterior distribution is approximated by a mixture model on $\Theta$, where each model corresponds to a mode. Such methods sample from the posterior using two kernels: a local kernel that explores around a mode, and a jump kernel that permits jumps between already identified modes. Such methods include darting MC (DMC) (Andricioaei et al., 2001), jumping adaptative multimodal sampler (JAMS) (Pompe et al., 2020), regeneration darting MC (RDMC) (Ahn et al., 2013) and wormhole Hamiltonian Monte Carlo (WHMC) (Lan et al., 2014). WHMC is a particular case of the Riemannian manifold HMC algorithm introduced in Section 2.2.2.5. The metric of the corresponding manifold combines the standard Euclidean distance and a wormhole metric that shortens the distances between already identified modes, which simplifies transitions from one to another.

DMC and JAMS require a prior identification of the distribution modes by some optimization methods. RDMC and WHMC allow running optimization methods in parallel to the sampler and

<footnote n="5">5. The presented tempering approach is called likelihood tempering. Other types of tempering exist, such as data tempering or model tempering (Mlikota and Schorfheide, 2022).</footnote>
update the distribution and the sampler parameters at so-called random regeneration times (Gilks et al., 1998). However, in high-dimensional settings and with a non-linear forward model, the posterior has potentially many modes with only a few of significant weight in the mixture. The identification of the relevant modes with standard optimization methods remains difficult in this context.

### 2.A.3 Nested sampling

Nested sampling (Skilling, 2004; Skilling, 2006) is a class of algorithms whose primary goal is to compute the Bayesian evidence and that produces posterior samples as a by-product. Nested sampling is therefore not an MCMC algorithm. To evaluate the evidence $\pi(\mathbf{Y}|\mathcal{M})$, here noted $Z$, it considers the prior mass associated with a likelihood value greater than $v \ge 0$

$$x(v) = \int_{\pi(\mathbf{Y}|\mathbf{\Theta})>v} \pi(\mathbf{\Theta}) d\mathbf{\Theta}. \qquad (2.64)$$

By construction, the function $x$ decreases on $[0, \pi(\mathbf{Y}|\hat{\mathbf{\Theta}}_{MLE})]$. It verifies $x(0) = 1$ and $x(\pi(\mathbf{Y}|\hat{\mathbf{\Theta}}_{MLE})) = 0$. Besides, a range $[x, x + dx]$ corresponds to the prior weight associated with likelihood values in $[v, v + dv]$. The evidence $Z$ can then be rewritten as the uni-dimensional integral

$$Z = \int_{0}^{1} v(x) dx, \qquad (2.65)$$

where $v: x \mapsto v(x)$ is the inverse function of $x: v \mapsto x(v)$. As $x$ is decreasing on $[0, \pi(\mathbf{Y}|\hat{\mathbf{\Theta}}_{MLE})]$, so is $v$ on $[0, 1]$.

In its first version Skilling (2004); Skilling (2006), nested sampling considers a population of $K$ points $\mathbf{\Theta}^{(k,t)}$. At each iteration $t$, the point $\mathbf{\Theta}^{(i,t-1)}$ with lowest likelihood $v_t$ is exploited to update the Bayesian evidence estimator $\hat{Z}^{(t)}$. This estimator implements the rectangle approximation of an integral:

$$\begin{aligned} \hat{Z}^{(t)} &= \sum_{\tau=1}^{t} v_{\tau} \left[ \exp \left( -\frac{\tau-1}{T_{\mathrm{MC}}} \right) - \exp \left( -\frac{\tau}{T_{\mathrm{MC}}} \right) \right] \\ &= \hat{Z}^{(t-1)} + v_t \left[ \exp \left( -\frac{t-1}{T_{\mathrm{MC}}} \right) - \exp \left( -\frac{t}{T_{\mathrm{MC}}} \right) \right]. \end{aligned} \qquad (2.66)$$

Then, the point $\mathbf{\Theta}^{(i,t-1)}$ is updated to $\mathbf{\Theta}^{(i,t)}$ with MCMC steps. The other points in the population are passed as is to the next iteration, i.e., $\mathbf{\Theta}^{(k,t)} = \mathbf{\Theta}^{(k,t-1)}$, for all $k \neq i$. In these steps, the transition kernel $K$ samples from the prior distribution with the constraint on the likelihood value $\pi(\mathbf{Y}|\mathbf{\Theta}^{(k)}) > v_t$. Therefore, in nested sampling, the main difficulty lies in sampling from the prior with a hard lower-level constraint on the likelihood values. Algorithm 2.7 summarizes the procedure of nested sampling. Figure 2.4 illustrates the principle of this original version of nested sampling in a two-dimensional case. In particular, note that higher population sizes $K$ lead to horizontally finer rectangles, and that higher number of iterations $T_{MC}$ permit to explore high likelihood regions.

In more recent papers, better sampling strategies were proposed including ellipsoidal rejection sampling (Feroz and Hobson, 2008) - and the associated code MULTINEST (Feroz et al., 2009), diffusive sampling (Brewer et al., 2011) and slice sampling (Handley et al., 2015) - and the associated code POLYCHORD. For a general review on nested sampling, see e.g., Buchner (2023).

Nested sampling is applicable for distributions with dimension up to $\mathcal{O}(10^3)$ for the latest methods. One proximal nested sampling algorithm was proposed in Cai et al. (2022) to be applied to distributions on very high dimensional spaces - up to $\mathcal{O}(10^6)$ dimensions. However, this algorithm was built for imaging analysis and thus is only valid for log-concave posterior distributions. This property is seldom verified in astrophysical studies, due to the non-linearity of the forward model.
<figure><img src="image_20.png" /><figcaption>Figure 2.4 – Illustration of nested sampling principle in a 2-dimensional parameter space with $T_{\text{MC}} = 4$ steps and a population of size $K = 2$. (Left) Parameter space with a population of $K = 2$ points $\theta^{(k,t)} \in \mathbb{R}^2$. The gray lines space represent the associated likelihood isocontours in the parameter space, with darker lines corresponding to higher likelihood. At each step $t$, one of the two points is selected to update the evidence estimator $\hat{Z}^{(t)}$. The selected point is then updated: the dashed arrows show their trajectories. (Right) illustration of the integral in Eq. 2.65, with $x_t$ the prior mass associated with a likelihood greater than $v_t$, as defined in Eq. 2.64. The true evidence is the area under the dashed curve. The evidence estimator $\hat{Z}^{(t)}$ from Eq. 2.66 is the grayed area.</figcaption></figure>

**Algorithm 2.7: Nested sampling**

**Input:** population size $K$, number of iterations $T_{\text{MC}}$, sampling kernel $\mathcal{K}$, number of MCMC iterations $N_{\text{MCMC}}$

1 **Initialization:** Sample $K$ points from the prior $\Theta^{(k,0)} \sim \pi(\Theta)$, set of samples $S^{(0)} = \emptyset$, Bayesian evidence estimator $\hat{Z}^{(0)} = 0$

2 **for** $t = 1, \dots, T_{\text{MC}}$ **do**
3 &nbsp;&nbsp;&nbsp; $i = \arg \min_k \pi(\mathbf{Y}|\Theta^{(k,t-1)})$ // select lowest likelihood
4 &nbsp;&nbsp;&nbsp; $v_t = \pi(\mathbf{Y}|\Theta^{(i,t-1)})$
5 &nbsp;&nbsp;&nbsp; $x_t = \exp(-t/T_{\text{MC}})$
6 &nbsp;&nbsp;&nbsp; $w_t = x_{t-1} - x_t$ // weight of sample $\Theta^{(k)}$ in the set of samples
7 &nbsp;&nbsp;&nbsp; $\hat{Z}^{(t)} = \hat{Z}^{(t-1)} + v_t w_t$ // update evidence estimator (Eq. 2.66)
8 &nbsp;&nbsp;&nbsp; $S^{(t)} = S^{(t-1)} \cup \{(\Theta^{(i,t-1)}, w_t)\}$ // update set of samples
9 &nbsp;&nbsp;&nbsp; $\Theta^{(i,t)} \sim \mathcal{K}^{N_{\text{MCMC}}}(\Theta^{(i,t-1)}, \cdot)$ // update $\Theta^{(k)}$ to increase its likelihood
10 &nbsp;&nbsp;&nbsp; $\Theta^{(k,t)} = \Theta^{(k,t-1)}$ for all $k \neq i$ // the other points are not updated

**Output:** Bayesian evidence estimator $\hat{Z}^{(T_{\text{MC}})}$, set of samples $S^{(T_{\text{MC}})}$
# Chapter 3

## Review of statistical inference in ISM studies and position of our problem

"No computer is ever going to ask a new, reasonable question. It takes trained people to do that."

Grace Hopper

### Contents

<table><tr><td>3.1 Statistical modeling</td><td>62</td></tr><tr><td>3.1.1 Forward model: handling a numerical model</td><td>62</td></tr><tr><td>3.1.1.1 Direct evaluations or using precomputed models</td><td>62</td></tr><tr><td>3.1.1.2 Datasets structure</td><td>63</td></tr><tr><td>3.1.1.3 Learning an emulator</td><td>64</td></tr><tr><td>3.1.2 Noise model and forward model misspecification</td><td>65</td></tr><tr><td>3.1.2.1 Notation</td><td>65</td></tr><tr><td>3.1.2.2 Purely Gaussian noise models</td><td>66</td></tr><tr><td>3.1.2.3 Combining additive and multiplicative noises</td><td>66</td></tr><tr><td>3.1.2.4 Including censorship in the likelihood function</td><td>68</td></tr><tr><td>3.1.2.5 Avoiding the manual specification of a noise model</td><td>70</td></tr><tr><td>3.1.3 Prior distributions and regularization functions</td><td>70</td></tr><tr><td>3.2 Statistical inference</td><td>72</td></tr><tr><td>3.2.1 Optimization-based inference</td><td>72</td></tr><tr><td>3.2.2 Sampling-based inference</td><td>72</td></tr><tr><td>3.3 Comparing and checking observation models</td><td>74</td></tr><tr><td>3.3.1 Model comparison</td><td>74</td></tr><tr><td>3.3.2 Model checking with posterior predictive assessment</td><td>75</td></tr><tr><td>3.4 Our observation model</td><td>76</td></tr><tr><td>3.4.1 Noise model on the original hyperspectral cube</td><td>76</td></tr><tr><td>3.4.2 From the hyperspectral cube to maps of integrated intensities</td><td>79</td></tr><tr><td>3.4.3 Forward model, observational effects and censorship</td><td>80</td></tr><tr><td>3.5 Modeling, inference and model assessment choices</td><td>81</td></tr><tr><td>3.6 Conclusion</td><td>82</td></tr></table>
In Chapter 2, we described some statistical modeling and inference methods, including model selection and model assessment. This chapter reviews the main statistical approaches adopted in ISM studies. As we will show, ISM studies involving integrated intensities of atomic or molecular emission lines are mostly based on grid search. Few sampling-based studies have been proposed, limited to low-dimensional settings. Conversely, cosmology was among the first astrophysics fields to widely adopt sampling methods and to use advanced models. For instance, many works published in the early 2000's infer cosmological parameters from observations of the cosmological microwave background, sometimes combined with other observations. Therefore, we expand the review to stellar physics and cosmology, where the degree of adoption of Bayesian statistical modeling and inference is higher.

The goal of this chapter is to detail existing applications in ISM studies at each step of the modeling or inference to be able to motivate our contributions. It isolates the different steps and components of the Bayesian methodology. Section 3.1 describes the statistical modeling tendencies, including how numerical forward models are handled and the choice of noise model and of prior distribution. Section 3.2 lists popular methods for both optimization-based and sampling-based inference. Section 3.3 reviews applications of model selection and model checking. Finally, Section 3.6, briefly compares our contributions with the literature.

## 3.1 Statistical modeling

As presented in Chapter 2 (Section 2.1), statistical modeling aims at defining a posterior distribution $\pi(\Theta|Y)$. The posterior involves a prior distribution $\pi(\Theta)$ and a likelihood function $\pi(Y|\Theta)$. In turn, the likelihood function is defined from a forward model – usually a numerical model in astrophysics – and a noise model. In this section, we discuss choices and standards in the ISM community and some related communities on the forward model, the noise model and the prior distribution.

### 3.1.1 Forward model: handling a numerical model

For any numerical model, small changes in the physical parameters can lead to very different predicted observables. Adopting a realistic forward model f that predicts observables from physical parameters is critical for inference. Codes that model the observed environment more realistically lead to more meaningful estimations. However, the complexity of the physics considered in a code directly impacts its evaluation time and hence its usability. This is especially true for inference, which requires many evaluations of the forward model.

#### 3.1.1.1 Direct evaluations or using precomputed models

There are two main approaches to handle a numerical forward model in inference, depending on their speed. The first approach consists in running the code during the inversion procedure. This is only possible for fast numerical codes. The second approach targets slower codes. It consists in generating a grid of precomputed models prior to running the inference procedure.

**Running the code during the inversion** – On the one hand, the lighter codes can run in a few seconds and are often used directly in optimization or Bayesian inference. In ISM studies, many inverse problems rely on a fast model, such as those listed in Chapter 1 (Section 1.2). Combinations of these codes can also be used. For instance, Behrens et al. (2022) and Keil et al. (2022) rely on a combination of RADEX and ULCCHEM – see Chapter 1 (Sections 1.2.1 and 1.2.2). Even for the fastest codes, the evaluation of the numerical code is usually the slowest step in the inversion procedure. In astrochemistry, many efforts were spent on reducing the execution time. Some large chemical networks were simplified to smaller ones (Holdship et al., 2018). Alternatively, the original numerical model is used directly in the inversion, but the number of
evaluations is limited by progressively reusing past evaluations (Keil et al., 2022).

**Generating a dataset of precomputed models prior to the analysis** – On the other hand, a more comprehensive model such as the Meudon PDR code Le Petit et al. (2006) – introduced in Chapter 1 – that handles multiple physical processes on a 1D spatial grid typically requires several hours of computations. These durations are prohibitively long for inferring physical parameters on large observation maps. These heavier numerical codes are frequently evaluated on a grid {$\theta_i$} on the parameter space prior to any analysis. Following the ISM community standard, we call precomputed model one evaluation $f(\theta_i)$ of a numerical model $f$ prior to any analysis. Besides, we call a set of $(\theta, f(\theta))$ pairs a dataset of precomputed models or grid of precomputed models interchangeably, whether the dataset has a regular grid structure or not.

Datasets of precomputed models can be used for grid search in optimization approach or to approximate integrals with Bayesian-like approaches – see Section 3.2. They can also be used to derive faster emulators of the numerical model.

We now discuss the choice of the structure of the dataset, which can impact both uses, and strategies to derive emulators.

### 3.1.1.2 Datasets structure

A dataset of precomputed models can either be built at once, prior to any analysis, or iteratively until an accuracy criterion is satisfied. In Chapter 4 (Section 4.3), we resort to the first option. We build a dataset of precomputed models at once using a lattice structure.

**Generation of the dataset prior to any analysis** – The first approach is widespread in ISM studies. Most datasets of precomputed models are structured as lattices, i.e., as uniform regular grids, or as near-uniform grids with few and light differences in spacings. For instance, Sheffer et al. (2011) resorts to a lattice structure to perform a grid search in an optimization approach. Similarly, Pérez-Montero (2014) and Blanc et al. (2015) use this structure to approximate integrals in a Bayesian approach. Finally, Ramambason et al. (2022) and Smirnov-Pinchukov et al. (2022) exploit this structure to train emulators of their numerical models based on a nearest-neighbor interpolation and a $k$-nearest-neighbor regression, respectively.

Although very common, uniform or near-uniform grids are not systematically used. For instance, to derive an emulator of the Meudon PDR code, Wu et al. (2018) resorts to a non-uniform grid of models with more points in the region of interest in the parameter space. This non-uniformity ensures good accuracy of the emulator in this region. Besides, Bron et al. (2021) generates two datasets using independent and identically distributed (i.i.d.) draws from the uniform distribution on a cube defined with lower and upper values in the physical parameter space.

In cosmology, datasets are often generated using the latin hypercube sampling (LHS) algorithm – to be introduced in Chapter 4. For instance, see Spurio Mancini et al. (2022); Mootoovaloo et al. (2022); Agarwal et al. (2012). Some of its variants are also popular Heitmann et al. (2009); Ramachandra et al. (2021).

**Iterative construction of the dataset** – There are only two cases in ISM studies, to the best of our knowledge:

- In Vale Asari et al. (2016), the authors start from a coarse lattice for a first approximation, and refine it with an octree. An octree is a tree structure that represents the first grid as its first "layer". Regions of the parameter space in which the emulator is not accurate enough are divided into smaller regions and explored with finer grids. See Saftly et al. (2013) for an introduction to octrees with an application in a Monte Carlo radiative transfer model.

- In Galliano (2018), the authors start from a coarse grid and progressively refine it using a mid-point strategy until they satisfy an accuracy criterion.
For a given size of the dataset, the iterative approach usually takes much more time compared to the first approach, but can lead to better performance by focusing on challenging regions of the parameter space. In other words, a non-iterative approach can lead to oversampled parameter spaces, i.e., to more evaluations than needed to obtain a good approximation of the original model. For large grids, this oversampling can result in high storage and computational costs. In particular, the grid might be too large to be loaded in memory.

**Our proposal** – In Chapter 4 (Section 4.3), we present the results obtained with a lattice structured dataset of precomputed models of the Meudon PDR code. This dataset was generated to derive an emulator of the Meudon PDR code.

### 3.1.1.3 Learning an emulator

When the numerical model is too slow to be used directly in an inversion procedure, it is common to derive faster emulators. In this thesis, an *emulator* of a numerical code is an approximation of the relation between its inputs $\theta$ and its outputs $f(\theta)$, derived from a dataset of precomputed models. Emulating the numerical model, e.g., with interpolation methods, permits to predict observables for new points $\theta$ with lower evaluation time. In particular, an emulator does not perform the same intermediate computations as the original code.

Emulating the full likelihood function, which includes both the numerical model and the observation uncertainty model, is possible but quite rare. Indeed, it requires specifying a noise model, and is therefore either observation-specific or very generic. For instance, the SKYNET artificial neural network (ANN) (Graff et al., 2014) emulates the full likelihood function by approximating a numerical model and by assuming a noise model. The noise model is considered Gaussian and uncorrelated with fixed variance for continuous variable inference. The likelihood function is set to the cross-entropy (Bishop, 2006, chapter 4) for classification tasks. For a short introduction to ANNs, see Chapter 4 (Section 4.2.1). In this thesis, we focus on numerical model emulation to handle separately the forward model and the noise model.

In ISM studies, many numerical model emulations have already been proposed, with different choices of class of functions to define the emulator:

- Interpolation methods are well spread in the ISM community. For instance, linear interpolation was used in Galliano (2018) to emulate a combination of simple dust codes. Linear interpolation was also used in Ramambason et al. (2022), combined with nearest-neighbor interpolation to emulate CLOUDY, a HII region model. Similarly, radial basis function interpolation are used in Wu et al. (2018) to emulate the Meudon PDR code. Spline interpolation methods are sometimes used to compute integrals in Bayesian approaches (Blanc et al., 2015), but not to emulate numerical models in the ISM community. However, they are used for emulator derivation in other communities. For instance, Bailer-Jones (2011) emulates a simulator using a thin-plate spline to infer a set of stellar parameters.

- Classic machine learning regression approaches are punctually used. For instance, in Smirnov-Pinchukov et al. (2022), a *k*-nearest-neighbor regression emulates a protoplanetary disks model. In Bron et al. (2021), a random forest emulates a chemistry model.

- ANNs form a versatile class of functions that is often preferred to address the complexity of comprehensive ISM models. For instance, ANN emulators of astrochemical models are derived in de Mijolla et al. (2019); Holdship et al. (2021); Grassi et al. (2022). In addition, in Grassi et al. (2011), the authors derive a new simulation code and an associated ANN emulator.

In cosmology, deriving emulators of numerical codes such as CAMB that compute power spectra has received a lot of attention. Polynomial regression was used at first (Jimenez et al., 2004; Fendt and Wandelt, 2007). Since 2007, Gaussian process regression and ANNs are the
two main approaches. Gaussian processes are used, e.g., in Heitmann et al. (2009); Ramachandra et al. (2021); Mootoovaloo et al. (2022), and ANNs e.g., in Auld et al. (2007); Agarwal et al. (2012); Manrique-Yus and Sellentin (2019). The current state-of-the-art emulators, COSMOPower (Spurio Mancini et al., 2022), are fast ANNs that achieve an estimated mean error of 0.4% on power spectra. These state-of-the-art emulators were derived to be used as fast surrogate models in MCMC algorithms for Bayesian parameter inference.

**Our proposal** – We resort to ANNs to approximate the Meudon PDR code. We chose this family of functions as 1) it enjoys strong theoretical properties, 2) it proved state-of-the-art in a variety of regression problems, including in ISM studies, and 3) it allows to easily compute first and second order derivatives using auto-differentiation. Chapter 4 (Section 4.6) presents the design and training strategies we proposed to address specificities of complex ISM numerical models.

### 3.1.2 Noise model and forward model misspecification

Once the forward model is selected, an uncertainty model on the observations needs to be specified. Indeed, it is highly unlikely for a forward model to reproduce exactly the observations, unless over-parameterized. The noise model defines the likelihood function, which quantifies how distant a prediction $f(\theta)$ is from the corresponding observation $y$. In this subsection, we describe existing noise models used in ISM studies.

The noise model of the inverse problem considered in this thesis will be fully introduced in Section 3.4. It involves an additive thermal Gaussian noise and a multiplicative lognormal noise. The multiplicative uncertainty source encodes both calibration errors and forward model mis-specification. In addition to noise, observations may contain censored observations due to the sensitivity limits of the telescopes.

#### 3.1.2.1 Notation

For this subsection on noise models, we extend the notation, in particular $\mathbf{Y}$ and $\mathbf{\Theta}$. In this subsection, the $L$ observables are not restricted to the integrated intensity of emission lines, and may correspond to the flux in spectral channels. Similarly, the $N$ components are not restricted to pixels of a map, but correspond to the more general notion of observation beams, which does not assume a map structure. Maps of $N$ pixels and with integrated intensities of $L$ ionic, atomic or molecular emission lines are thus a special case of this extended notation. In general, the observation model can be written

$$ \forall n \in [1, N], \quad \forall \ell \in [1, L], \quad y_{n\ell} = \mathcal{A}(I_{n\ell}), \qquad (3.1) $$

where $I_{n\ell}$ is the original true and unaltered value for observable $\ell$ and beam $n$ and $\mathcal{A}$ is a general observation operator. The observation operator $\mathcal{A}$ can include measurement noise, calibration error, modeling error from the choice of forward model, or censorship, i.e., upper bounds on observables. The number of beams ranges from $N = 1$, e.g., in Pérez-Montero (2014); Keil et al. (2022), while some others consider multiple beams, e.g., $N = 176$ pixels in Wu et al. (2018) and $N = 798$ galaxies in Galliano et al. (2021). Similarly, depending on the observed signal, the number $L$ of observables per beam can be $\mathcal{O}(10^2)$ to $\mathcal{O}(10^3)$ for non spectrally integrated signals and $\mathcal{O}(10^0)$ to $\mathcal{O}(10^1)$ for line integrated intensities.

Inverse problem usually rely on the assumption that the true signals $I_{n\ell}$ can be reproduced exactly and simultaneously for all $n$ and $\ell$ by a forward model $f$ for a physical parameter $\mathbf{\Theta} = (\theta_n)_{n=1}^N$. In other words, it is assumed that there exists a physical parameter $\mathbf{\Theta}$ such that for all $n$ and $\ell$, $I_{n\ell} = f_\ell(\theta_n)$. The observation model (Eq. 3.1) is therefore rewritten

$$ \forall n \in [1, N], \quad \forall \ell \in [1, L], \quad y_{n\ell} = \mathcal{A}(f_\ell(\theta_n)). \qquad (3.2) $$
In the following, we list some observation operators A used in ISM studies, including noise models and censorship.

### 3.1.2.2 Purely Gaussian noise models

Gaussian and additive uncertainty models are widespread in ISM studies:

$$y_{nl} = f_{l}(\theta_{n}) + \varepsilon_{nl}^{(a)}, \quad \varepsilon^{(a)} \sim \mathcal{N}(0, \Sigma^{(a)}), \qquad (3.3)$$

where $\varepsilon^{(a)}$ is a measurement noise with zero mean and covariance matrix $\Sigma^{(a)}$. Such a model is often a good approximation for noise models thanks to the central limit theorem. Besides, it is extremely simple to manipulate. For these two reasons, Gaussian distributions are a by-default choice in statistical modeling in ISM inverse problems (see e.g., Galliano et al. (2003); Chevallard et al. (2013); Pérez-Montero (2014); Chevance et al. (2016); Wu et al. (2018); Lee et al. (2019); Roueff et al. (2021); Keil et al. (2022); Behrens et al. (2022)). In all these works, the covariance matrix $\Sigma^{(a)}$ is diagonal, i.e., all the noise components are assumed independent. In Galliano et al. (2003); Chevance et al. (2016); Lee et al. (2019), this observation model underlies the choice of the $\chi^2$ as a loss function:

$$\mathcal{L}(\Theta) = \chi_{\nu}^{2}(\Theta) = \sum_{n=1}^{N} \sum_{l=1}^{L} \frac{(f_{l}(\theta_{n}) - y_{nl})^{2}}{\sigma_{nl}^{2}}, \qquad (3.4)$$

where $\nu$ is a so-called degree of freedom. This choice of notation comes from the fact that a sum of $\nu$ i.i.d. squared Gaussian random variables follows a $\chi^2$ distribution. The degree of freedom is often distinct from the dimension of $\Theta$, even for linear models f (Andrae et al., 2010). In other words, in general, $\nu \neq NL$. In some papers (e.g., Vale Asari et al. (2016)), the Gaussian noise is considered on the log of the observations $\log y_{nl}$, which is equivalent to a multiplicative noise $\varepsilon^{(m)}$ following a lognormal distribution. A non-diagonal covariance matrix $\Sigma^{(a)}$ permits to account for correlation in Gaussian noise models. This assumption appears, e.g., in star property inference (Bailer-Jones, 2011) where $\Sigma^{(a)}$ is known.

### 3.1.2.3 Combining additive and multiplicative noises

It is common in ISM studies to combine two sources of noise to account for thermal noise, calibration noise, or model misspecification. Recall that in the inverse problem considered in this thesis, we consider a Gaussian additive thermal noise and a multiplicative lognormal noise. This latter noise includes both calibration error and model misspecification.

**Gaussian additive and Gaussian multiplicative distributions** – Some ISM studies combine a Gaussian additive noise $\varepsilon^{(a)}$ and a Gaussian multiplicative noise $\varepsilon^{(m)}$. This multiplicative error can represent calibration noise (Gordon et al., 2014; Ciurlo et al., 2016; Galliano, 2018; Galliano et al., 2021) or model misspecification (Blanc et al., 2015; Vale Asari et al., 2016; Jóhannesson et al., 2016). The observation model then reads

$$y_{nl} = \varepsilon_{nl}^{(m)} f_{l}(\theta_{n}) + \varepsilon_{nl}^{(a)}, \quad \varepsilon^{(m)} \sim \mathcal{N}(1, \Sigma^{(m)}), \quad \varepsilon^{(a)} \sim \mathcal{N}(0, \Sigma^{(a)}). \qquad (3.5)$$

The covariance matrices $\Sigma^{(a)}$ are assumed known and diagonal in all these papers. In Vale Asari et al. (2016) and Jóhannesson et al. (2016), the multiplicative error $\varepsilon_{nl}^{(m)}$ is a single scalar value equal for all observations, and handled as a nuisance parameter, i.e., inferred and marginalized. In Blanc et al. (2015), the multiplicative uncertainty source $\varepsilon_{nl}^{(m)}$ is also a unique scalar with fixed variance. When the observation $\mathbf{Y}$ is a multispectral map, a hyperspectral map, or a set of multichannel observations, the covariance matrix $\Sigma^{(m)}$ is assumed known. In such cases, it needs to be specified both spatially and spectrally. Ciurlo et al. (2016) assumes no correlation,
i.e., independent realizations for all $n$ and $\ell$. Gordon et al. (2014) considers a partial spectral correlation and no spatial correlation. Galliano (2018); Galliano et al. (2021) consider a partial spectral correlation and a spatial correlation of exactly 1.

Using a Gaussian model for multiplicative noise simplifies computation, as the overall uncertainty model is Gaussian. However, depending on the standard deviation of the multiplicative error, a Gaussian model may not accurately account for multiplicative uncertainty. A typical model for multiplicative noise is the lognormal distribution, i.e., a Gaussian distribution for the multiplying factor on the log scale. The lognormal distribution is symmetric in log scale, i.e., the chances of multiplying or dividing by a given factor are equal, which is not the case with a Gaussian model.

Figure 3.1 illustrates this property. It shows that a Gaussian approximation is relevant for a low mean calibration error (left), but inappropriate for larger multiplicative errors (right). In particular, in the latter case, a Gaussian distribution allows negative multiplicative factors, while they cannot be accommodated by the lognormal distribution.

<figure><img src="image_21.png" /><figcaption>Figure 3.1 – Quality of Gaussian approximation of a lognormal distribution for low and high variance. Left: the multiplicative noise corresponds to an error of 10% in average. The associated lognormal distribution is quite symmetric and is thus well approximated by a Gaussian distribution. Right: the multiplicative noise corresponds to an error of a factor 2 in average. The associated lognormal distribution is significantly asymmetric and is thus poorly approximated by a Gaussian distribution.</figcaption></figure>

**Student's t additive and Student's t multiplicative distributions** – Using a non-Gaussian uncertainty model usually leads to a complex likelihood function with no closed-form expression. For instance, in Kelly et al. (2012), the observation model is similar to Eq. 3.5, but with different noise models: the additive noise $\varepsilon^{(a)}$ and the log of the multiplicative noise $\log \varepsilon^{(m)}$ are assumed to be independent and to follow a Student's t distribution, whose pdf is

$$ \forall n, \ell, \quad \pi \left( \varepsilon_{n\ell}^{(a)} \right) = \frac{\Gamma\left(\frac{v+1}{2}\right)}{\sqrt{v\pi} \Gamma\left(\frac{v}{2}\right)} \left( 1 + \frac{\left(\varepsilon_{n\ell}^{(a)}\right)^2}{v} \right)^{-\frac{v+1}{2}}, \quad (3.6) $$

with $v$ the degree of freedom of the distribution, set to $v = 8$ in Kelly et al. (2012). The Student's t distribution boils down to the Cauchy distribution for $v = 1$ and to the Gaussian distribution for $v = \infty$. In the article, the covariance matrix $\Sigma^{(a)}$ of the additive noise is known and diagonal, and the calibration errors $\varepsilon_{n\ell}^{(m)}$ are assumed independent spectrally and with spatial correlation of exactly 1.

Setting a Student's t distribution on the additive noise $\varepsilon_{n\ell}^{(a)}$ and on the log of the multiplicative noise $\log \varepsilon_{n\ell}^{(m)}$ is equivalent to assuming a Gaussian distribution where the standard deviations are unknown but described with an inverse gamma prior. As the estimation of the standard deviations
is performed prior to the inverse problem, using a Student's t distribution is a statistically relevant way to encode uncertainty on the noise standard deviations.

With these two Student's t noise sources, the likelihood cannot be written in closed-form. The authors obtained a closed-form likelihood by using a hierarchical model, i.e., by including $\log \epsilon^{(m)}$ in the set of parameters to infer. This auxiliary variable was thus sampled, and ignored when deriving estimators.

**Other combinations of noises** – Mixtures of noise models (Robert and Casella, 2004, chapter 1) are sometimes exploited to set one model for standard observations and one for outliers. In astrochemistry, Gratier et al. (2016) uses a Gaussian mixture model (GMM) with two components

$$y_{n\ell}|\theta_n \sim (1-p_o)\mathcal{N}(\mu_s, \sigma_s^2) + p_o\mathcal{N}(\mu_o, \sigma_o^2), \quad p_o \in [0, 1], \qquad (3.7)$$

where $\mathcal{N}(\mu_s, \sigma_s^2)$ is the noise model for standard observations, $\mathcal{N}(\mu_o, \sigma_o^2)$ is a model for outliers, typically with $\sigma_o \gg \sigma_s$, and $p_o$ is the probability for the observation to be an outlier. Similarly, to evaluate distances to local molecular clouds, Zucker et al. (2019) mixes a Gaussian model and a uniform one. In both cases, the mixture parameter $p_o$ is also inferred, which makes the model hierarchical. Gratier et al. (2016) goes farther and also infers most of the noise model parameters from the data.

### 3.1.2.4 Including censorship in the likelihood function

Upper bounds on observations are quite common in ISM studies. The corresponding observables are sometimes discarded in the inversion process. Omitting observables leads to a loss of information that could damage the inference results. Including these upper bounds on observations in the likelihood function permits to account for all constraints provided by the observations. In statistics, censorship permits to include such upper limits $\omega \ge 0$. In case of Gaussian noise, the observation model becomes

$$y_{n\ell} = \max \{\omega, f_{\ell}(\theta_n) + \varepsilon_{n\ell}^{(a)}\} = \begin{cases} \omega & \text{if } f_{\ell}(\theta_n) + \varepsilon_{n\ell}^{(a)} \le \omega \\ f_{\ell}(\theta_n) + \varepsilon_{n\ell}^{(a)} & \text{otherwise} \end{cases} \qquad (3.8)$$

In (Ramambason et al., 2022), censorship is modeled with a half-normal distribution $\mathcal{N}_-(\omega, \sigma^2)$. Up to an additive constant, the associated likelihood function is

$$-\ln \pi(y_{n\ell}|\theta_n) = \frac{1}{2\sigma^2} (f_{\ell}(\theta_n) - \omega)^2 \mathbf{I}_{[-\infty, \omega]}(f_{\ell}(\theta_n)), \qquad (3.9)$$

where $\mathbf{I}_{[-\infty, \omega]}(f_{\ell}(\theta_n)) = 0$ if $f_{\ell}(\theta_n) \le \omega$ and $+\infty$ otherwise. This approach suffers from two main drawbacks. First, it prohibits values that are above the censoring threshold $\omega$, though values slightly above $\omega$ should not be unrealistic. Second, it biases $f_{\ell}(\theta_n)$ towards values close to the threshold $\omega$, while values of $f_{\ell}(\theta_n) \ll \omega$ are also compatible with the upper bound constraint.

In statistics, the standard approach to encode censorship in the likelihood is (Robert and Casella, 2004, chapter 1)

$$-\ln \pi(\omega|\theta_n) = -\ln \int_{-\infty}^{\omega} \pi(y_{n\ell}|\theta_n) dy_{n\ell} \qquad (3.10)$$

The integral covers all the integrated intensities $y_{n\ell} \le \omega$, and sums the associated uncensored likelihoods $\pi(y_{n\ell}|\theta_n)$.

Figure 3.2 illustrates this censorship modeling for a Gaussian additive noise model $\varepsilon_{n\ell}^{(a)} \sim \mathcal{N}(0, \sigma^2)$. In this case, the integral in Eq. 3.10 can be evaluated efficiently for Gaussian noise using the Gaussian cumulative density function (cdf) $\Phi(\omega|f_{\ell}(\theta), \sigma^2)$. Unlike the half Gaussian mentioned above, it does not forbid all physical parameters $\theta$ such that $f_{\ell}(\theta) > \omega$. Instead, it
<figure><img src="image_22.png" /><figcaption>Figure 3.2 – Illustration of censorship in likelihood for Gaussian additive noise models $\epsilon_{nl}^{(a)} \sim N(0, \sigma^2)$. The first four panels show four values of $f(\theta^{(i)}), i = 0, \dots, 3$ and the associated Gaussian additive pdf. The censored portion, below the censorship threshold $\omega$, is highlighted. The fifth panel shows the likelihood function taking censorship into account. The sixth panel shows the corresponding negative log-likelihood.</figcaption></figure>
smoothly and increasingly penalizes them as $f_{\ell}(\theta)$ gets farther to the threshold limit, i.e., as the
probability of $y_{nl} \le \omega$ becomes smaller.

This censorship modeling is already used in some ISM inversion codes (Blanc et al., 2015; Holdship et al., 2018; Thomas et al., 2018).

3.1.2.5 Avoiding the manual specification of a noise model

Finally, two types of approaches do not require specifying a noise model to perform the in-
ference. The first one abandons the noise model specification, while the second learns it from
observation without any signal.

**Likelihood-free approach** – The most famous approaches that avoid the noise specification are called *likelihood-free*, with for instance the approximate Bayes computing (ABC) methods, introduced in Beaumont et al. (2002) for statistical inference in genetics. See, e.g., Cranmer et al. (2020) for a review on likelihood-free methods and Beaumont (2010) for a review on ABC methods. Put in simple words, ABC methods replace the likelihood function by a distance between summary statistics of true and simulated observations. They were first used in cosmology in Cameron and Pettitt (2012). Then, they became quite popular, e.g., with Weyant et al. (2013), COSMOABC (Ishida et al., 2015), ABCPMC (Akeret et al., 2015)<fnref n="1" />, and ASTROABC (Jennings and Madigan, 2017). ABC sampling is much rarer in ISM studies. It is used e.g., in Robin et al. (2014) to study the thick disk in the Milky Way.

**Learning the noise from observations** – Abandoning all hope of specifying an accurate noise model and replacing it with a distance on some statistics seems suboptimal. A recent alternative consists in learning a correct noise model from observations of noise realizations without any signal. For instance, in Hahn et al. (2019), the authors combine GMMs and independent component analysis (ICA) (see, e.g., Hastie et al. (n.d., chapter 14)) to better approximate the likelihood from datasets usually used for covariance analysis. In Legin et al. (2023), the noise distribution of some observations of the Hubble space telescope and of the James Webb spatial telescope (JWST) are learned with a score-based diffusion generative model (see, e.g., Song and Ermon (2020)). This type of approach fully learns the noise distribution. Hierarchical models can also learn from observations to estimate noise parameters. However, these parameters depend on a noise model set a priori, which is more restrictive. To the best of our knowledge, no published work in the ISM community tries to fully learn a noise distribution from observations.

Our proposal

In the general inverse problem considered in this thesis, we consider a Gaussian additive noise
$\varepsilon^{(a)} \sim \mathcal{N}(0, \Sigma^{(a)})$, a lognormal multiplicative error $\varepsilon^{(m)} \sim \mathcal{N}(0, \Sigma^{(m)})$ and a censorship level
$\omega \ge 0$. The multiplicative noise includes calibration error and, depending on the observation,
a model misspecification error. The two covariance matrices $\Sigma^{(a)}$ and $\Sigma^{(m)}$ are assumed to be
known and diagonal. Censorship is modeled as in Eq. 3.10. As the multiplicative noise is not
Gaussian, the overall likelihood function has no closed-form expression. In Section 3.4, we detail
the uncertainty model, its origin and the corresponding assumptions. Chapter 5 (Section 5.1.1)
proposes an approximation with controlled error of the associated intractable likelihood function.

3.1.3 Prior distributions and regularization functions

A prior distribution, or equivalently a regularization function, can encode any physical knowl-
edge on the parameters to infer Θ. There are two main philosophies regarding inference. The
first considers that the statistical model should only exploit information contained in the observed
data, and thus sets the prior as non-informative as possible, if any. The second prefers to add

<footnote n="1">1. https://abcpmc.readthedocs.io/en/latest/</footnote>
physical prior knowledge, e.g., to favor desirable properties in the posterior or to rule out non-
physical solutions.

**Non-informative and weakly-informative priors** are extremely widespread in ISM studies. The most common is a uniform prior on $\Theta$, with validity intervals on its components. This prior is either stated explicitly (e.g., in Behrens et al. (2022); Blanc et al. (2015); Thomas et al. (2018); Holdship et al. (2018)) or implicitly, when working on a lattice dataset (e.g., in Joblin et al. (2018); Sheffer et al. (2011); Sheffer and Wolfire (2013)). This uniform prior can be set on the linear or log scale, depending on the physical parameter and its dynamic range. For this prior to be weakly informative, the lower and upper bounds are set to obtain wide intervals. These bounds typically correspond to the orders of magnitude spanned by the physical parameter – see, e.g., Wu et al. (2018); Joblin et al. (2018).

**Informative priors** are punctually exploited to encode additional knowledge. There are two main classes of informative priors in the considered astrophysics communities:

- The first class exploits the structure of the physical parameter $\Theta$. For instance, Paumard et al. (2014); Ciurlo et al. (2016); Paumard et al. (2022) fit Gaussian line profiles on hyperspectral observations. Exploiting the map structure of the physical parameter, they enforce spatial smoothness in the maps by adding an $L_1-L_2$ spatial regularization. The ROHSA fitting algorithm (Marchal et al., 2019) also fits Gaussian line profiles on hyperspectral maps using a spatial regularization. However, ROHSA relies on the $L_2$ norm of the map Laplacian as a regularization. Spatial regularization improves reconstructions particularly in the low signal-to-noise ratio regions, where the observations mostly contain noise.

- The second class directly encodes physical information. For instance, Wu et al. (2018) rules out non-physical solutions with a physically informed prior on predicted observations $f(\theta)$. This prior constrained the reconstructed integrated intensities to be decreasing for the CO lines $J = 11-10$, $J = 12-11$, and $J = 13-12$, as this condition seems always satisfied in observations. In stellar properties studies, more advanced physics-informed priors are sometimes used, such as the Hertzsprung-Russell Diagram prior on population of stars in Bailer-Jones (2011). This prior ensures that inference results are physically consistent.

In the case of informative priors, one may know a relevant family of distributions for a con-
sidered inverse problem, but face uncertainty on the values of its parameters, e.g., its mean or
covariance. This hyperparameter tuning problem is crucial for inference, as different values yield
different priors and thus different trade-offs. For instance, Ciurlo et al. (2016) presents a tedious
manual setting of the six hyperparameters of its $L_1-L_2$ spatial regularization. Hierarchical Bayesian
models, also called multi-level models, represent an intermediate between weakly-informative pri-
ors and wrong over-informative priors. They offer a way to decompose the prior knowledge into
several elementary components, while accounting for uncertainties over each included item of
knowledge. The parameters of the prior distribution can then be inferred from the data along
with the physical parameters of interest $\Theta$.

In the ISM community, Kelly et al. (2012) popularized hierarchical models for dust studies. For instance, Galliano (2018) – which introduces the dust Bayesian inversion code HERBIE – and Galliano et al. (2021) rely on a multivariate Student’s t distribution prior on all physical parameters. Using a hierarchical statistical model, the mean and covariance parameters of this prior distribution are inferred along with the physical parameters. To simplify the inference of the covariance matrix, they exploit the separation method (Barnard et al., 2000) to divide covariances into correlations and standard deviations. Hierarchical models have been also used in other dust studies e.g., in Juvela et al. (2013); Veneziani et al. (2013). To the best of our knowledge, this approach still remains to be applied in other ISM fields.

**Our proposal** – The proposed approach will exploit the spatial regularity of the parameters
to be inferred through the $L_2$ norm of the Laplacian operator, while restraining their values to
physically consistent ranges.

## 3.2 Statistical inference

In this section, we review the inference methods widespread in ISM studies to evaluate estimators defined either as minimum of a loss function or with an integral.

### 3.2.1 Optimization-based inference

The MLE and MAP are widespread estimators in ISM studies. The methods listed below are introduced in Chapter 2 (Section 2.2.1).

**Grid search** – Searching for the point in the dataset that best reproduces the observations is very common in the inversion of integrated intensities of ionic, atomic or molecular lines e.g., in Sheffer et al. (2011); Sheffer and Wolfire (2013); Joblin et al. (2018); Lee et al. (2019). This grid search approach only works in low dimensional settings, as the number of required precomputed models grows exponentially with the number of physical parameters. Besides, it requires very fine grids to return accurate estimations, which means many evaluations of the forward model.

**Gradient descent (GD)** – GD algorithms are also quite common in ISM studies. They lead to more accurate results and work in higher dimensional settings. One of the most common GD method is the Levenberg-Marquardt algorithm, used e.g., in Schilke et al. (2010) to fit Herschel/HIFI spectra of an ortho-H₂O⁺ line with a radiative transfer code, and in Galliano et al. (2003) in a dust study. Other GD algorithms are sometimes used in ISM studies, such as conjugate gradient descent in Paumard et al. (2022) and preconditioned gradient descent – using the limited memory BFGS (L-BFGS) preconditioner – in Wu et al. (2018). These gradient descent algorithms perform well with unimodal posteriors. However, unlike grid search methods, they can get trapped in local minima in case of a multimodal posterior.

**Meta-heuristics** – In multimodal cases, meta-heuristics such as simulated annealing and genetic algorithms are sometimes preferred to escape from local minima, e.g., in the MAGIX inversion code (Möller et al., 2013). These algorithms are also used in cosmology (Hannestad, 1999) or in the inference of star properties (Sarro et al., 2018).

The main limitation of optimization-based methods is that they only provide point estimates, with no information about the associated uncertainties. This is highly problematic, in particular with non-informative or weakly-informative priors and ill-posed problems, i.e., badly constrained problems. For instance, low SNR observations or tracers with limited sensitivity with respect to physical parameters can lead to non-physical results. Additional studies are necessary to quantify uncertainty. For instance, Roueff et al. (2021) quantifies the uncertainty associated with the MLE with the Cramér-Rao bound. As noted by Panter et al. (2003), this approach is only relevant when the posterior is well approximated by a Gaussian at its mode, which is generally not the case in astrophysics. In contrast, sampling-based approaches provide credibility intervals on the physical parameter $\Theta$ along with point estimates for general posterior distributions.

### 3.2.2 Sampling-based inference

A sampling-based approach gives access to estimators defined with an integral – see Chapter 2 (Section 2.1.2). Several methods were already used in the ISM community.

**Riemann integration with a grid of pdf evaluations** – A simple method to evaluate an integral performs a Riemann integration with a grid in the physical parameter space. The nega-
tive log-posterior is evaluated at each point of the grid. This method was applied in some dust physics studies (Da Cunha et al., 2008; Pacifici et al., 2012). The Bayesian-like HII-CHI-MISTRY code (Pérez-Montero, 2014) uses a grid of CLOUDY models to study HII regions. From a grid of models, it returns a mean estimator computed with a weighted $\chi^2$. This grid approach was implemented for other ISM inverse problems with extremely low-dimensional settings ($N = 1, D \le 10$), e.g., in the codes IZI (Blanc et al., 2015), BOND (Vale Asari et al., 2016), NEBULABAYES (Thomas et al., 2018) or in Villa-Vélez et al. (2021). As the required number of precomputed models grows exponentially with the dimension, this approach does not scale to higher dimensions.

**Markov chain Monte Carlo (MCMC) algorithms**, introduced in Chapter 2 (Section 2.2.2.2), scale better to high dimensions. MCMC algorithms were popularized in astronomy through cosmology in Christensen et al. (2001). The first public MCMC code, COSMOMC was published in Lewis and Bridle (2002). These articles used random walk Metropolis-Hastings (RWMH) to generate posterior samples, arguably the most widespread sampling kernel. COSMOMC then became common in cosmological parameter inference, see, e.g., Tegmark et al. (2004), and was also applied in stellar formation (Acquaviva et al., 2011), dust studies in nearby galaxies (Chevallard et al., 2013) and fundamental properties of galaxies (Serra et al., 2011). All these applications had low dimensionality – up to 11 in Lewis and Bridle (2002). RWMH was also applied in ISM studies e.g., in dark interstellar clouds (Makrymallis and Viti, 2014), in dust (Paradis et al., 2010) or in astrochemistry (Makrymallis and Viti, 2014).

Some more advanced Bayesian methods, introduced in Chapter 2 (Section 2.A), are also already popular in the ISM community:

- The affine-invariant MCMC sampler and the associated `EMCEE` package is sometimes considered as the most popular MCMC algorithm in astronomy (Thrane and Talbot, 2019). It was applied in astrochemistry (Gratier et al., 2016; Holdship et al., 2018; Keil et al., 2022), extragalactic molecular gas (Yang et al., 2017). It is also popular out of the ISM community, e.g., in stellar physics (Jackiewicz, 2020; Kashyap et al., 2021). This sampler can exclusively address low dimensional problems, and requires to be initialized close to a mode (Foreman-Mackey et al., 2013).

- Nested sampling is also a very popular Bayesian framework in astrophysics and ISM studies. Popular algorithms include `MULTINEST` (Feroz et al., 2009), `DYNESTY` (Speagle, 2020) and `ULTRANEST` (Buchner, 2016b; Buchner, 2021), that were developed by and for astrophysicists. `MULTINEST` is applied e.g., to study extragalactic molecular clouds (Kamenetzky et al., 2014; Chevallard and Charlot, 2016) and cosmic ray propagation (Jóhannesson et al., 2016). `DYNESTY` is applied e.g., to evaluate distances between the Earth and nearby molecular clouds (Zucker et al., 2018; Zucker et al., 2019) or to reconstruct their 3D structure (Zucker et al., 2021). `ULTRANEST` is applied e.g., on extragalactic SEDs (Behrens et al., 2022). In particular, these codes can handle multimodal distribution, but remain limited to low dimensional distributions, e.g., up to 30 in a toy case for `MULTINEST` Feroz et al. (2009) and in the inverse problem considered in Jóhannesson et al. (2016).

- Sequential MC (SMC) is quite rare in ISM but is nonetheless applied in Ramambason et al. (2022) in an inversion of integrated intensities of ionic, atomic and molecular emission lines – the problem we are interested in in this thesis. In that work, it is chosen for its ability to sample from multimodal distributions in low dimension settings – up to 28 in the paper, for three sectors – and to compute the Bayesian evidence. SMC is also used when the likelihood is not specified with an ABC approach, especially in cosmology (see, e.g., Cameron and Pettitt (2012); Weyant et al. (2013); Jennings and Madigan (2017)).

- One sampling algorithm based on a combination of MH, Gibbs sampling and ancillarity-sufficiency interweaving strategy (ASIS) (Yu and Meng, 2011) is widespread in dust studies
that involve a hierarchical model. ASIS is a sampling strategy that relies on complementary data augmentation schemes to better explore the posterior, including along strongly correlated directions. This sampler was first used in a dust study in Kelly et al. (2012), with implementation details provided in Kelly (2011). It was then coupled with a Gibbs block coordinate approach to exploit the structure in the physical parameter $\Theta$ in very high dimensional applications. For instance, Galliano (2018) involve a few thousand dimensions, and (Galliano et al., 2021) about 10 000.

**Our proposal** – Recall that in the inverse problem considered in this thesis, the posterior distribution is non-log-concave and thus potentially multimodal. The negative log-posterior is also assumed to be twice differentiable and non gradient Lipschitz continuous. Besides, the physical parameters we aim at inferring live in high dimensions – up to $O(10^4)$. The RWMH or the affine-invariant sampler do not scale well with dimensionality and thus are not able to sample from the considered posterior. Similarly, nested sampling and sequential MC could bypass the multimodality and non gradient Lipschitz issues, but do not scale well to high dimensions. Finally, ASIS is based on two complementary data augmentation schemes. There exists other data augmentation schemes in the signal processing community dedicated to accelerating MCMC algorithms, e.g., (Vono et al., 2019). In this thesis, we do not resort to data augmentation strategies. In Chapter 5 (Section 5.2), we define a dedicated sampler. Section 5.4 demonstrates that it yields state-of-the art performance on two multimodal distributions. Future work may include a comparison of the proposed sampler with MCMC algorithms involving data augmentation.

## 3.3 Comparing and checking observation models

As discussed in Chapter 2 (Section 2.3), once the inference is performed, one may want to evaluate the quality of the fit and compare it to other fits. Model selection is a quite common approach in astrophysics. However, it only compares models, as the values of the associated criteria are not interpretable. Conversely, model checking assesses a model individually, and yields interpretable diagnoses.

In this section, we first list applications of model comparison with Bayesian evidence and information criteria. Then, applications of model checking methods are reviewed.

### 3.3.1 Model comparison

Model comparison is an active research topic in astrophysics and cosmology that motivated the development of multiple statistical methods. As discussed in Chapter 2, model comparison can be performed by estimating the expected log-predictive density (elpd) or the marginal likelihood, also called Bayesian evidence.

The **elpd and information criteria** are rarely used in ISM studies, except for some simple methods. For instance, in a study of star formation history in other galaxies, Acquaviva et al. (2011) compares models $M_i$ using the $\chi^2$ values of their respective best-fit. As it uses an additive uncorrelated Gaussian noise model for a point estimate $\hat{\Theta}$, this criterion is equivalent to the log-predictive density (lpd), introduced in Chapter 2 (Section 2.3.1). The extragalactic SED analysis from Villa-Vélez et al. (2021) and the study of the thick disc in the Milky Way in Robin et al. (2014) both use Bayesian approaches and approximate the full posterior distribution. Both evaluate their results with the Bayesian information criterion (BIC) (Section 2.3.1) using the maximum a posteriori $\hat{\Theta}_{MAP}$. The ABC-based approach from Robin et al. (2014) resorts to a surrogate likelihood function to evaluate the BIC. Finally, Lebouteiller and Ramambason (2022) evaluates the widely applicable information criterion (WAIC) and leave-one-out cross-validation approximation from Vehtari et al. (2017), in addition to the Bayesian evidence.
The Bayesian evidence is a more widespread model comparison criterion in astrophysics. The need for methods that address multimodal distributions in high dimension led to many dedicated statistical developments. For instance, new statistical methods were proposed by astrophysicists to evaluate the marginal likelihood from posterior samples (McEwen et al., 2022) or with nested sampling, with e.g., development of MULTINEST (Feroz and Hobson, 2008), ULTRANEST (Buchner, 2016b; Buchner, 2021) and DYNESTY (Speagle, 2020). Evaluating the Bayesian evidence quickly became popular in cosmology to compare models – see Trotta (2008) for a review. It is also widespread in astrophysical fields such as stellar physics, e.g., in Hatta et al. (2022). In ISM studies, the evidence is sometimes evaluated:

- Ramambason et al. (2022) computes the evidence directly with the SMC sampler in a study of extragalactic HII regions.

- The inversion code BEAGLE (Chevallard and Charlot, 2016), which analyzes the SEDs of galaxies, evaluates the evidence using MULTINEST.

- Zucker et al. (2021) evaluates the evidence using DYNESTY to study the 3D structure of nearby molecular clouds.

- Kamenetzky et al. (2014), which studies molecular clouds properties in nearby galaxies, resorts to model selection to compare local modes of one posterior distribution, using “local evidences”. This evaluation of local evidence of each mode is performed with MULTINEST.

### 3.3.2 Model checking with posterior predictive assessment

In ISM studies, the obtained loss function value in a non-linear least squares problem, often noted $\chi^2$, is used in a three-case interpretation (Ivezić et al., 2020, chapter 4). When $\chi^2 > 1$, the estimated parameters are judged not able to reproduce observations. A $\chi^2 \ll 1$ suggests an overfit or an overestimation of uncertainties on observations. Finally, $\chi^2 \simeq 1$ is the ideal case, indicating physical parameters that reasonably reproduce the observations. This interpretation was applied in ISM studies e.g., in Chevance et al. (2016); Joblin et al. (2018); Villa-Vélez et al. (2021).

This interpretation comes with some limitations and was already criticized in the astrophysics community (Andrae et al., 2010). In particular, the degree of freedom is challenging to estimate for non-linear forward models and when a prior distribution is considered. Besides, this $\chi^2$ rule only applies to Gaussian noise models. Finally, when $\Theta$ is described by a posterior distribution and not by a point estimate $\hat{\Theta}$, the obtained $\chi^2$ value is approximated with a Monte Carlo (MC) estimator. Replacing the true $\chi^2$ with an Monte Carlo (MC) estimator introduces an error that is seldom accounted for in astrophysics.

In Chapter 2 (Section 2.3), we introduced the so-called Bayesian p-value (Eq. 2.57). This p-value permits to check the ability of the forward model to reproduce observations. The aforementioned $\chi^2$ rule can be transformed into a Bayesian p-value approach for a point estimate $\hat{\Theta}$ using the $L_2$ test statistic from Eq. 2.55, reminded below

$$T(\tilde{\mathbf{Y}}, \mathbf{\Theta}) = \sum_{n=1}^{N} \sum_{\ell=1}^{L} \frac{(y_{n\ell} - f_{\ell}(\mathbf{\Theta}_n))^2}{\sigma_{n\ell}^2}. \qquad (3.11)$$

The inversion code BEAGLE (Chevallard and Charlot, 2016) aims at modeling and interpreting spectral energy densities (SEDs) of galaxies. In BEAGLE, $y_\ell$ is the observed flux in band $\ell$, $f_\ell(\Theta)$ the flux predicted by the model in band $\ell$, and $\sigma_\ell$ the standard deviation for band $\ell$. This code uses the $\chi^2$ both as a likelihood function and as a test statistic T for the Bayesian p-value. As the authors consider a posterior distribution and not a point estimate $\hat{\Theta}$, this p-value is intractable. The authors approximate it using the MC estimator presented in Gelman et al. (1996) and in Eq. 2.59.
Galliano et al. (2021) infers dust properties from the SED of nearby galaxies using the HER-BIE inversion code (Galliano, 2018). In this article, $y_\ell$ is the observed flux in band $\ell$, $f_\ell(\Theta)$ the flux predicted by the model in band $\ell$, and $\sigma_\ell$ the standard deviation for band $\ell$. Lebouteiller and Ramambason (2022) fits physical parameters to integrated intensity of atomic and molecular emission lines using CLOUDY. In this work, $y_\ell$ and $f_\ell(\Theta)$ are integrated intensities of an atomic or molecular emission line $\ell$, and $\sigma_\ell$ the corresponding observation standard deviation. In both these articles, the authors use a likelihood more complex than uncorrelated Gaussian noise. For instance, both include censored observations, i.e., upper limits on some low intensity observations. With a censored likelihood, the posterior enforces observation reproductions $\tilde{y}_\ell$ to be below an upper limit. Both these articles also use the $\chi^2$ as a test statistic $T$. With this statistic $T$, reproductions that are far below this upper limit are heavily penalized by the test statistic, while relevant for the likelihood function. A test statistic that favors the same behaviors as the likelihood should be preferred. Finally, both high and low $p$-values are rejected in Galliano et al. (2021). While rejection for high $p$-values is indicated in Gelman et al. (2015) for a simple test statistic, $p$-values close to 1 do not indicate a problematic fit with the $\chi^2$ statistic. On the contrary, with this discrepancy measure $T$, $p$-values close to 1 indicate that the model can reproduce observations with smaller errors than the uncertainties in the noise model.

**Our proposal** – In Chapter 5 (Section 5.3), we resort to hypothesis testing to assess the compatibility of a forward model with the observations. We evaluate the Bayesian $p$-value from Gelman et al. (1996), with three specificities:

1. We do not set the test statistic to the $\chi^2$. To ensure that the $p$-value is consistent with the likelihood function in general cases, we set the test statistic to the negative log-likelihood $T(Y, \Theta) = -\ln \pi(Y|\Theta)$. This choice generalizes the case of BEAGLE to arbitrary observation models.

2. We compute maps of $p$-values instead of one global $p$-value. This approach helps to identify regions that are badly modeled by the Meudon PDR code in analyses.

3. We account for uncertainties on the $p$-value that come from the MC evaluation. This avoids the rejections that are due to chance and caused by insufficient number of samples.

## 3.4 Our observation model

In this thesis, we consider hyperspectral observations such as the large maps IRAM-30m "ORION-B" Large Program. Such maps are usually transformed into maps of integrated intensities of $L \sim 5-30$ ionic, atomic or molecular emission lines. In this section, we describe the reduction process considered in this thesis and the associated uncertainty model.

### 3.4.1 Noise model on the original hyperspectral cube

We start with an observed hyperspectral cube $O = (o_{nk}) \in \mathbb{R}^{N \times K}$, where $N$ is the number of pixels (observed spatial positions) and $K$ the number of spectral channels (observed frequencies). For instance, the Orion-B observation described in Pety et al. (2017) contains $N = 10^6$ pixels and $K = 200\,000$ spectral channels. Usually, observation maps contain between $N = O(1)$ and $O(10^4)$ pixels.

Figure 3.3 shows examples of spectra of molecular emission lines from the Orion-B observations presented in Pety et al. (2017). Each emission line is limited to an effective small range of spectral channels.

As in Kelly et al. (2012); Galliano (2018); Galliano et al. (2021), we assume a mixture of multiplicative and additive noises on the voxels, i.e., the elements of the cube:

$$ \forall n \in [1, N], \quad \forall k \in [1, K], \quad o_{nk} = \tilde{\varepsilon}_{nk}^{(c)} I_{nk} + \tilde{\varepsilon}_{nk}^{(a)}, \qquad (3.12) $$
with *n* a pixel index, *k* a spectral channel, and *I*<sub>nk</sub> the unaltered signal of interest. The term

$\tilde{\varepsilon}^{(c)} = (\tilde{\varepsilon}_{nk}^{(c)})_{nk}$ is a multiplicative noise due to calibration error, and $\tilde{\varepsilon}^{(a)} = (\tilde{\varepsilon}_{nk}^{(a)})_{nk}$ an additive

measurement noise.

**The measurement noise** $\tilde{\varepsilon}^{(a)}$ combines two sources of uncertainty. A common observation procedure defines $o_{nk}$ as the difference between two measurements, one “off” (observation of a dark portion of the sky), and one “on” (observation of the position of interest $n$). Each measurement counts photons on a small area during a fixed duration. As photons are considered independent, the number of observed photons is modeled as independent realizations of a Poisson distribution. When many photons are observed, the Poisson distribution can be well approximated by a Gaussian distribution with known variance. Finally, the difference of the two measurements yields a difference of two independent Gaussian distributions, which is itself Gaussian. Besides, thermal noise also affects the observations. This thermal noise can be accurately modeled by an additive Gaussian noise as well.

As a sum of these two Gaussian noises, $\tilde{\varepsilon}_{nk}^{(a)}$ is modeled by a Gaussian distribution. In this thesis, its covariance matrix $\tilde{\Sigma}^{(a)}$ is assumed diagonal, i.e., $\tilde{\Sigma}^{(a)} = \text{diag}(\tilde{\sigma}_{a,n\ell}^2)$, for simplicity. Realizations on different voxels $(n, k)$ are thus assumed independent. This assumption is generally not verified on real observations. For instance, Einig et al. (2023) the structure of the noise covariance in the IRAM-30m Orion-B observations (Pety et al., 2017).

**The calibration noise** $\tilde{\varepsilon}^{(c)}$ represents a relative error with respect to the true signal $I_{nk}$ of 10% or lower (Einig et al., 2023). As described in Chapter 3 (Section 3.1.2.3), an additive Gaussian noise cannot capture accurately such relative uncertainty, contrary to a multiplicative noise. As in Kelly et al. (2012), we assume a lognormal model, i.e., a Gaussian distribution on $\ln \tilde{\varepsilon}^{(c)}$. The spectral correlation is partly taken into account. We assume that the multiplicative noise realizations $\tilde{\varepsilon}_{nk}^{(c)}$ in all the spectral channels associated with one emission line $\ell$ are equal. The noise realization value is denoted $\varepsilon_{n\ell}^{(c)}$. Conversely, we neglect the correlation between the spectral channels of distinct lines. As in the additive noise, the noise spatial correlation is neglected. Overall, for a line $\ell$ and a pixel $n$, we assume $\varepsilon_{n\ell}^{(c)} \sim \text{Lognormal}(-\sigma_c^2/2, \sigma_c^2)$, where the non-zero mean on the log scale ensures $\mathbb{E}[\varepsilon_{n\ell}^{(c)}] = 1$.

**Intensities** $I_{nk}$ can cover multiple decades within one hyperspectral cube – see, e.g., Pety et al. (2017). The dominant noise for a voxel $o_{nk}$ thus depends on the intensity $I_{nk}$. Indeed, if $I_{nk} \gg \tilde{\sigma}_{a,n\ell}$, the additive noise becomes negligible and the multiplicative dominates. Conversely, if $I_{nk} \ll \tilde{\sigma}_{a,n\ell}$, the additive noise dominates. Therefore, one cannot neglect one of the two sources of uncertainty.
<figure><img src="image_23.png" /><figcaption>Figure 3.3 – Examples of molecular emission line spectra for different molecules in the Orion-B cloud. The line profiles displayed here are averaged over all pixels in the observed map. The velocities in the horizontal axes are translated for each line $ℓ$ such that the 0 corresponds to its characteristic frequency $ν_ℓ$. Adapted from Pety et al. (2017).</figcaption></figure>

### Details for statisticians 3.1: Units for specific and integrated intensities

In Figure 3.3, the spectral channels are called velocity channels and expressed in km s⁻¹, and the intensities are expressed as temperatures. Integrated intensities are thus expressed in K km s⁻¹. These conventions are widespread in radio observations. In physics, the integrated intensities are expressed in erg cm⁻² s⁻¹ sr⁻¹, using the centimetre-gram-second (cgs) unit system. The conversion is detailed here.

The Doppler effect provides a simple bijective relation between frequency $ν$ and velocity $v$

$$d\nu \, [\text{Hz}] = \frac{\nu}{c} d\nu \, [\text{m s}^{-1}] = \left(10^3 \frac{\nu}{c}\right) d\nu \, [\text{km s}^{-1}]. \quad (3.13)$$

The measured specific intensity $I_ν$ in a spectral channel of frequency $ν$ is converted to a temperature $T$ using a black body model and the Rayleigh-Jeans (RJ) approximation:

$$I_{\nu} \, [\text{W m}^{-2} \text{ Hz}^{-1} \text{ sr}^{-1}] = \frac{2h\nu^3}{c^2} \frac{1}{\exp\left\{\frac{h\nu}{k_B T [\text{K}]\right\} - 1} \stackrel{\text{RJ}}{\approx} \left(\frac{2\nu^2}{c^2} k_B\right) T [\text{K}], \quad (3.14)$$

with $k_B$ the Boltzmann constant and $h$ the Planck constant. Combining Eq. 3.13 and Eq. 3.14, the integrated intensity $I_ℓ$ reads

$$I_{\ell} \, [\text{erg cm}^{-2} \text{ s}^{-1} \text{ sr}^{-1}] = 10^3 I_{\ell} \, [\text{W m}^{-2} \text{ sr}^{-1}] \quad (3.15)$$

$$= 10^3 \int I_{\ell} \, [\text{W m}^{-2} \text{ Hz}^{-1} \text{ sr}^{-1}] \times d\nu \, [\text{Hz}] \quad (3.16)$$

$$= \left( 2 \times 10^6 \frac{\nu^3}{c^3} k_B \right) \int T [\text{K}] \times d\nu \, [\text{km s}^{-1}]. \quad (3.17)$$
### 3.4.2 From the hyperspectral cube to maps of integrated intensities

There are multiple methods to evaluate the integrated intensity of an emission line $ℓ$ at a given pixel $n$ from the hyperspectral cube. A widespread approach is to assume a line profile shape, usually Voigt or Gaussian, and to fit it to the observations.

#### Details for statisticians 3.2: Evaluating integrated intensity with a line profile

For a static cloud, an emission line $ℓ$ that corresponds to a quantum de-excitation is characterized by a central frequency $ν_ℓ$ and a line profile $φ$ such that $∫ φ(ν)dν = 1$. The integrated intensity $I_ℓ$ associated with an emission line $ℓ$ is therefore obtained with

$$I_ℓ = \int I_νφ(ν - ν_ℓ)dν. \qquad (3.18)$$

In astrophysics, the line profile $φ$ is usually a Voigt profile, i.e., the convolution of a Gaussian profile and of a Lorentzian profile. The Gaussian profile accounts for the speed distribution of the gas particles and for microturbulence. The Lorentzian profile, also called Cauchy profile, models a quantum effect that causes uncertainty on the energy level values. The Voigt profile accounts for the three aforementioned sources of uncertainty. It is commonly used for absorption observations, e.g., for H or the H$_2$ absorption lines in the UV domain. For emission lines, it is often simplified to a Gaussian profile, as the Lorentzian wings get drowned in thermal noise. See Draine (2011, chapter 6) for more information on standard line profiles.

In real observations, a cloud is generally not static, which causes perturbations in the characteristic frequency $ν_ℓ$. Codes such as CUBEFIT (Paumard et al., 2022) perform a Gaussian profile fit to each spectrum to account for this perturbation. Other codes such as ROHSA (Paumard et al., 2022) perform a Gaussian mixture fit to account for the existence of multiple components in the cloud, each with a specific speed.

Figure 3.3 shows spectra of line emissions. In the high SNR regime, fitting a profile can provide an accurate fit and yield an accurate estimate of the integrated intensity $y_{nℓ}$. For instance, fitting a Gaussian profile seems relevant HCO$^+$ $J = 1-0$ and C$^{18}$O $J = 1-0$ emission lines (top right). Similarly, fitting a Gaussian mixture on the $^{12}$CO $J = 1-0$ profile (top left) seems relevant. However, it is unclear how to define an uncertainty model on the resulting integrated intensities. Besides, in the low SNR regime (bottom left) or in the absence of signal (bottom right), fitting a Gaussian profile may either provide a non-physical result or fail. In the latter case, the pixel is usually considered censored, but it is unclear with which upper bound – see Chapter 3 (Section 3.1.2.4) for an introduction to likelihood model in presence of observation censoring. See e.g., Gaudel et al. (2023) for another application to the Orion-B cloud to study the gas kinematics.

A second approach consists in summing the spectral channels associated with a given line $ℓ$. Some approaches such as Einig et al. (2023) may denoise the spectrum beforehand. Such denoising relies on assumptions on the spectrum regularity (e.g., continuity or differentiability) and on the noise (e.g., zero mean), but does not assume any shape for the profile. We do not resort to denoising to avoid distorting the signal of interest and to exploit a relevant noise model. For a given emission line $ℓ$ at pixel $n$, the integrated intensity $y_{nℓ}$ is defined as the sum of a set of spectral channels $[k_{min}^{(ℓ)}, k_{max}^{(ℓ)}]$:

$$y_{nℓ} = \sum_{k=k_{min}^{(ℓ)}}^{k_{max}^{(ℓ)}} o_{nk} = \sum_{k=k_{min}^{(ℓ)}}^{k_{max}^{(ℓ)}} (\varepsilon_{nℓ}^{(c)} I_{nk} + \tilde{\varepsilon}_{nk}^{(a)}) = \varepsilon_{nℓ}^{(c)} \sum_{k=k_{min}^{(ℓ)}}^{k_{max}^{(ℓ)}} I_{nk} + \sum_{k=k_{min}^{(ℓ)}}^{k_{max}^{(ℓ)}} \tilde{\varepsilon}_{nk}^{(a)}. \qquad (3.19)$$
The first sum defines the integrated intensity $I_{nk}$ before any alteration by noise. The second sum is a sum of independent Gaussian distributions with known standard deviations $\tilde{\sigma}_{a,nk}$, and is therefore itself a Gaussian distribution $\varepsilon_{n\ell}^{(a)} \sim N(0, \sigma_{a,n\ell}^2)$, with $\sigma_{a,n\ell}^2 = \sum_{k=k_{\min}^{(\ell)}}^{k_{\max}^{(\ell)}} \tilde{\sigma}_{a,nk}^2$. We rewrite Eq. 3.19 as

$$y_{n\ell} = \varepsilon_{n\ell}^{(c)} I_{n\ell} + \varepsilon_{n\ell}^{(a)}. \qquad (3.20)$$

The set of spectral channels $[k_{\min}^{(\ell)}, k_{\max}^{(\ell)}]$ may be fixed for all pixels $n$ or adaptive. An adaptive set may result in reduced variance on additive noise, but include a bias on the integrated intensities (Pety et al., 2017). In particular, an adaptive approach may cut the wings of the line profile and underestimate the integrated intensity $I_{n\ell}$. We use a fixed interval $[k_{\min}^{(\ell)}, k_{\max}^{(\ell)}]$ for all pixels $n$, at the cost of potentially higher noise variance. The intervals were set conservatively, reasonably large to include the whole line profile, even in case of Doppler shifting.

### 3.4.3 Forward model, observational effects and censorship

**The Meudon PDR code as a forward model and model misspecification** – A usual hypothesis in inversion procedures is that the forward model $\mathbf{f}$ can predict exactly this integrated intensity, i.e., that there exists a physical parameter $\theta_n$ such that $I_{n\ell} = f_{\ell}(\theta_n)$. In this thesis, the forward model is set to the Meudon PDR code, introduced in Chapter 1 (Section 1.3). As already mentioned, the Meudon PDR code might not be able to exactly reconstruct the observations. For instance, the reaction constants, collision rates and photo-reaction cross-sections used in the code come from laboratory measurements or from theoretical computations. The uncertainty on these values induces uncertainty on the code predictions. Besides, the encoded physics are based on physical assumptions and simplifications, such as the parallel slab geometry. Finally, in inference procedures the Meudon PDR code $\mathbf{f}$ is replaced by the neural network emulator $\tilde{\mathbf{f}}$ proposed in Chapter 4, which introduces another approximation error.

We introduce a second multiplicative error source $\varepsilon_{n\ell}^{(v)}$ to encode uncertainty on the model validity such that $I_{n\ell} = \varepsilon_{n\ell}^{(v)} \tilde{f}_{\ell}(\theta_n)$. For simplicity, we assume $\varepsilon_{n\ell}^{(v)} \sim \text{Lognormal}(-\sigma_{\text{mod}}^2/2, \sigma_{\text{mod}}^2)$ with independent and identically distributed (i.i.d.) realizations for different pixels $n$ or lines $\ell$. In this case, the full multiplicative uncertainty is $\varepsilon_{n\ell}^{(m)} = \varepsilon_{n\ell}^{(v)} \varepsilon_{n\ell}^{(c)}$. As a product of log-normal distributions, this multiplicative uncertainty is itself a lognormal distribution $\varepsilon_{n\ell}^{(m)} \sim \text{Lognormal}(-\sigma_m^2/2, \sigma_m^2)$, with $\sigma_m^2 = \sigma_c^2 + \sigma_{\text{mod}}^2$. The observation model in Eq. 3.20 thus becomes

$$y_{n\ell} = \varepsilon_{n\ell}^{(c)} (\varepsilon_{n\ell}^{(v)} \tilde{f}_{\ell}(\theta_n)) + \varepsilon_{n\ell}^{(a)} = \varepsilon_{n\ell}^{(m)} \tilde{f}_{\ell}(\theta_n) + \varepsilon_{n\ell}^{(a)}. \qquad (3.21)$$

**Observational geometry effects** – To account for beam dilution and for signal amplification or attenuation effects due to the angle with the surface of the cloud, we introduce a last multiplicative parameter $\kappa = (\kappa_n)_{n=1}^N$. Unlike the other uncertainty parameters, this parameter is to be inferred. It scales all lines identically so that Eq. 3.21 is rewritten

$$y_{n\ell} = \varepsilon_{n\ell}^{(m)} \kappa_n \tilde{f}_{\ell}(\theta_n) + \varepsilon_{n\ell}^{(a)}. \qquad (3.22)$$

Such a parameter was also introduced e.g., in Sheffer and Wolfire (2013), where $\kappa_n$ is defined as a product of 4 terms. Two of these terms are smaller than 1, one of which being the beam dilution factor. The two remaining are greater than 1. Similar multiplicative parameters were also considered in Joblin et al. (2018). We assume that for all $n$, $\log_{10} \kappa_n \in [-1, 1]$. To simplify notation, this parameter $\kappa_n$ is added to the parameter vector $\theta_n$ such that $\tilde{\mathbf{f}}(\theta_n)$. Adding this parameter does not affect the twice differentiable property of the forward model $\tilde{\mathbf{f}}$, as it is simple to compute $\partial \tilde{\mathbf{f}} / \partial \kappa_n(\theta_n)$ and higher order derivatives.
**Censorship** – For real data applications as in Chapter 6, we only have access to reduced multiline integrated intensity maps instead of the raw hyperspectral map. In these cases, a simpler observation model may be imposed. Upper values $\omega_{nl} \ge 0$ on observations are sometimes provided when $I_{nl} \ll \sigma_{a,nl}$. These upper values $\omega_{nl}$ on observations correspond to lower sensitivity bounds of the telescope. In these cases, we consider censorship as described in Chapter 3 (Eq. 3.10). Overall, the general observation model considered in this thesis is

$$y_{nl} = \max \left\{ \omega_{nl}, \varepsilon_{nl}^{(m)} \tilde{f}_{l}(\theta_{n}) + \varepsilon_{nl}^{(a)} \right\} = \begin{cases} \omega_{nl} & \text{if } \varepsilon_{nl}^{(m)} \tilde{f}_{l}(\theta_{n}) + \varepsilon_{nl}^{(a)} \le \omega_{nl} \\ \varepsilon_{nl}^{(m)} \tilde{f}_{l}(\theta_{n}) + \varepsilon_{nl}^{(a)} & \text{otherwise} \end{cases}. \quad (3.23)$$

The absence of censorship corresponds to the case where $\omega_{nl} = -\infty$.

## 3.5 Modeling, inference and model assessment choices

This section summarizes the choices we made in this thesis regarding statistical modeling, inference and model assessment. These choices will be detailed in Part II.

### Statistical modeling:

1. The forward model is set to the Meudon PDR code, as we assume that it simulates photodissociation regions accurately. Unlike lighter codes, this numerical simulator requires a few hours per evaluation. Evaluating it during the inference would result in an extremely slow inversion procedure. To accelerate inference and be able to compute gradients efficiently, we derive an accurate, fast, and light ANN emulator of the Meudon PDR code. This emulator is trained from a training dataset with a regular grid structure on the parameter space. The derivation of this emulator is the core of Chapter 4.

2. The uncertainty model from Eq. 3.23 combines two sources of noise: one additive and Gaussian, associated with thermal noise, and one multiplicative that includes both calibration and model misspecification errors. The additive Gaussian noise dominates in low SNR observations and becomes negligible in the high SNR regime. Conversely, the multiplicative lognormal noise dominates at high SNR observations and is negligible in the low SNR regime. Therefore, one cannot neglect one of the two noise sources. The resulting likelihood function has an expression challenging to address as is. In Chapter 5 (Section 5.1.1), we propose a closed-form and non-hierarchical approximation of the likelihood function with controlled error.

3. The proposed prior distribution combines a weakly-informative smooth uniform prior on a cube and an informative spatial regularization. The former encodes validity intervals for the components of $\Theta$. The latter exploits the map structure in $\Theta$. It is based on the $L_2$ norm of the maps Laplacian, as in Marchal et al. (2019). The resulting negative log-prior is also twice differentiable. The proposed prior is fully introduced in Chapter 5 (Section 5.1).

**Statistical inference** – Chapter 5 (Section 5.2) presents a new MCMC algorithm dedicated to map-structured data and high-dimensional multimodal distributions. The proposed sampler combines two kernels:

1. The first kernel addresses the regularity issue of the posterior. It performs efficient local exploration of the posterior by exploiting its local geometry, encoded by a RMSProp preconditioner.

2. The second kernel addresses the multimodality of the posterior. It permits escapes from the local minima by exploiting the map structure.

**Model assessment** – We use Bayesian hypothesis testing as Chevallard and Charlot (2016); Galliano et al. (2021); Lebouteiller and Ramambason (2022). In Chapter 5 (Section 5.3), we extend the test with three specificities:
1. The discrepancy measure $T$ is set to the negative log-likelihood for generality. The resulting $p$-value thus remains relevant for non-Gaussian noise models.

2. We compute one $p$-value per pixel to help identify potential regions where the Meudon PDR code does not simulate the physics accurately.

3. Uncertainties associated with the Monte Carlo estimators of the $p$-values are accounted for. The test is therefore made more robust to wrong decisions caused by the error inherent to MC estimation.

## 3.6 Conclusion

In this chapter, we reviewed statistical modeling and inference approaches adopted in the ISM community. We showed that the choices of model and inference method heavily depend on the community. Cosmology seems to be a pioneer in using state-of-the-art statistical approaches, e.g., with early uses of ANNs to emulate heavy numerical models and of advanced MCMC algorithms to perform inference. Dust studies in the ISM community are also quite advanced, with e.g., the use of sophisticated uncertainty models and the inference of prior hyperparameters with hierarchical models. In particular, works like Galliano (2018) or Galliano et al. (2021) infer physical parameters in very high dimensions – $O(10^4)$ – and with advanced noise models. To the best of our knowledge, these works are currently the only sampling-based applications with such dimensionality in the ISM community. The only other mentioned works relying on such dimensions are Paumard et al. (2014); Ciurlo et al. (2016); Paumard et al. (2022), which infer maps of parameters by exploiting spatial regularization within an optimization algorithm.

In comparison, inverse problems based on maps of integrated intensities of ionic, atomic and molecular lines have mainly been addressed for maps with a limited number of pixels. For instance, Sheffer et al. (2011); Sheffer and Wolfire (2013); Joblin et al. (2018); Ramambason et al. (2022) work on single-pixel observations. In Wu et al. (2018), the 176-pixel observation is handled with a pixel-by-pixel approach. In the coming years, with the upcoming JWST data and large observation surveys such as the IRAM-30m "ORION-B" Large Program, large observation maps with varying SNR are expected to be more common. The inverse problems we are interested in are therefore expected to increase in dimensionality.

In this thesis, we aim at inferring maps of physical parameter $\Theta$ with up to $O(10^4)$ pixels. We perform this inference using a state-of-the-art numerical model, the Meudon PDR code, and state-of-the-art statistical modeling and inference methods. The resulting likelihood function is highly non-linear, which causes the posterior distribution to be potentially multimodal. Usual sampling algorithms such as RWMH or the affine-invariant MCMC usually fail to explore such distributions and remain stuck in one local minimum. In addition, the code Meudon PDR is assumed to be twice differentiable, but not gradient Lipschitz continuous. This latter property makes the local exploration of the posterior complicated for classic gradient-based sampling algorithms such as MALA or HMC.

Chapter 4 derives a fast, light and accurate approximation of the Meudon PDR code. Chapter 5 (Section 5.1) details a likelihood function approximation, proposes a spatial regularization prior and the resulting posterior distribution. Chapter 5 (Section 5.2) presents the proposed MCMC algorithm used for inference. Chapter 5 (Section 5.4.2) presents a first very high dimensional synthetic case. Chapter 6 presents multiple applications of the proposed approach to real observations.
# Part II

## Solving inverse problems on ISM
multiline maps
<figure><img src="image_24.png" /></figure>
# Chapter 4

## Fast simulations of photodissociation region models

> "Learning with neural networks was proposed in the mid-20th century. It yields an effective learning paradigm and has recently been shown to achieve cutting-edge performance on several learning tasks."

Shalev-Shwartz and Ben-David (2014, chapter 20)

## Contents

<table><tr><td>4.1 Deriving emulators with interpolation or regression</td><td>86</td></tr><tr><td>4.1.1 Interpolation methods</td><td>87</td></tr><tr><td>4.1.2 Regression methods</td><td>87</td></tr><tr><td>4.2 Neural networks for regression</td><td>88</td></tr><tr><td>4.2.1 Generalities on neural networks</td><td>88</td></tr><tr><td>4.2.2 Fitting a neural network to a dataset</td><td>89</td></tr><tr><td>4.3 Dataset structure</td><td>90</td></tr><tr><td>4.4 Illustration: the log Rosenbrock function</td><td>92</td></tr><tr><td>4.5 Towards an emulator of the Meudon PDR code</td><td>96</td></tr><tr><td>4.5.1 Datasets</td><td>96</td></tr><tr><td>4.5.2 Comparison metrics</td><td>98</td></tr><tr><td>4.6 Designing and training adapted ANNs</td><td>100</td></tr><tr><td>4.6.1 Removing outliers from the training set</td><td>100</td></tr><tr><td>4.6.2 Exploiting correlations between line intensities</td><td>102</td></tr><tr><td>4.6.3 A polynomial transform to learn nonlinearities</td><td>104</td></tr><tr><td>4.6.4 Dense networks to reuse intermediate computations</td><td>106</td></tr><tr><td>4.7 Experiments: application to the Meudon PDR code</td><td>107</td></tr><tr><td>4.7.1 Performance analysis</td><td>107</td></tr><tr><td>4.7.2 Removing outliers is crucial</td><td>109</td></tr><tr><td>4.7.3 The importance of the polynomial feature augmentation</td><td>109</td></tr><tr><td>4.8 Conclusion</td><td>110</td></tr><tr><td>Appendix 4.A Automatic outlier detection procedure</td><td>111</td></tr><tr><td>Appendix 4.B Content of clusters of lines</td><td>112</td></tr></table>
In Chapter 3 (Section 3.1.1), we reviewed how numerical forward models are handled in studies of the interstellar medium (ISM). In case of a simple and fast astrophysical model, the numerical simulation is sometimes directly used within the inversion, e.g., in Holdship et al. (2018). The inverse problem studied in this thesis – to be fully introduced in Chapter 5 – involves large observation maps. Each pixel of these maps is modeled using the Meudon PDR code. As the Meudon PDR code requires a few hours for each evaluation, evaluating the associated likelihood during the inversion procedure would be prohibitively slow. We farther described some simulation-based inference approaches, including approximate Bayes computing methods and likelihood emulation.

This chapter is focused on emulating the Meudon PDR code from a set of evaluations instead of the full likelihood function. This approach avoids re-training an emulator whenever the nature or level of the noise changes. The resulting surrogate model will be included in the observation model detailed in Chapter 5. A similar strategy was adopted e.g., in Wu et al. (2018) and Ramambason et al. (2022). Using an emulator instead of the original numerical code induces approximation errors that should be minimized.

In this chapter, we propose a neural network-based emulator, and compare it with standard interpolation methods in terms of speed, memory requirements and accuracy. We eventually obtain a fast, light and accurate artificial neural network (ANN) emulator. This emulator is used as a forward model in Chapter 5.

This chapter is based on the journal article Palud et al. (2023c) and the associated Gretsi conference article Einig et al. (n.d.). This article is the product of an equal collaboration with Lucas Einig, another PhD candidate of the ORION-B consortium.

The code associated with this chapter and its experiments was published in a `PYTHON` package called `NNBMA` (for Neural Network-Based Model Approximation). The code can be found at https://github.com/einigl/ism-model-nn-approximation. It can be installed via `PIP` – see https://pypi.org/project/nnbma/. The associated documentation can be found at https://ism-model-nn-approximation.readthedocs.io/en/latest/? badge=latest.

Section 4.1 provides a description of the considered interpolation and regression methods. Section 4.2, introduces ANNs and how to fit them to a grid of precomputed models. Section 4.3 presents some methods to generate datasets of precomputed models with a good coverage of the parameter space. The illustrative Section 4.4 compares all considered methods and dataset structures on a simple case, the log Rosenbrock function. Section 4.5 describes the dataset of precomputed models and introduces the metrics used to compare surrogate models. In Section 4.6, we design ANNs that address the specificities of ISM numerical codes. Section 4.7 compares these ANNs with classic interpolation methods in terms of speed, memory requirements and accuracy.

## 4.1 Deriving emulators with interpolation or regression

Deriving an emulator consists in estimating the function $\tilde{f}: \mathbb{R}^D \to \mathbb{R}^L$ that maps input vectors $\theta$ to output vectors $y = f(\theta)$ as closely as possible. This function $\tilde{f}$ is *learned*, or equivalently, *fitted*, from a dataset of precomputed models $\mathcal{D} = \{(\theta_n, y_n) \in \mathbb{R}^D \times \mathbb{R}^L, n = 1, \dots, N\}$. We remind that in this thesis, the input vector $\theta$ corresponds to a vector of physical parameters, e.g., temperature, thermal pressure, volume density, and the output vector $y$ to observables computed by a numerical code, e.g., integrated intensities of specific emission lines.

In this Section, we first summarize four interpolation methods we will use in our comparison. Then, we present ANNs and the associated regression methods. For a more detailed introduction to ANNs, see Shalev-Shwartz and Ben-David (2014, chapter 20).
4.1.1 Interpolation methods

Interpolation methods have become a common approach to build surrogates of comprehensive ISM models over the last years thanks to their conceptual and implementation simplicity. Four families of interpolation methods are usually considered:

- Nearest-neighbor interpolation assigns to a new point the value of the closest point in the dataset. It is fast but generally performs poorly in terms of accuracy. It was used e.g., in Ramambason et al. (2022).

- Piecewise linear interpolation generally performs better, while remaining quite fast. It first triangulates the dataset, so that a new point is associated with a cell of the triangulation. Then, it returns a weighted average of the cell points values. It was used e.g., in Thomas et al. (2018).

- Spline interpolation (Bojanov et al., 1993) methods are based on piecewise polynomials, yielding an even more accurate and still fast surrogate model. It was used e.g., in Blanc et al. (2015).

- Radial basis function (RBF) interpolation (Fasshauer, 2007, chapter 6) uses the full dataset for each evaluation. For a new point, it returns a weighted sum of the values of all the dataset points, where the weights depend on the distance to this new point. Surrogate models defined with RBF interpolation are generally very accurate but slow. It was used e.g., in Wu et al. (2018).

Interpolation methods suffer from some drawbacks. By definition, a surrogate model defined with an interpolation method passes exactly through the points of the dataset. This constraint does not guarantee good accuracy on points not used during the fit. Besides, evaluating a surrogate model defined with an interpolation method requires loading the whole dataset, which can induce a very heavy memory cost when it contains many precomputed models or many quantities associated with each model. Finally, although they are generally faster than the original numerical codes, interpolation methods handle outputs (i.e., observables) independently. They are thus quite slow when the number of outputs is large.

Accuracy can be improved by relaxing the constraint of passing through the points of the dataset. The memory and speed drawbacks can be addressed using an approach that allows predicting all outputs at once.

4.1.2 Regression methods

Relaxing the constraint of passing through the points of the dataset turns the interpolation problem into a regression problem. The first step towards emulating a numerical model $\mathbf{f}$ is to restrict the search of $\tilde{\mathbf{f}}$ to a class of functions. Most classes of functions are parametrized with vectors $\psi$. In the following, potential surrogate functions are sometimes denoted $\tilde{\mathbf{f}}_{\psi}$ to emphasize this point. For instance, in linear regression, an affine function $\tilde{\mathbf{f}}_{\psi} : \theta \mapsto \mathbf{W}\theta + \mathbf{b}$ is uniquely described by $\psi = (\mathbf{W}, \mathbf{b})$. Given the complexity of ISM numerical models, the class of affine functions is too restrictive to produce accurate surrogate models, and richer classes are required.

Multiple classes of functions and the associated regression algorithms enable the emulation of complex non-linear functions from data of precomputed models, such as polynomial functions, k-nearest-neighbor regression (Shalev-Shwartz and Ben-David, 2014, chapter 19) – used e.g., in Smirnov-Pinchukov et al. (2022) –, Gaussian process regression (Rasmussen and Williams, 2006), decision trees and the associated ensemble methods such as random forest (RF) (Shalev-Shwartz and Ben-David, 2014, chapter 18) – used e.g., in Bron et al. (2021) – or XGBoost (Chen and Guestrin, 2016), ANNs (Shalev-Shwartz and Ben-David, 2014, chapter 20). All methods based on decision trees or nearest-neighbors yield piecewise functions, which prevents from enforcing desirable regularity property in the surrogate model such as continuity or differentiability. Besides, all the listed algorithms, except ANNs and nearest-neighbor interpolation, handle multiple outputs independently, which slows prediction when the number of outputs is high. Conversely,
ANNs predict all outputs at once from a common sequence of intermediate computations, which is considerably faster. In addition, ANNs are known to yield very accurate surrogate models both in theory and in practice (Shalev-Shwartz and Ben-David, 2014, chapter 20). Finally, an ANN comes with the ability to automatically and efficiently compute the first or second order derivative of its outputs with respect to its inputs $\theta$, using automatic differentiation (Paszke et al., 2017). In other words, for any output $\ell$ of an ANN, one can efficiently access the gradient $\nabla \tilde{f}_{\ell}$ or the Hessian matrix $\nabla^2 \tilde{f}_{\ell}$. This property is desirable to efficiently explore high dimensional parameter space in inverse problems, as detailed in Chapter 3 (Section 3.1.1). We consequently adapt the rich and versatile class of ANNs to address the complexity of ISM numerical models, to exploit prior knowledge on the regularity of the function to approximate, and to efficiently predict all outputs at once.

## 4.2 Neural networks for regression

The class of ANNs is first introduced. We then explain how to fit ANNs to a dataset.

### 4.2.1 Generalities on neural networks

Artificial neural networks (ANNs) form a class of mathematical models inspired from biological neural systems. The first ANN was proposed in McCulloch and Pitts (1943) to perform logical operations. Since then, multiple hardware and algorithmic developments such as GPU computing and back-propagation (Rumelhart et al., 1986) made them capable of learning more complex patterns and relationships in data. They enjoy strong theoretical results. For different sets of assumptions on the architecture, universal approximation theorems establish that ANNs can approximate almost any continuous function (Hornik et al., 1989; Leshno et al., 1993). They gained a widespread popularity after the 2012 ImageNet Challenge, an image classification competition in which an ANN significantly outperformed competing methods. They are nowadays state-of-the-art for a variety of tasks in vector, image, sound or text processing in multiple scientific and industrial fields. As we showed in Chapter 3 (Section 3.1.1), ANNs are already used in astrophysics, including ISM studies, to emulate numerical models (Grassi et al., 2011; de Mijolla et al., 2019; Holdship et al., 2021; Grassi et al., 2022).

**Formal description of an ANN** – Throughout this thesis, an ANN is considered as a function $\tilde{f}: \theta \in \mathbb{R}^D \to \tilde{f}(\theta) \in \mathbb{R}^L$, where $D$ and $L$ are input and output dimensions, respectively. For a numerical model, $D$ is the number of considered physical parameters, e.g., thermal pressure or visual extinction, and $L$ is the number of predicted observables, e.g., line intensities. An ANN is made of $h+1$ intermediate functions $\tilde{f}^{(j)}$, called *layers*. Intermediate layers $1 \le j \le h$ are called the *hidden layers* and the final layer is the *output layer*. The $j^{th}$ layer takes an intermediate vector $\theta^{(j)} \in \mathbb{R}^{i_j}$ as input and computes an intermediate output $y^{(j)} \in \mathbb{R}^{o_j}$. Intermediate dimensions $i_j$ and $o_j$ can be chosen arbitrarily, except for $i_1 = D$ and $o_{h+1} = L$. In a feedforward neural network, layers interact following an acyclic graph. The output of a layer $j$ feeds one or more of the next layers $j' > j$, hence the notion of direction in a feedforward neural network.

Figure 4.1 shows the structure of a simple ANN that contains $h = 2$ hidden layers and one output layer. This ANN takes in input $D = 2$ physical parameters and predicts $L = 10$ observables. It is an instance of a feedforward neural network as each layer only feeds layers closer to the output layer, as shown on the left-hand side of the figure. More precisely, its layer graph is linear: the output of one of its layers $j$ is the input of the next layer $j + 1$, as the $\theta^{(j+1)} = y^{(j)}$ and $i_{j+1} = o_j$ equalities show. Alternative feedforward architectures with non-linear layer graph exist, such as residual networks (He et al., 2016) and dense networks (Huang et al., 2017). These architectures include skip connections between layers that bypass the activation function to preserve original input information and intermediate computations. However, linear layer graphs remain the simplest and most widespread multi-layer architectures for vector classification or regression
tasks. In the remainder of this chapter, ANNs are considered with such an architecture unless mentioned otherwise.

<figure><img src="image_25.png" /><figcaption>Figure 4.1 – Structure of a simple feedforward neural network with h = 2 hidden layers and a linear layer graph, shown on the left.</figcaption></figure>

A hidden layer combines an affine transformation and a non-linear scalar function g(j) applied element-wise:

$$
\tilde{\mathbf{f}}^{(j)} : \theta^{(j)} \mapsto \mathbf{y}^{(j)} = g^{(j)}(\mathbf{W}^{(j)}\theta^{(j)} + \mathbf{b}^{(j)}), \quad (4.1)
$$

with $\mathbf{W}^{(j)} \in \mathbb{R}^{o_j \times i_j}$ and $\mathbf{b}^{(j)} \in \mathbb{R}^{o_j}$ the weight matrix and bias vector of the affine transformation, respectively. The non-linear scalar function $g^{(j)}$ is called an activation function. Common activation functions include the sigmoid, hyperbolic tangent, rectified error linear unit (ReLU), and multiple variants – see Nwankpa et al. (2021) for a review. Choosing a specific activation function $g^{(j)}$ for each of the h hidden layers might lead to better performance but would require training many ANNs. A unique g is therefore generally set for all hidden layers. As an ANN is a composition of affine transformations and activation functions, the regularity properties of the activation function g apply to the full ANN. For instance, if g is differentiable, so is the full ANN $\tilde{f}$. The same implication holds for infinitely continuously differentiable functions.

The output layer transforms the outputs of one or more hidden layers into the desired pre-
diction with an affine transformation and a second activation function. This second activation
function depends on the considered problem. The sigmoid and softmax functions are usually
employed to return probabilities in binary and multi-class classification, respectively. In regression
tasks, the identity function is generally used.

**Architecture and learnable parameters** – Overall, in a regression context, an ANN is uniquely defined by its layer graph, an activation function g, a number of hidden layers h ≥ 0, a sequence of sizes of its layers (i_j, o_j)_{j=1}^{h+1}, and the sequence of weight matrices and bias vectors ψ = (W^{(j)}, b^{(j)})_{j=1}^{h+1}. All but ψ are manually set by the user prior to the fitting procedure. They define the so-called architecture of the ANN. The only learnable parameters are the weight matrices and bias vectors ψ = (W^{(j)}, b^{(j)})_{j=1}^{h+1}. The fitting procedure adjusts them so that $\tilde{f}_{\psi} \simeq f$.

4.2.2 Fitting a neural network to a dataset

In regression, once the class of function is set, the associated parameter $\psi$ is adjusted so that $\tilde{f}_{\psi}$ fits the dataset $\mathcal{D}$ of precomputed models. In our case, the class of functions is set by selecting an ANN architecture, and the associated parameter $\psi$ contains the weights and biases of the network's layers. A loss function $\mathcal{L}(\tilde{f}; \mathcal{D})$ quantifies the distance between predictions $\tilde{f}_{\ell}(\theta)$
and the corresponding true values $y_\ell$. It is based on an error function, e.g., the absolute error (AE)

$$AE(\tilde{\mathbf{f}}; (\boldsymbol{\theta}, y_\ell)) = |\tilde{f}_\ell(\boldsymbol{\theta}) - y_\ell|, \qquad (4.2)$$

or the squared error (SE)

$$SE(\tilde{\mathbf{f}}; (\boldsymbol{\theta}, y_\ell)) = (\tilde{f}_\ell(\boldsymbol{\theta}) - y_\ell)^2. \qquad (4.3)$$

The loss function summarizes the set of $N \times L$ errors obtained on the dataset $\mathcal{D}$. The mean is often used for computational efficiency of evaluation and differentiation, yielding e.g., the mean squared error (MSE) or the mean absolute error (MAE). Obtaining the best function $\tilde{\mathbf{f}}$ boils down to minimizing the loss function with respect to the parameter $\boldsymbol{\psi}$

$$\tilde{\mathbf{f}} \in \arg \min_{\boldsymbol{\psi}} \mathcal{L}(\tilde{\mathbf{f}}_{\boldsymbol{\psi}}; \mathcal{D}). \qquad (4.4)$$

In general, problems of the form of Eq. 4.4 do not admit a closed-form solution. Furthermore, with ANNs, the loss function $\mathcal{L}(\tilde{\mathbf{f}}_{\boldsymbol{\psi}}; \mathcal{D})$ is generally not convex, with multiple saddle points and local minima (Shalev-Shwartz and Ben-David, 2014, chapter 20). When $\boldsymbol{\psi}$ is low dimensional, such problems can be solved approximately using a meta-heuristic such as a genetic algorithm or simulated annealing – see Section 2.A.1. As ANNs typically contain at least hundreds of parameters to tune, meta-heuristic methods are prohibitively slow.

In contrast, gradient descent (GD) methods are computationally very efficient. They rely on auto-differentiation to efficiently evaluate the gradient of the loss function $\nabla_{\boldsymbol{\psi}}\mathcal{L}$ and on back-propagation (Rumelhart et al., 1986) to efficiently update $\boldsymbol{\psi}$. Stochastic gradient descent (SGD) Shalev-Shwartz and Ben-David, 2014, chapter 14 accelerates the search by using batches instead of the full dataset in gradient evaluations. Preconditioned variants such as RM-SProp (Tieleman and Hinton, 2012) or Adam (Kingma and Ba, 2017) exploit the local geometry of the loss function to escape from saddle points and farther accelerate convergence to a good local minimum. This optimization procedure is often called training or learning phase with ANNs, because the network progressively learns from the data as the loss function decreases.

## 4.3 Dataset structure

Like the choice of the class of function, the structure of the training dataset $\mathcal{D}$ in the $\boldsymbol{\theta}$-space can have a dramatic impact on the accuracy of the obtained emulator. To approximate the Meudon PDR code, we propose to resort to a lattice structure on a cube of dimension $D$ generated prior to any fit. This approach is very common in ISM studies (Joblin et al., 2018; Wu et al., 2018; Ramambason et al., 2022) – see Chapter 3.

**Lattice structured datasets** – A regular grid structure has many advantages. First, if the parameter space is a $D$-dimensional cube and the number of points $N$ is such that $N = K^D$ for some $K \ge 1$, the lattice structure maximizes the minimum distance between two points of the dataset. Therefore, this structure enforces good coverage of the parameter space. Second, it is often more convenient to manually inspect a dataset with such a structure, especially for $D \ge 3$, when the full dataset cannot be visualized at once. Besides, it allows the use of efficient interpolation methods such as splines, for which regular grid structures are mandatory. Also, the regularity of the grid can be exploited to accelerate nearest-neighbor and linear interpolations. Indeed, it facilitates the localization of the new point within the dataset and thus the identification of the points to use for prediction. However, this structure is not necessary for RBF interpolation methods or in regression approaches.

Despite these advantages, the lattice structure may not yield optimal emulators, and other structures might yield better accuracy. For instance, for coarse grids, a lattice structure creates
wide regions in the parameter space without any point in the dataset $D$. An emulator fitted in this context could struggle to reproduce accurately brutal changes occurring in these regions. Besides, it creates a strong anisotropy in the cube, which can reduce the accuracy of an emulator. For instance, large gradient variations occurring between two slices of the lattice would be difficult to emulate accurately.

**Alternative methods to generate datasets with good coverage** – Many alternative methods enable sampling a $D$-dimensional cube. The simplest would be to sample independent and identically distributed (i.i.d.) points from a uniform distribution on the cube. However, this method would leave some regions uncovered, and some pairs of points could be arbitrarily close, which limits the information brought by each. Better methods enforce repulsion between points to ensure a good coverage of the cube. Some of these methods are deterministic, with low discrepancy sequences such as the Halton sequences used in quasi-Monte Carlo (QMC) methods (Asmussen and Glynn, 2007, chapter 9). Some others are partially random, such as latin hypercube sampling (LHS) (McKay et al., 1979), its variants orthogonal LHS (OLHS) (Tang, 1993) and symmetric LHS (SLHS) (Ye et al., 2000), or stratified Monte Carlo (Haber, 1966). These methods are simple to implement. To generate $N$ points in the unit cube $[0, 1]^D$, LHS handles each dimension independently. For each dimension, it divides the $[0, 1]$ segment in $N$ smaller segments of equal size, draws one sample per segment with a uniform distribution, and randomly permutes the draws. Stratified Monte Carlo divides the cube in $N$ smaller regions of equal volume and draws one sample from a uniform distribution on each smaller region. Orthogonal LHS combines both LHS and stratified MC to ensure better coverage. Symmetric LHS is a generalization of OLHS that enforces symmetry properties to farther improve the coverage of the cube. Other LHS algorithms that optimize criteria such as the entropy or the minimum distance between two points also exist - see Joseph (2016) for a review.

Figure 4.2 illustrates how each method would cover the unit square, i.e., for $D = 2$. The i.i.d. draws with uniform distribution on square yield the visually worst datasets. The lattice structure best covers corners and edges, but leaves large gaps between points. LHS and stratified MC lead to datasets with better coverage of the square. By combining them, orthogonal LHS yields even better datasets.

As illustrated in the next section, these methods may yield better accuracy than the lattice structure. However, the associated accuracy has a high variance, as it depends on the presence or absence of points in regions with high gradients. Using the lattice structure is therefore a good heuristic.

**Iterative methods** resort e.g., to Gaussian process regression – also known as Kriging – to progressively recommend new points to evaluate. Such methods could lead to even more cleverly structured datasets, and yield better accuracy surrogate models with fewer data points. However, they require re-fitting the surrogate model for each new point or batch of points. It can be easily implemented for interpolation methods. It is computationally more expensive for regression approaches, where the surrogate model is defined as solution of an optimization problem. In addition, it is not straightforward how to use such an approach in the case of the Meudon PDR code, considered in this chapter. The code predicts numerous outputs, which adds to the computational cost for Gaussian processes. Some outputs of the code cover more decades than others: a naive iterative approach would favor those outputs, which are not necessarily the most physically interesting. Finally, the Meudon PDR code may produce un-physical outliers – see Section 4.5.1. An iterative method favoring sampling large error regions could result in an over-representation of input combinations leading to outliers, which are to be eliminated in Section 4.6.1. For these reasons, iterative approaches are not considered in this chapter.
<figure><img src="image_26.png" /><figcaption>Figure 4.2 – Illustration of methods to generate non-lattice datasets of $N = 9$ points. The sampling principle of LHS, stratified MC and orthogonal LHS are indicated with the grids.</figcaption></figure>

## 4.4 Illustration: the log Rosenbrock function

In this Section, we illustrate the four interpolation methods mentioned in Section 4.1.1 and two regression methods, namely random forest (RF) and ANNs, on the log Rosenbrock function (Rosenbrock, 1960):

$$ \ln R : \theta \in \mathbb{R}^2 \mapsto \ln \left[ 1 + (1 - \theta_1)^2 + 100 \left( \theta_2 - \theta_1^2 \right)^2 \right], \quad (4.5) $$

which is positive and admits a minimum at $(1, 1)$ such that $\ln R(1, 1) = 0$. We restrict ourselves to the $[-1, 1] \times [-0.5, 1.5]$ square.

The SCIPY<fnref n="1" /> (Virtanen et al., 2020) implementation is used for the interpolation methods. For RBF interpolation, the cubic kernel is used as it yields the best results. For random forest and ANNs, the SCIKIT-LEARN<fnref n="2" /> (Pedregosa et al., 2011) implementations are used. Note that the PYTORCH implementation is used in the remainder of this thesis. The SCIKIT-LEARN implementation is used in this example for simplicity. The ANN architecture is set to a perceptron with $H = 2$ hidden layers of 32 neurons each. The Rosenbrock function is smooth, i.e., infinitely continuously differentiable. The tanh activation function is used to enforce this smoothness property in the emulator.

Figure 4.3 shows the values taken by this function on this subspace. Its banana shape being relatively thin, it is challenging to reproduce accurately with a limited amount of points. The results are compared over six dataset structures, namely lattice, Halton sequence, i.i.d. uniform, LHS, stratified MC and orthogonal LHS. Each training dataset contains $N = 49$ points. The accuracies of emulators are evaluated on a $151 \times 151$ grid on the $[-1, 1] \times [-0.5, 1.5]$ square, using the maximum absolute error, the MAE and the MSE.

Figure 4.4 shows the emulator obtained for each method and the associated errors. All methods struggle to reproduce accurately the banana shape of the Rosenbrock function, as there are few points in the training set. In ISM numerical models, such strong and fast variations can

<footnote n="1">1. https://scipy.org/</footnote>

<footnote n="2">2. https://scikit-learn.org/stable/index.html</footnote>
<figure><img src="image_27.png" /><figcaption>Figure 4.3 – log Rosenbrock function $\ln R$ (Eq. 4.5).</figcaption></figure>

correspond to a change of regime and are thus of crucial importance. The overall performance of an emulator highly relies on the presence of a few points in such regions. The presence or absence of these key points in datasets generated by random methods induces a high variance in the accuracy of emulators built from these datasets.

Table 4.1 quantifies the results of this comparison with three metrics: the maximum absolute error, the MAE and the MSE. In all scenarios, ANNs and RBF interpolation yield the most accurate emulators. Besides, ANNs provide the best emulator with respect to each of the three metrics. However, the best results are not always obtained with a lattice structure.

Overall, this illustrative example shows three important aspects of the emulator derivation of the Meudon PDR code:

1. ANNs often yield much more accurate emulators than interpolation methods.

2. The lattice structure might not always lead to optimal emulators. However, as it covers the square very well, it still leads to good emulators. Besides, as the other considered methods are random by nature and as only a few training points have a great impact on the emulator accuracy, two datasets generated with the same method can lead to emulators with very different accuracies. To limit uncertainties on causes of the quality of the results, we prefer to exploit a lattice dataset in our emulation of the Meudon PDR code.

3. The result of this comparison thus depends on the error function. Choosing a relevant accuracy metric is therefore crucial to derive a good emulator. In Section 4.5.2, we introduce a new accuracy metric, the Error Factor, tailored to the Meudon PDR code.
<figure><img src="image_28.png" /><figcaption>Figure 4.4 – Comparison of four popular interpolation methods and two regression methods, namely random forest (RF) and ANN, on the log Rosenbrock function ln R (Eq. 4.5). For each dataset structure, the first row shows the emulator obtained with the considered algorithm, and the second row shows the absolute error. Cubic splines are only defined on a lattice. For linear interpolation, white values indicate that the emulator is not defined for the corresponding value of θ.</figcaption></figure>
<table><caption>Table 4.1 – Comparison of popular interpolation and regression methods on the log Rosenbrock function with datasets of N = 49 points. The values correspond to the datasets and emulators displayed in Figure 4.4. Bold values indicate the best emulator class for each training dataset structure, and blue values indicate the best couple (structure, emulator class) for each accuracy metric. As cubic splines can only be defined on a lattice, there is no value for the other dataset structures. As linear interpolation is only defined inside the convex envelope of the training dataset, its error can be computed on a fraction of the test set for all structures except the lattice. For this reason, its values are indicated in parentheses.</caption><thead><tr><th rowspan="2">Structure</th><th colspan="4">Interpolation</th><th colspan="2">Regression</th></tr><tr><th>nearest neighbor</th><th>linear</th><th>cubic spline</th><th>RBF</th><th>random forest</th><th>neural network</th></tr></thead><tbody><tr><td colspan="7" style="text-align:center; font-weight:bold;">Max absolute error</td></tr><tr><td>Lattice</td><td>2.974</td><td>1.880</td><td>1.611</td><td>1.568</td><td>2.754</td><td style="color:blue;">0.787</td></tr><tr><td>i.i.d. uniform</td><td>4.404</td><td>(3.120)</td><td>-</td><td>2.742</td><td>4.001</td><td style="color:blue;">2.104</td></tr><tr><td>QMC</td><td>2.585</td><td>(1.909)</td><td>-</td><td style="color:blue;">1.766</td><td>2.850</td><td>2.742</td></tr><tr><td>LHS</td><td>3.740</td><td>(3.259)</td><td>-</td><td>1.860</td><td>3.371</td><td style="color:blue;">1.098</td></tr><tr><td>Stratified MC</td><td>3.148</td><td>(3.555)</td><td>-</td><td style="color:blue;">1.853</td><td>3.613</td><td>3.106</td></tr><tr><td>Orthogonal LHS</td><td>3.173</td><td>(3.378)</td><td>-</td><td>2.127</td><td>3.225</td><td style="color:blue;">1.736</td></tr><tr><td colspan="7" style="text-align:center; font-weight:bold;">mean absolute error (MAE)</td></tr><tr><td>Lattice</td><td>0.419</td><td>0.185</td><td>0.151</td><td>0.131</td><td>0.399</td><td style="color:blue;">0.082</td></tr><tr><td>i.i.d. uniform</td><td>0.516</td><td>(0.332)</td><td>-</td><td>0.225</td><td>0.653</td><td style="color:blue;">0.163</td></tr><tr><td>QMC</td><td>0.391</td><td>(0.167)</td><td>-</td><td>0.138</td><td>0.528</td><td style="color:blue;">0.133</td></tr><tr><td>LHS</td><td>0.431</td><td>(0.302)</td><td>-</td><td>0.139</td><td>0.544</td><td style="color:blue;">0.072</td></tr><tr><td>Stratified MC</td><td>0.406</td><td>(0.162)</td><td>-</td><td style="color:blue;">0.138</td><td>0.566</td><td>0.177</td></tr><tr><td>Orthogonal LHS</td><td>0.407</td><td>(0.190)</td><td>-</td><td>0.136</td><td>0.536</td><td style="color:blue;">0.083</td></tr><tr><td colspan="7" style="text-align:center; font-weight:bold;">mean squared error (MSE)</td></tr><tr><td>Lattice</td><td>0.324</td><td>0.094</td><td>0.060</td><td>0.052</td><td>0.280</td><td style="color:blue;">0.022</td></tr><tr><td>i.i.d. uniform</td><td>0.547</td><td>(0.343)</td><td>-</td><td>0.177</td><td>0.788</td><td style="color:blue;">0.092</td></tr><tr><td>QMC</td><td>0.276</td><td>(0.079)</td><td>-</td><td style="color:blue;">0.062</td><td>0.480</td><td>0.098</td></tr><tr><td>LHS</td><td>0.351</td><td>(0.276)</td><td>-</td><td>0.080</td><td>0.549</td><td style="color:blue;">0.025</td></tr><tr><td>Stratified MC</td><td>0.313</td><td>(0.074)</td><td>-</td><td style="color:blue;">0.065</td><td>0.529</td><td>0.152</td></tr><tr><td>Orthogonal LHS</td><td>0.307</td><td>(0.122)</td><td>-</td><td>0.069</td><td>0.526</td><td style="color:blue;">0.034</td></tr></tbody></table>
## 4.5 Towards an emulator of the Meudon PDR code

In this Section, we describe the setup used to define and fit emulators of the Meudon PDR code (Le Petit et al., 2006), introduced in Chapter 1. We first present the train and test sets, and then introduce the metrics used to compare emulators.

### 4.5.1 Datasets

**Physical parameters** – We approximate the Meudon PDR code with respect to the *D* = 4 input parameters that are most relevant for inference (Wu et al., 2018). The three main ones are the thermal pressure *P*<sub>th</sub>, the scaling factor *G*<sub>0</sub> of the interstellar standard radiation field and the size of the slab of gas measured in total visual extinction *A*<sub>V</sub><sup>tot</sup>. As in Wu et al. (2018), we consider a wide variety of environments with *P*<sub>th</sub> ∈ [10<sup>5</sup>, 10<sup>9</sup>] K cm<sup>-3</sup>, *G*<sub>0</sub> ∈ [1, 10<sup>5</sup>] and *A*<sub>V</sub><sup>tot</sup> ∈ [1, 40] mag. The Meudon PDR code computes line intensities for multiple angles φ between the cloud surface and the line of sight. In the Meudon PDR code, this angle φ can cover a [0, 60] deg interval. A face-on geometry corresponds to φ = 0 deg, and φ = 60 deg is the closest to an edge-on geometry. To enable analyses of PDRs with known edge-on geometry such as the Orion Bar (Joblin et al., 2018), this angle is added to the considered physical parameters. Table 4.2 details the considered ranges of the main input parameters. The secondary parameters of the code, listed in Section 1.3, are set to their default values – see Table 1.3.

<table><caption>Table 4.2 – Input parameters in the Meudon PDR code and structure of training dataset.</caption><thead><tr><th rowspan="2">Parameter</th><th rowspan="2">Range of values</th><th rowspan="2">Unit</th><th colspan="2">Training set lattice structure</th></tr><tr><th>Number of points</th><th>Grid</th></tr></thead><tbody><tr><td>Thermal pressure P<sub>th</sub></td><td>[10<sup>5</sup>, 10<sup>9</sup>]</td><td>K cm<sup>-3</sup></td><td>14</td><td>log spacing</td></tr><tr><td>Radiative intensity G<sub>0</sub></td><td>[1, 10<sup>5</sup>]</td><td>Mathis</td><td>14</td><td>log spacing</td></tr><tr><td>Total visual extinction A<sub>V</sub><sup>tot</sup></td><td>[1, 40]</td><td>mag</td><td>14</td><td>log spacing</td></tr><tr><td>Line-of-sight angle φ</td><td>[0, 60]</td><td>deg</td><td>7</td><td>linear spacing</td></tr></tbody></table>

**Datasets and first preprocessing** – Two datasets of Meudon PDR code evaluations are defined: a training set and a test set<sup>3</sup>.

The **training set** is used to fit all surrogate models. It contains $N_{train} = 19208$ points, structured as a 14 × 14 × 14 × 7 uniform regular grid on ($\log_{10} P_{th}$, $\log_{10} G_0$, $\log_{10} A_V^{tot}$, $\phi$). Table 4.2 shows its lattice structure. This uniform grid structure is chosen to include spline interpolation in the comparison and to simplify the outlier identification procedure presented in Section 4.6.1. The Meudon PDR code predicts line intensities that are strictly positive and span multiple decades, from $10^{-50}$ to $10^{-2}$ erg cm<sup>-2</sup> sr<sup>-1</sup> s<sup>-1</sup>. To avoid giving more weight to high line intensities in the regression and disregarding the lowest ones, we consider the intensities in log scale. In other words, we derive emulators to reproduce the log-intensities $\ln y \in \mathbb{R}^L$, i.e., we derive emulators of $\ln f$.

Similarly, $P_{th}$, $G_0$, and $A_V^{tot}$ are considered in log scale. Even in log scale, the parameters of interest cover intervals with quite different sizes. For instance, $\log_{10} G_0 \in [0, 5]$ while $\log_{10} A_V^{tot} \in [0, 1.602]$, i.e., $A_V^{tot}$ covers an interval more than 3 times smaller than $G_0$. Both interpolation methods and neural network-based regression typically suffer from this difference. All *D* parameters are thus standardized to have a zero mean and a unit standard deviation. This simple transformation generally improves accuracy for both interpolation methods and ANNs (Shalev-Shwartz and Ben-David, 2014, chapter 25).

The **test set** is used to assess the accuracy of surrogate models on data not used in the training step. It contains $N_{test} = 3192$ points. These points are generated with 456 independent

<footnote n="3">3. Both datasets can be found in https://ism.obspm.fr/files/ArticleData/2023_Palud_Einig/2023_Palud_Einig_data.zip.</footnote>
random draws from a uniform distribution on the (log₁₀ *P*<sub>th</sub>, log₁₀ *G*₀, log₁₀ *A*<sub>V</sub><sup>tot</sup>) cube and with a regular grid of 7 values on φ. To ensure consistent preprocessing between the two sets, both the input values θ and output values y of the test set undergo the same transformations as the training set. In particular, the standardization applied to its input values θ relies on the means and standard deviations obtained on the training set, and its output values y are considered in log₁₀-scale.

**Outliers** – Numerical codes may yield numerical instabilities. Such behavior was reported, e.g., in Auld et al. (2007) about the WMAP 3-yr likelihood code. In its domain of validity, the Meudon PDR code produces few of them. However, the considered complex non-linear physics can also lead to physical bistabilities or multistabilities. For example, the H₂ heating process can produce bistable solutions (Burton et al., 1990; Röllig and Ossenkopf-Okada, 2022). In such a case, spatial profiles computed by the code, e.g., of a species density or of the gas temperature, can oscillate between the possible solutions at each position in the modeled cloud. The line integrated intensities computed from these profiles can contain errors of up to a factor of 100 and thus are highly unreliable. Figure 4.5 shows one example of bistability and of its impact on one of the integrated intensities. This case can be easily detected as the bistability affects the temperature profile of the cloud, which is used to compute the integrated intensity of many emission lines. More often, only a handful out of the *L* lines show such pathological behavior.

<figure><img src="image_29.png" /><figcaption>(a) Temperature profiles in five simulated clouds, including the one with clear bistability.</figcaption></figure>

<caption>(b) Two-dimensional slice of the training set for one line *ℓ*. The point indicated by a black square corresponds to the bistable profile on the left.</caption>

Figure 4.5 – Illustration of an invalid integrated intensity due to a bistability in temperature profile. The corresponding point should be withdrawn from the training set.

The code being deterministic, an input vector θ consistently leads to a unique output vector y. However, in the regions of the parameter space with such multistabilities, variations of intensities can be very chaotic and challenging for a surrogate model to reproduce. Such chaotic values thus deteriorate the accuracy of any surrogate model, interpolation or ANN, and thus should not be used. Unfortunately, as of today there exists no simple or complete procedure to check the physical validity of a precomputed model of the Meudon PDR code. With a first scan of the datasets, we remove a few lines that are particularly affected. The total number of considered lines is therefore reduced from 5409 to *L* = 5375. This simple filter leaves other outliers in the training and test datasets. Though we observe that these outliers are rare – we expect less than 1%, we do not have any specific a priori knowledge on their location nor on their exact proportion. Manually checking the validity of each value is unrealistic given the sizes of the two datasets. The most informative hypothesis we can make on outliers is that if one line in a precomputed model is identified as an outlier, then it is likely for this precomputed model to contain other outliers,
especially in lines of the same species or of isotopologues. This hypothesis is exploited in the more thorough outlier detection method exploiting an ANN, introduced in Section 4.6.1.

**Summary of the setup** – The Meudon PDR code version to emulate is a function $\mathbf{f}: \theta \in \mathbb{R}^D \to y \in \mathbb{R}^L$, with $D = 4$ and $L = 5375$. As the integrated intensities predicted by the Meudon PDR code cover many decades, we train emulators of $\ln \mathbf{f}$, denoted $\ln \tilde{\mathbf{f}}$, to reproduce log-intensities $\ln y$. We assume the log-Meudon PDR code, $\ln \mathbf{f}$, to be twice continuously differentiable, except in the case of outliers that should be disregarded. This hypothesis $\ln \mathbf{f} \in \mathcal{C}^2$ is a reasonable physical assumption exploited in Chapter 5 (Section 5.2.1). In Section 4.6, we build emulators $\ln \tilde{\mathbf{f}}$ so that they satisfy this regularity property.

## 4.5.2 Comparison metrics

Interpolation methods and ANNs are compared on evaluation speed, memory requirements and approximation accuracy. We describe here the metrics used for the comparison, regardless of how the surrogate models are defined or trained.

The **evaluation speed** is measured on the full set of $L$ lines for 1000 random points. The measurements are performed on a personal laptop equipped with a 11th Gen Intel(R) Core(TM) i7-1185G7, with 8 logical cores running at 3.00 GHz. ANNs and interpolation methods are run on CPU to obtain a meaningful comparison. Running ANNs on a GPU could farther reduce their evaluation times. The implementations of interpolation methods are from the SciPy PYTHON package, popular in ISM studies (Wu et al., 2018). Nearest-neighbor, linear and RBF interpolation implementations allow for the evaluation of a vector function at once. Conversely, the spline interpolation implementation requires an explicit loop over the $L$ lines, which is slow in PYTHON. To avoid an unfair comparison with the other methods, the spline interpolation speeds are not evaluated.

The **memory requirements** are quantified with the number of parameters necessary to fully describe the surrogate model. Interpolation methods, for instance, require storing the full training set. It corresponds to $N_{train}(D+L) \approx 1.03 \times 10^8$ parameters. In PYTHON, these parameters are stored using 64-bits floating-point numbers. Storing the full grid requires about 1.65 GB.

The **accuracies** of surrogate models are evaluated on the test set, composed of points not used during training. To quantify accuracies, we define a new metric called the error factor (EF). This metric is a symmetric version of relative errors. Unlike the AE or the SE on $\ln \tilde{\mathbf{f}}$ and $\ln y$, this metric does not favor predictions lower than the true value. Indeed, as line intensities are considered in log-scale, the absolute error (Eq. 4.2) corresponds to the log-ratio of the predicted and true line intensities. The error factor is this ratio transformed back to linear scale. For a surrogate model $\tilde{\mathbf{f}}$ on a given tuple $(\theta, \ln y)$ and line $\ell$, it reads

$$ \mathrm{EF}(\ln \tilde{\mathbf{f}}; (\theta, \ln y_{\ell})) = \exp\left\{\left|\ln \tilde{f}_{\ell}(\theta) - \ln y_{\ell}\right|\right\} = \max\left\{\frac{\tilde{f}_{\ell}(\theta)}{y_{\ell}}, \frac{y_{\ell}}{\tilde{f}_{\ell}(\theta)}\right\}. \quad (4.6) $$

As the absolute value ensures positivity in log scale, an error factor is always superior or equal to 1. It can be expressed in percent using a $100 \times (EF - 1)$ transformation. For readability, the error factor is displayed in percent when $EF < 2$, i.e., 100%. An error factor that is not in percent is indicated by the multiplication sign, e.g., "$\times 3$" corresponds to $EF = 3$.

Figure 4.6 compares the relative error and the error factor. The error factor is a symmetrized relative error, as the absolute value also ensures symmetry in log scale. For small errors, i.e., $EF \approx 1$, it is similar to the standard relative error. However, for larger errors, the error factor is more relevant in our case. A standard relative error would return 100% for a factor of 2 too high and 50% for a factor of 2 too low, while in both cases $EF = 2$. In the worst case, a relative
error of 100% corresponds to a factor of 2 too high or a prediction of exactly 0, while EF = 2 in the former case and EF = +∞ in the latter. Minimizing a standard relative error would therefore lead to an under-estimation tendency, which is not the case for the proposed error factor.

<figure><img src="image_30.png" /><figcaption>Figure 4.6 – Comparison of the error factor and the relative error. To simplify the visualization, the errors are plotted in decimal log scale, while the surrogate models train on the natural log scale. We add 1 to the relative error to simplify the comparison and to avoid a divergence of the log-relative error towards -∞ when $\tilde{f}_{\ell}(\theta) \rightarrow y_{\ell}$. The relative error + 1 tends to 2 for large negative errors, i.e., the relative error tends to and is upper bounded by 100%. This effect favors underestimations compared to overestimations. Using the error factor cancels this undesirable effect.</figcaption></figure>

When applied to the full test set, the error factor yields a distribution of errors. This distribution is summarized by its mean, its 99th percentile and its maximum. The mean provides an estimation of the average error to expect. The 99th percentile and maximum provide upper bounds on the error. The maximum is very sensitive to outliers while the 99th percentile is more robust. Figure 4.7 illustrates this distribution summary on a fictional dataset of 20 000 points including 0.5% of outliers. The maximum is affected by the outliers, which induces a pessimistic bias for the corresponding error upper bound estimation. The 99th percentile is not significantly affected by the outliers, and provides a more relevant estimator of the actual upper bound of the error factor for this fictional dataset. This example shows that the choice of percentile is a trade-off based on the expected proportion of outliers. Lower percentiles such as 90th or 95th underestimate the upper bound on the error factor, and percentiles higher than 99.5th would in turn be sensitive to outliers like the max. The training and test sets generated with the Meudon PDR case are expected to contain less than 1% of outliers. The 99th percentile is therefore expected to be a relevant estimator of the error upper bound, robust to outliers.

In current IRAM-30m observations, the relative day-to-day calibration accuracy ranges from 3% to 10% – see e.g., Einig et al. (2023). The absolute flux calibration accuracy for ground based observations is more difficult to estimate but cannot be better than the relative calibration accuracy. For a surrogate model to be relevant for observations analysis and physical parameter inference, we set the constraint that satisfactory surrogate models must have a mean error factor below 10%.
<figure><img src="image_31.png" /><figcaption>Figure 4.7 – Illustration of the impact of outliers on the mean, the 99th percentile and the max estimators of the error factor. Note that the dataset used for this illustration is fictitious.</figcaption></figure>

## 4.6 Designing and training adapted ANNs

The choice of architecture and training approach of ANNs are now discussed. In the following, ANNs are trained with the MSE loss function. We set the activation function $g$ to the Gaussian error linear unit (GELU) (Hendrycks and Gimpel, 2023) to ensure that the resulting ANN is twice continuously differentiable. Unless explicitly mentioned, our ANNs have $H = 3$ hidden layers of equal size. This choice may not be optimal. A hyperparameter optimization step could improve the network performance, but would require a validation dataset and the training of many networks. As the results of Section 4.7 will show, this step is not necessary to obtain satisfactory results.

ISM models such as the Meudon PDR code present specificities, namely the presence of outliers and the unusual dimensions of the problem – very few inputs to predict many outputs. To address these specificities, dedicated strategies are needed and are described in the subsections below:

1. we apply an outlier removal procedure;

2. with a clustering method, we derive homogeneous line groups simpler to emulate with separate networks;

3. to select an adequate size for hidden layers, we resort to a dimension reduction method;

4. we apply a polynomial transform to augment the input data and thus ease the learning of nonlinearities;

5. we replace the standard ANN architecture by a dense architecture that exploits values in intermediate layers to re-use intermediate computations.

### 4.6.1 Removing outliers from the training set

Outliers that come from either numerical instabilities or physical bistabilities or multistabilities can be found in both the training and test sets, as described in Section 4.5.1. With a loss function such as the MSE, outliers in the training set greatly deteriorate the quality of a fitted ANN. A similar problem was reported in Auld et al. (2007) about the WMAP 3-yr likelihood code, and addressed by removing them from the dataset. In Auld et al. (2007), the outliers could be identified with a simple binary test. Such a test is not available in our case. More sophisticated methods are therefore required.

Performing regression in presence of outliers is a crucial topic in machine learning. Multiple methods exist for non-linear regression (Rousseeuw and Leroy, 1987). We resort to the method proposed in Motulsky and Brown (2006), which consists of 3 steps:

1. A statistical model is fitted to the training set with a strategy robust to outliers.
2. The training points with largest errors are reviewed.

3. Identified outliers are removed and a new statistical model is fitted to the cleaned training set.

To avoid any risk of biasing the analysis towards optimistic results, we do not remove any value from the test set. We now detail these three steps.

**Step 1: first fit robust to outliers** – For this first fit, we resort to an ANN designed as described at the introduction of Section 4.6. The size of hidden layers is fixed with the dimension reduction strategy that will be described in Section 4.6.2. We also include the polynomial transform of the input, to be described presented in Section 4.6.3. For specific outlier removal step, this fit is performed using the Cauchy loss function (CL). For a surrogate function $\tilde{f}$ and a pair $(\theta, \ln y)$, it reads

$$ \mathrm{CL}(\ln \tilde{f}; (\theta, \ln y_{\ell})) = \ln \left[ 1 + (\ln \tilde{f}_{\ell}(\theta) - \ln y_{\ell})^2 \right]. \quad (4.7) $$

Figure 4.8 shows how the squared error (Eq. 4.3), the absolute error (Eq. 4.2), the error factor (Eq. 4.6), and the Cauchy loss function penalize errors. The Cauchy function gives less weight to very large errors, which makes it more robust to outliers.

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

Figure 4.8 – Evolution of errors penalization by the considered loss functions. To simplify the visualization, the errors are plotted in decimal log scale, while the surrogate models train on the natural log scale. An error of 30 thus corresponds to a factor of $10^{30}$. Since some line intensities range from $10^{-50}$ to $10^{-2}$, very high error can occur early in the training phase.

**Step 2: review of training points with largest errors** – The review of training points with high errors is performed with a manual procedure. An instability in a given model of the grid does not affect all lines, as all lines are not emitted in the same spatial regions of the model. Therefore, we only remove affected lines instead of the full model. To accelerate this procedure, we exploit similarities between lines. For instance, when one water line intensity is identified as an outlier, it is highly likely that most of the water line intensities of the corresponding precomputed model are outliers. We emphasize that outliers are associated with instabilities or multistabilities. Physically consistent intensities that are challenging to reproduce, e.g., due to fast variations in a change of regimes, are not considered as outliers and maintained in the training set. In total, 71 239 values were identified as outliers, i.e., 0.069% of the training set. Note that this outlier
identification step is very informative, as it reveals regions of the parameter space that lead to multistabilities. However, studying these regions is beyond the scope of this thesis.

**Step 3: fit on the cleaned dataset** – Once outliers are identified, they are removed from the dataset, through a binary mask matrix $M = (m_{nl})_{nl}$. This mask permits to disregard only the identified outliers instead of removing all $L$ lines of precomputed models with at least one outlier. In this binary mask, $m_{nl} = 1$ indicates that $y_{nl}$ is an outlier and should not be taken into account, and $m_{nl} = 0$ indicates that $y_{nl}$ is not an outlier. Elements of the training set $(\theta_n, y_n) \in \mathbb{R}^D \times \mathbb{R}^L$ are augmented with corresponding binary mask vectors $m_n \in \{0, 1\}^L$. An ANN, on the one hand, can easily take this mask into account for training by computing the loss function and its gradient on non-masked values only. In the following, a masked version of the MSE relying on the binary mask $M$ is used when this outlier removal step is taken into account.

Existing implementations of interpolation methods, on the other hand, lack flexibility to handle such a mask during the fit. As some points of the grid are removed for some lines, the spline interpolation can not be applied on the masked training set. Nearest-neighbor, piecewise linear and RBF interpolation methods can be applied but would require line by line fits and predictions, as outliers don't occur for the same training points $\theta$ for all lines. Such line by line manipulation would be extremely slow with a PYTHON implementation. To present a somewhat meaningful comparison between ANNs and interpolation methods on the masked dataset, the masked values are imputed. This imputation relies on a line by line fit of an RBF interpolator with linear kernel. Masked values are replaced by interpolations computed from available data points. Interpolation methods are then fitted with this imputed training set.

Note also that it would be possible to avoid the manual review by training a network and automatically identifying outliers at once. Appendix 4.A outlines a possible approach as well as the reasons why we chose not to apply it in this particular case.

## 4.6.2 Exploiting correlations between line intensities

Line intensities computed by the Meudon PDR code come from the radiative de-excitation of energy levels. While non-local effects are accounted for in the radiative transfer, the excitation of many lines is affected to a large extent by local variables such as the gas temperature or density. As a result, high correlations between some lines are expected. Figure 4.9 shows the $L \times L$ matrix of absolute Pearson correlation coefficients, with lines grouped by molecule. We indeed find some strong correlations. In particular, lines from the same species are often highly correlated, especially for water isotopologues and molecular hydrogen. However, some species produce lines that are not correlated. For instance, high energy lines of SO have a very small correlation with low energy lines, as the corresponding sub-matrix has a diagonal shape. Finally, some lines from different species are highly correlated, e.g., $OH^+$, $SH^+$ and $H_2$. Handling the $L$ lines independently, as in interpolation methods, ignores these correlations in the line intensities. We exploit these correlations with two strategies: line clustering and dimension reduction.

### Line clustering to divide and conquer

Figure 4.9 highlights some clusters of highly correlated lines. These clusters are not simply related to the line carrier. We derive clusters of lines automatically from the correlation matrix using the spectral clustering algorithm (Shalev-Shwartz and Ben-David, 2014, chapter 22). Spectral clustering defines clusters such that lines from a same cluster are as similar as possible and such that lines from different clusters are as different as possible. It relies on a similarity measure such as the Pearson correlation, while most clustering algorithms are distance-based. We set the number of clusters to the value that maximizes the ratio of intra to inter-cluster mean correlations.

Figure 4.10 presents the 4 clusters we obtain. The mean intra and inter-cluster correlations are 0.895 et 0.462 respectively, while the mean correlation among all lines is 0.73. The obtained clusters contain 3 712, 1 272, 241 and 150 lines, respectively. This imbalance between clusters
<figure><img src="image_33.png" /><figcaption>Figure 4.9 – Meudon PDR code correlations among the *L* = 5 375 predicted lines from 27 chemical species, shown with the *L* × *L* matrix of absolute Pearson correlation coefficients. A value of exactly 1 for two lines means that there exists an exact affine relationship between their log-intensities. The black squares on the diagonal group lines from a common chemical species. For readability reasons, only the names of species with more than 100 lines predicted by the Meudon PDR code are displayed.</figcaption></figure>
comes from the imbalance between molecules. For instance, H₂ corresponds to 3282 lines, i.e., 61% of the lines computed by the Meudon PDR code. Figure 4.9 shows that all these lines are highly correlated. Appendix 4.B provides a more complete description of the content of these 4 clusters. With this approach, an ANN is trained for each cluster.

### Using PCA to set the size of the last hidden layer

A second and complementary approach to exploit these correlations is to assume that a vector y with the *L* line intensities can be compressed to a vector of size $\tilde{L} < L$ with limited loss of information. Formally, we assume that the line intensities y live in a subspace of dimension $\tilde{L} < L$ where $\tilde{L}$ can be estimated using a dimension reduction algorithm. We resort to a principal component analysis (PCA) (Shalev-Shwartz and Ben-David, 2014, chapter 23) on the training set, which performs compressions using only affine transformations. We obtain that compressing all $L = 5375$ lines with only $\tilde{L} \approx 1000$ principal components and then decompressing leads to a mean error factor below 0.1% on the training set, which confirms our hypothesis.

Figure 4.1 shows that in an ANN such that $D \ll L$, most parameters belong to the last hidden layer. The size of this layer is thus critical to obtain a good accuracy. Too large a layer might lead to overfitting; too small a layer could not capture the nonlinearities of the dataset. Using PCA, we determine the dimension $\tilde{L} < L$ of a subspace in which the code outputs can be described with limited accuracy loss. In regression, this last hidden layer also applies an affine transformation. We therefore set the size of the last hidden layer to the estimated dimension $\tilde{L}$. To predict the $\tilde{L}$ intermediate values of the last hidden layer, which are then used to predict all $L$ line intensities, the first two hidden layers are set with the same size.

**Link with principal component regression** – This approach shares similarities with principal component regression, used e.g., in Spurio Mancini et al. (2022), but is not equivalent. Using PCA, one obtains the optimal affine transformation to a subspace of dimension $\tilde{L} \ll L$. The proposed approach uses PCA to estimate $\tilde{L}$ to set the size for the hidden layers to a relevant value. Principal component regression (PCR) goes farther by fixing the output layer of the network to the learned transformation. This approach might seem very advantageous in the considered case. Since the output layer is the largest layer in our networks, it would reduce the number of parameters to learn. In this work, we tried both approaches. We chose to leave the output layer as a learnable parameter as PCR leads to worse results than those presented in the manuscript.

**Application to the four clusters of lines** – For the networks trained on the four clusters of lines obtained in Section 4.6.2, the size of the last hidden layer is also set to the minimum number of principal components that ensures a decompression with mean error factor below 0.1% on the training set. The obtained sizes $\tilde{L}$ are approximately 500 (about 13% of the $L = 3712$ lines of the cluster), 350 (about 28%), 100 (about 41%) and 75 (50%), respectively. As the bigger clusters are the most homogeneous, they have the smallest ratio $\tilde{L}/L$ of subspace dimension $\tilde{L}$ with the total number of lines $L$. The number of parameters necessary to describe four small specialized ANNs is thus greatly reduced in comparison to a single larger general network.

### 4.6.3 A polynomial transform to learn nonlinearities

The nonlinearities in the Meudon PDR code make the approximation task challenging. In an ANN, nonlinearities come from the activation function *g*. However, learning meaningful and diversified nonlinearities is difficult with few hidden layers. Conversely, an ANN with numerous layers can lead to overfitting and requires more time for evaluations and memory for storage. Preprocessing the physical parameters *θ* with a variety of predefined non-linear functions eases this learning task while maintaining a small network architecture. We choose to apply a polynomial transform $P_p$ which replaces the input vector *θ* of dimension *D* with an input vector containing all monomials computed from the *D* entries of degree up to *p*. For instance, for $D = 3$ and
<figure><img src="image_34.png" /></figure>

Figure 4.10 – On both top and bottom plots: each column corresponds to one of the four line clusters obtained in Section 4.6.2.

Top row: composition of each cluster. Each bar indicates the proportion of lines of a species in a cluster. The red crosses correspond to exactly zero line. For each cluster, the 3 species with the most lines are highlighted.

Bottom row: Pearson correlation of the most representative line of each cluster with the 3 main physical parameters. The most representative line of a cluster is defined as the line with the highest average correlation with the other lines. A round marker at a vertex indicates a negative correlation.
$p = 2, \theta = (\theta_1, \theta_2, \theta_3)$ is replaced with $P_2(\theta) = (\theta_1, \theta_2, \theta_3, \theta_1^2, \theta_2^2, \theta_3^2, \theta_1\theta_2, \theta_1\theta_3, \theta_2\theta_3) \in \mathbb{R}^9$.
For $D = 4$ and $p = 3$, $P_3(\theta) \in \mathbb{R}^{34}$. This approach is classic in regression (Ostertagová, 2012),
but less common to train ANNs.

It is well known in polynomial regression that a high maximum degree $p$ can lead to overfitting (Shalev-Shwartz and Ben-David, 2014, chapter 11). The analysis of the physical processes indicates that the gas structure and emission properties depend on control quantities combining $G_0, n_H$ (or $P_{\text{th}}$) and $A_V^{\text{tot}}$. For instance, $G_0/n_H$ is known to play an important role in PDRs (Sternberg et al., 2014). It is therefore important to consider monomials combining these 3 physical parameters. In contrast, the angle $\varphi$ is assumed to have a simpler role in the model. To avoid overfitting, we choose the minimum value that combines the three parameters, $p = 3$, and thus consider the polynomial transforms $P_3$. This transformation is applied to the input variables after the preprocessing step described in Section 4.5.1 – log scale for $P_{\text{th}}$, $G_0$ and $A_V^{\text{tot}}$, and standardization of the $D = 4$ parameters. It is implemented as an additional first fixed hidden layer. The gradient $\nabla_{\theta} \tilde{f}$ can thus be efficiently evaluated with auto-differentiation methods.

### 4.6.4 Dense networks to reuse intermediate computations

**Introduction to the dense architecture** – The fully connected ANNs architecture considered so far, shown in Figure 4.1, is widely used in the deep learning community. However, this architecture struggles to maintain input information in hidden layers, as it is transformed with non-linear activation functions. It might therefore fail to reproduce very simple relationships. For instance, the intensity of UV-pumped lines of H₂ is highly correlated with $G_0$. Using $G_0$ directly to predict intensities of such lines thus might be more relevant than passing it through non-linear transformations. This architecture also struggles to pass gradient information all the way back to the first hidden layers. This phenomenon, called *gradient vanishing*, might lead to largely sub-optimal trained networks. The recent residual (He et al., 2016) and dense architectures (Huang et al., 2017) address these two issues. We use the dense architecture for our regression problem.

A dense architecture is a special type of feedforward architecture where the input of a layer $j+1$ is the concatenation of the input and output vectors of the previous layer $j$: $\theta^{(j+1)} = [\theta^{(j)}, y^{(j)}]$. This architecture focuses on reusing intermediate values in hidden layers and can thus reduce the number of parameters to train.

Figure 4.11 illustrates this dense architecture for a simple ANN with $H = 2$ hidden layers and the same sequence of layer input sizes ($i_j$)$_{j=1}^{H+1}$ used to illustrate the standard feedforward architecture in Figure 4.1. The output sizes $o_j$ of hidden layers are much smaller with the dense architecture, as the input of layer $j$ concatenates the input and output of layer $j-1$. The weight matrices $W^{(j)}$ of hidden layers are thus much smaller as well, which reduces the total number of parameters to train. By lowering the number of parameters to learn while providing the same number of inputs to the output layer, this architecture limits overfitting risks.

**Considered dense architecture** – As the number of parameters per layer is reduced, we define ANNs with $H = 9$ hidden layers, i.e., 6 more layers than for the proposed networks with the standard architecture, and yet with a similar total number of parameters. By definition, the size of the hidden layers in a dense architecture is strictly increasing, as the size $i_{j+1}$ of a layer input is the sum $i_j + o_j$ of the input and output sizes of the previous layer. The network is set so that the input $i_{j+1}$ of a layer $j+1$ is 50% larger than the input of the previous layer $i_j$. With this geometric progression and the polynomial transform $P_3$, the input of the output layer contains 1296 neurons, which is 29.6% larger than the recommendation from PCA obtained in Section 4.6.2. However, out of these 1296 neurons, 34 correspond to the input values, 17 to the output of the first hidden layer, 25 to the output of the second hidden layer, etc. In other words, though the input of the output layer contains more neurons for the considered dense neural network than the PCA recommendation, a majority of these neurons are the result of fewer transformations.

When using this dense architecture strategy with the clustering approach, four dense networks
<figure><img src="image_35.png" /><figcaption>Figure 4.11 – Structure of a dense neural network with $H = 2$ hidden layers and the same sequence of layer input sizes ($i_j$)$_{j=1}^{H+1}$ used to illustrate the feedforward architecture in Figure 4.1.</figcaption></figure>

with $H = 9$ hidden layers are designed. The size of the last hidden layer is also set to a slightly larger value than the PCA recommendation. The geometric progressions of these 4 networks are set to 35%, 30%, 15% and 10%, respectively.

## 4.7 Experiments: application to the Meudon PDR code

The ANNs designed and trained with the proposed strategies are compared with interpolation methods with respect to accuracy, memory and speed. Table 4.3 shows the results of the comparison. It is divided in two halves. The first presents models trained on the raw training set, while the second contains models trained on the cleaned training set. The cleaned dataset is based on the mask defined with the outlier detection procedure of Section 4.6.1. In each half, the results of interpolation methods are first listed, followed by ANNs combining one or more of the presented strategies.

The code and data produced for this experiment are public. The training set, the test sets, and the mask on the train set can be found at https://ism.obspm.fr/files/ArticleData/2023_Palud_Einig/2023_Palud_Einig_trained_ANN.zip. The code can be found at https://github.com/einigl/ism-model-nn-approximation. All proposed ANNs were implemented using the PYTORCH 4 PYTHON library (Paszke et al., 2019). The most accurate ANN is available at https://ism.obspm.fr/files/ArticleData/2023_Palud_Einig/2023_Palud_Einig_data.zip.

### 4.7.1 Performance analysis

The proposed ANNs outperform all interpolation methods on all aspects by a large margin: they are between 100 and 1000 times faster than reasonably accurate interpolation methods, and between 14 and 38 times lighter in terms of memory. Interpolation methods handle the prediction of $L$ lines as $L$ independent operations, while ANNs handle the $L$ lines at once, which is much faster. Interpolation methods require storing the full training set that contains 103 million 64-bit floating point numbers, i.e., 1.65 GB in size. In contrast, ANNs use shared intermediate values in hidden layers to predict all lines, which limits redundant computations and effectively compresses the dataset. They can thus be fully described with between 2.7 and 7.8 million parameters, i.e., between 43 MB and 118 MB. Finally, the proposed ANNs are roughly twice as accurate as the best interpolation methods on average and between two and three times as accurate with respect to

<footnote n="4">4. https://pytorch.org/</footnote>
Table 4.3 – Performance of interpolation methods and of the proposed ANNs with and without removal of outlier from the training set. Evaluation speeds are measured on the full set of *L* lines for 1 000 random points. The measurements are performed on a personal laptop equipped with 8 logical cores running at 3.00 GHz. Error factors are evaluated on the test set. Spline interpolation methods lack vectorized implementation in Python and are thus slow: their speed is not evaluated for comparison fairness. For ANN architectures:

- R: design of the last hidden layer using PCA,

- P: polynomial transform of the input up to degree 3,

- C: line clustering and parallel training of four specialist ANNs,

- D: dense architecture.

<table><thead><tr><th rowspan="2"></th><th rowspan="2">Method</th><th colspan="3">Error factor</th><th>Memory</th><th>Speed</th></tr><tr><th>mean</th><th>99th per.</th><th>max</th><th>(MB)</th><th>(ms)</th></tr></thead><tbody><tr><td rowspan="10">No outlier removal</td><td rowspan="2">near. neighbor<br/>linear</td><td>×13.1</td><td>×11.3</td><td>×3e5</td><td>1650</td><td>62</td></tr><tr><td>15.7</td><td>×2.3</td><td>×143</td><td>1650</td><td>1.5e3</td></tr><tr><td rowspan="3">spline<br/>linear<br/>cubic<br/>quintic</td><td>15.7</td><td>×2.3</td><td>×144</td><td>1650</td><td>...</td></tr><tr><td>11.2</td><td>×2.2</td><td>×122</td><td>1650</td><td>...</td></tr><tr><td>19.1</td><td>×2.9</td><td>×304</td><td>1650</td><td>...</td></tr><tr><td rowspan="3">RBF<br/>linear<br/>cubic<br/>quintic</td><td>10.2</td><td>96.8</td><td>×99</td><td>1650</td><td>1.1e4</td></tr><tr><td>10.4</td><td>×2.1</td><td>×112</td><td>1650</td><td>1.1e4</td></tr><tr><td>10.9</td><td>×2.1</td><td>×118</td><td>1650</td><td>1.1e4</td></tr><tr><td rowspan="2">ANN<br/>R<br/>R+P</td><td>7.3</td><td>64.8</td><td>×81</td><td>118</td><td>12</td></tr><tr><td><b>6.2</b></td><td><b>49.7</b></td><td>×84</td><td>118</td><td>13</td></tr><tr><td rowspan="10">Outlier removal on training set</td><td rowspan="2">near. neighbor<br/>linear</td><td>×13.1</td><td>×11.6</td><td>×3e5</td><td>1650</td><td>62</td></tr><tr><td>15.9</td><td>×2.4</td><td>×143</td><td>1650</td><td>1.5e3</td></tr><tr><td rowspan="3">spline<br/>linear<br/>cubic<br/>quintic</td><td>15.9</td><td>×2.4</td><td>×144</td><td>1650</td><td>...</td></tr><tr><td>11.1</td><td>×2.2</td><td>×120</td><td>1650</td><td>...</td></tr><tr><td>20.0</td><td>×2.7</td><td>×285</td><td>1650</td><td>...</td></tr><tr><td rowspan="3">RBF<br/>linear<br/>cubic<br/>quintic</td><td>10.3</td><td>97.3</td><td>×97.5</td><td>1650</td><td>1.1e4</td></tr><tr><td>10.5</td><td>×2.0</td><td>×106</td><td>1650</td><td>1.1e4</td></tr><tr><td>10.9</td><td>×2.0</td><td>×114</td><td>1650</td><td>1.1e4</td></tr><tr><td rowspan="2">ANN<br/>R<br/>R+P</td><td>5.1</td><td>42.0</td><td>×32.8</td><td>118</td><td>12</td></tr><tr><td>5.5</td><td>42.3</td><td>×41</td><td>118</td><td>13</td></tr><tr><td rowspan="3">ANN<br/>R+P+C<br/>R+P+D<br/>R+P+C+D</td><td>4.9</td><td>44.5</td><td>×44</td><td>51</td><td>14</td></tr><tr><td><b>4.5</b></td><td><b>33.1</b></td><td>×33.8</td><td>125</td><td><b>11</b></td></tr><tr><td>4.8</td><td>37.9</td><td>×37.6</td><td><b>43</b></td><td>14</td></tr></tbody></table>