the percentile 99th. Overall, the proposed ANNs are the only surrogate models that yield a mean error factor lower than 10% and thus that are suited to a comparison with actual observations.

## 4.7.2 Removing outliers is crucial

When the outlier removal step is not applied, the distribution of errors is highly skewed for all surrogate models. For interpolation methods, the 99th percentile reveals that around 99% of the predictions correspond to errors lower than a factor of 2. For the best neural network (R+P), it reveals that 99% of the errors are lower than 49.7%. However, for all methods, the maximum error is at least 80 times higher than the 99th percentile and reaches unacceptable values. An inspection of the highest errors reveals that they are close to training points with outliers, which indicates that these outliers significantly deteriorate the fit.

After removing outliers from the training set, interpolation methods do not show average accuracy improvement. Only little improvements can be observed on the 99th percentile and maximum EF, especially for the RBF interpolation methods. Replacing outliers with interpolated values is therefore not relevant to derive surrogate models based on interpolation methods in this case. In contrast, the two ANNs trained both with and without the outlier removal step (R and R+P) show consequent improvements. With outlier removal, the mean EF decreased from 7.3% and 6.2% to 5.1% and 5.5%, respectively. Similarly, the 99th percentile dropped from 64.8% and 49.7% to 42% and 42.3%. Finally, the maximum error is reduced by more than a factor of 2. These important improvements demonstrate the interest of filtering outliers from the training set before training ANNs.

## 4.7.3 The importance of the polynomial feature augmentation

The polynomial transform improves the accuracy in presence of outliers in the training set, but deteriorates it after the outlier removal step. It provides flexibility to learn abrupt nonlinearities caused by outliers. However, with outliers removed, the function to learn is smoother. The EF on the masked training set is 1.44% without the polynomial transform (R) and 0.77% with it (R+P), while the EF on the test set is lower without the polynomial transform. This improvement on the training set does not lead to an improvement on the test set, suggesting an overfit. The polynomial transform therefore requires additional strategies to better reproduce data unused during the training phase.

Both the clustering step and dense architecture, used with the polynomial transform, led to better accuracy. The surrogate model that exploits the line clustering but not the dense architecture (R+P+C) improves the mean accuracy by 0.2 percentage points, while requiring 57% fewer parameters than the first two networks (R and R+P). A potential cause of the average error factor improvement is the separation of the trainings of each specialized ANN. Since $H_2$ lines represent 61% of all $L$ lines, they dominate the loss function and thus are learned in priority. To separate them from other clusters might have improved performance on those other clusters.

The surrogate model based on a single network with dense architecture (R+P+D) is the most accurate on average and provides the lowest error upper bound for the robust 99th percentile estimator. Even with more trainable parameters than the first two networks, it does not overfit. It is also the fastest model as reusing intermediate values reduces the number of computations.

Finally, combining both line clustering and dense architectures (R+P+C+D) yields the lowest memory usage with only 2.7 million parameters, i.e., 43.2 MB, which is 38 times lighter than for interpolation methods. It also provides very good accuracy, both on average and for upper bounds.

Overall, a dense architecture and the line clustering effectively limit overfitting and thus perform better on data not used during the training phase. The line clustering leads to the lightest models regarding memory requirements, and the dense architecture to the most accurate models.
## 4.8 Conclusion

The interpretation of observations of atomic and molecular tracers in galactic and extragalactic interstellar medium requires comparison with state-of-the-art astrophysical models to infer physical conditions. Such inference procedure requires numerous evaluations of the numerical model, particularly so for Bayesian approaches. Inference on large observation maps – that are becoming more and more common – farther relies on many evaluations. ISM models are often too slow to perform such inference and are generally emulated using grids of precomputed models. This emulation approach induces errors that are seldom quantified in the literature. Besides, these methods can have high evaluation time and memory costs.

In this chapter, the general problem of deriving a fast, accurate and memory-light surrogate model for a time-consuming ISM numerical model has been addressed. The proposed approach has been assessed in the case of the Meudon PDR code, a state-of-the-art ISM code. Four common families of interpolation methods – nearest-neighbor, linear, spline and RBF – have been compared to specifically designed ANNs. We found that ANNs outperform all interpolation methods by a large margin in terms of accuracy, speed and memory usage.

Attaining this performance level for an ISM model such as the Meudon PDR code requires addressing their specificities. First, ISM models usually predict many statistically independent observables – e.g., line intensities of many species – from few parameters – e.g., gas density or temperature. This setting is unusual in ANN applications – except for ANNs that generate structured data such as images, text or time series. Second, due to numerical instabilities or physical bistabilities or multistabilities, such models sometimes produce outliers that harm the training process. In this chapter, we proposed and combined five strategies to design and train adapted ANNs:

- To identify outliers, we train a first ANN with a loss function robust to large errors. Training points corresponding to large errors are manually reviewed. Identified outliers are removed from the training set.

- Lines are clustered into homogeneous subsets that are simpler to emulate: for each cluster one ANN is defined and trained.

- A dimension reduction technique, PCA, is used to determine an adequate size of hidden layers.

- A polynomial transform of the input physical parameters provides precomputed nonlinearities to the network, which permits the learning of nonlinearities with a limited number of hidden layers.

- A dense architecture exploits intermediate computations and thus limits redundant computations. Using such an architecture instead of the standard feedforward neural network architecture improves speed and avoids overfitting.

With the proposed strategies, ANNs achieve 4.5% average accuracy while the best interpolation method, RBF, attains 10.2%. Upper bound on the errors, quantified using their 99th percentile, reaches 33.1% for our ANNs compared to 97% for RBF. Besides, our ANNs are 1000 times faster than RBF and are more than 10 times lighter in terms of memory. This chapter focuses on an application to the Meudon PDR code, motivated by the inverse problem at the core of this thesis. However, the proposed strategies are sufficiently general to be applicable to many other ISM models.

As we will show in Chapter 5, the fast and accurate ANN emulators obtained in this chapter enables performing sampling-based inference on observation maps using the Meudon PDR code, a physically comprehensive model. As we will show in the Chapter 6, it will also permit efficient analyses of large observations maps produced by today's instruments such as the JWST, ALMA, or the ORION-B dataset observed by the IRAM 30m (Pety et al., 2017).
## Appendix 4.A Automatic outlier detection procedure

In this Section, we describe one possible automatic method for outlier identification along with the reasons why we chose not to use it. This method consists in simultaneously identifying outliers and training the ANN with the corresponding masked loss function. A similar approach was used in Gratier et al. (2016), where the likelihood combined two generative models: one for outliers values and one for non-outliers values. An additional parameter – to be tuned during training – determines how likely a value is to be an outlier. Formally, the mask M would be considered as a continuous variable in $[0, 1]^{N \times L}$ optimized along with the parameters $\psi$. Using a masked squared error function, one loss function option could be

$$ \mathcal{L}(\tilde{\mathbf{f}}, \mathbf{M}; \mathcal{D}) = \left[ \frac{1}{NL} \sum_{n=1}^{N} \sum_{\ell=1}^{L} (1 - m_{n\ell})(\ln \tilde{f}_{\ell}(\theta_n) - \ln y_{n\ell})^2 \right] + \lambda_1 r_1(\mathbf{M}) + \lambda_2 r_2(\mathbf{M}), \quad (4.8) $$

where $r_1$ is a regularization function that favors $m_{n\ell}$ values close to 0 or 1, where $r_2$ is a regularization function that favors masks M that agree with a priori knowledge on outliers, and where $\lambda_1, \lambda_2 > 0$ are regularization weights to be tuned. These properties are satisfied for instance with

$$ r_1(\mathbf{M}) = \iota_{[0,1]}(m_{n\ell}) + \frac{1}{NL} m_{n\ell}(1 - m_{n\ell}) \quad (4.9) $$

where $\iota_{[0,1]}$ is the indicator function on the $[0, 1]$ interval, i.e., where $\iota_{[0,1]}(x) = 0$ if $x \in [0, 1]$ and $+\infty$ otherwise, and

$$ r_2(\mathbf{M}) = \|\mathbf{M}\|_1 = \sum_{n=1}^{N} \sum_{\ell=1}^{L} |m_{n\ell}|. \quad (4.10) $$

Figure 4.12 shows the corresponding loss function for one value $m_{n\ell}$. The total regularization has a global minimum of value 0 at $m_{n\ell} = 0$, i.e., when $y_{n\ell}$ is not considered as an outlier. It also has a local minimum at $m_{n\ell} = 1$, i.e., when $y_{n\ell}$ is considered as an outlier. This local minimum has a higher value to limit the number of identified outliers. Higher values of $\lambda_2$ result in higher values attained at this local minimum, which in turn result in less identified outliers. Intermediate values are penalized to encourage the mask values towards 0 and 1.

The main drawback of this approach is that the model might classify important and physically consistent points that are challenging to emulate as outliers. The presented function $r_2$ considers all couples $(n, \ell)$ independently, which would cause any minimization algorithm to mask only the hardest values to reproduce with an ANN. These hardest values to reproduce might not be outliers and carry significant physical content that a relevant surrogate model must learn. For instance, values close to an abrupt change of regime may be hard to reproduce, but an emulator that disregards them would be irrelevant in observation analyses. More subtle regularization functions $r_2$ incorporating a priori knowledge on the outliers should thus be designed. For instance, $r_2$ could encode the hypothesis stated in Section 4.5.1, that is, if a line of a precomputed model is identified as an outlier, then it is likely that this precomputed model contains other outliers, especially among the lines emitted by the same species and its isotopologues. The design of such regularization functions is not trivial, and will then require tuning the associated parameters $\lambda_1$ and $\lambda_2$, which is difficult as well. For this reason, we chose to use the method described in Section 4.6.1 for our case, with its required manual review of a part of the training set. In cases where informative prior knowledge on outlier location or distribution is accessible, and where the mathematical formulation of such prior knowledge is manageable, this type of approach can bypass the need for a manual review.
<figure><img src="image_36.png" /></figure>

Figure 4.12 – Regularization on the mask M with the $r_1$ and $r_2$ functions proposed in Eq. 4.9 and Eq. 4.10, respectively, and $\lambda_1 = 1$. Values out of the [0, 1] interval are +∞ due to the indicator function.

## Appendix 4.B Content of clusters of lines

Figure 4.10 presents the four clusters obtained in Section 4.6.2. We now describe the line content of the four clusters. All species have lines distributed in at most 3 clusters, except for $^{12}$CO that has lines in each of the 4 clusters. CO lines are indexed with 2 quantum numbers: the rotational number *J* and the vibrational number *v*.

**Cluster 1** – The first cluster gathers lines that are emitted from the most external UV illuminated layers of the cloud and trace hot chemistry. It includes all H$_2$ lines but three, and is thus the largest. It also contains all lines from OH$^+$ (209), CH$^+$ (10), all 8 lines of S$^+$, all 5 lines of N, all 3 lines of O and the 158 µm line of C$^+$, and CO transitions with low *J* values for *v* = 1–0 and *v* = 1–1. The line intensities of this cluster are highly and positively correlated to *G*$_0$, and not correlated at all with $A_V^{\text{tot}}$.

**Cluster 2** – The second cluster contains 99% of the 655 lines from water H$_2$$^{16}$O and its isotopologue H$_2$$^{18}$O. It also contains lines from high energy levels for several molecules (HD, CO, $^{13}$CO, C$^{18}$O, $^{13}$C$^{18}$O, HNC, HCN, HCO$^+$, SO, CN, SH$^+$, C$_2$H, OH, and CS), and transitions from the neutral atoms C, Si and S. Line intensities in this cluster are positively correlated with $P_{\text{th}}$ and $G_0$ and not at all with $A_V^{\text{tot}}$.

**Cluster 3** – The third cluster contains mostly c-C$_3$H$_2$ lines, and some C$_2$H lines. It also includes two transitions of CO with moderate *J* values in the lowest vibrational level *v* = 0 (*J* = 3–2 and *J* = 4–3). Its line intensities are overall positively correlated with $P_{\text{th}}$ and $A_V^{\text{tot}}$ and negatively correlated with $G_0$.

**Cluster 4** – The fourth cluster contains the low energy lines of $^{13}$CO, C$^{18}$O, $^{13}$C$^{18}$O, HNC, HCN, HCO$^+$, SO and c-C$_3$H$_2$, and the lowest temperature transitions of CO (*J* = 1–0 and *J* = 2–1). Its line intensities are overall positively correlated with $P_{\text{th}}$, strongly positively correlated with $A_V^{\text{tot}}$ and negatively correlated with $G_0$.
# Chapter 5

# An MCMC algorithm for efficient inversion with quantified uncertainty

> "At first glance, it might appear surprising that a trivial mathematical result obtained by an obscure minister over 200 years ago ought still to excite so much interest across so many disciplines, from econometrics to biostatistics, from financial risk analysis to cosmology."

Trotta (2008)

## Contents

<table><tr><td>5.1 Proposed statistical model</td><td>114</td></tr><tr><td>5.1.1 Approximation of the likelihood function</td><td>114</td></tr><tr><td>5.1.2 Prior distribution</td><td>117</td></tr><tr><td>5.1.3 Posterior distribution</td><td>118</td></tr><tr><td>5.2 Proposed MCMC algorithm</td><td>118</td></tr><tr><td>5.2.1 PMALA transition kernel</td><td>118</td></tr><tr><td>5.2.2 MTM transition kernel</td><td>120</td></tr><tr><td>5.2.3 Proposed sampler and implementation details</td><td>121</td></tr><tr><td>5.2.4 Illustration: 2D Gaussian mixture model</td><td>123</td></tr><tr><td>5.3 Model checking using a predictive posterior p-value</td><td>126</td></tr><tr><td>5.3.1 Definition of selected p-value</td><td>126</td></tr><tr><td>5.3.2 p-value approximation and associated uncertainties</td><td>126</td></tr><tr><td>5.3.3 Illustration: 2D Gaussian distribution</td><td>128</td></tr><tr><td>5.4 Applications</td><td>130</td></tr><tr><td>5.4.1 Sensor localization</td><td>130</td></tr><tr><td>5.4.2 Realistic astrophysical data</td><td>132</td></tr><tr><td>5.5 Conclusion</td><td>136</td></tr><tr><td>Appendix 5.A Tuning automatically the prior hyperparameters</td><td>137</td></tr><tr><td>Appendix 5.B Optimization of the approximation parameter</td><td>137</td></tr><tr><td>Appendix 5.C Sampling from the smoothed indicator distribution</td><td>138</td></tr></table>

In Part I, we depicted the interstellar medium (ISM), presented the Bayesian approach to model and solve an inverse problem, and reviewed many statistical inference applications in ISM
studies. In Chapter 4, we derived a fast, light and accurate artificial neural network (ANN) emulator of the Meudon PDR code, i.e., of the forward model considered in this thesis inverse problem. By construction, this ANN is twice continuously differentiable.

As outlined in Chapter 3, the proposed observation model combines additive Gaussian noise, multiplicative lognormal calibration errors, and censorship. In this chapter, we formally set the inverse problem studied in this thesis. Our goal is to derive large maps of physical parameters $\Theta \in \mathbb{R}^{N \times D}$ from large maps of observations $\mathbf{Y} \in \mathbb{R}^{N \times L}$, with $N$ the number of beams or pixels in the maps, $D$ the number of considered physical parameters (e.g., the thermal pressure or the visual extinction) and $L$ the number of observed atomic or molecular integrated emission lines. This inverse problem is solved with a new MCMC sampler. We resort to an MCMC algorithm to provide credibility intervals associated with the estimations. This quantifies the uncertainty for applications such as astrophysics where no ground truth is available. We exploit the map structure and resort to a spatial regularization prior. As the resulting posterior distribution is potentially multimodal, we design a dedicated MCMC sampler to explore it. This sampler combines two sampling kernels: one responsible for global exploration, and one for efficient local exploration.

This chapter is based on Palud et al. (2023b), a journal article published in IEEE transactions on signal processing (TSP) that describes the model and the sampler. Two conference articles also introduce the model and the sampler: Palud et al. (2022a), published in the Gretsi conference, and Palud et al. (2022b). Finally, the Bayesian hypothesis testing was published in a Gretsi conference article, Palud et al. (2023a).

The code associated with this chapter and its experiments was published in a PYTHON package called BEETROOTS (for BayEsian infErence with spaTial Regularization of nOisy multi-line Ob-servaTion mapS). The code can be found at https://github.com/pierrePalud/beetroots. It can be installed via PIP – see https://pypi.org/project/beetroots/. The associated doc-umentation can be found at https://beetroots.readthedocs.io/en/latest/index.html.

In Section 5.1, we derive an approximation of the likelihood function with controlled error. The prior and posterior distributions are also introduced. Section 5.2 introduces the proposed sampler. In the illustrative Section 5.4, the proposed sampler is applied to two classical synthetic multimodal use cases, namely a Gaussian mixture and a sensor localization problem. Results are compared with other state-of-the-art samplers outlined in Chapter 2 (Section 2.A). In the illustrative Section 5.4.2, the sampler is validated on the astrophysical inverse problem of interest on synthetic data.

# 5.1 Proposed statistical model

This section first proposes an approximation of the likelihood with controlled error based on a combination of purely Gaussian additive and purely lognormal multiplicative approximations. Then, it introduces the considered prior distribution. Finally, we combine the prior distribution with the likelihood function to define the posterior distribution.

## 5.1.1 Approximation of the likelihood function

**Summary of the observation model** – The observation model presented in Chapter 3 (Section 3.4) is briefly recalled. Individual observations $\mathbf{y} = (y_\ell)_{\ell=1}^L$ gather $L$ lines. They are considered to be generated from some parameter $\boldsymbol{\theta} \in \mathbb{R}^D$ and the Meudon PDR code $\mathbf{f}: \mathbb{R}^D \to \mathbb{R}^L$, where $D$ is assumed to remain moderate, e.g. $D \lesssim 10$. The Meudon PDR code prediction for a channel $\ell$ is denoted by $f_\ell$, so that for any $\boldsymbol{\theta} \in \mathbb{R}^D$, $\mathbf{f}(\boldsymbol{\theta}) = (f_1(\boldsymbol{\theta}), \cdots, f_L(\boldsymbol{\theta}))$. As a single evaluation of the Meudon PDR code $\mathbf{f}$ is slow and computationally expensive, it is replaced by the neural network-based surrogate model $\tilde{\mathbf{f}}$ derived in Chapter 4. Individual observations and parameters are grouped in maps $\mathbf{Y} = (y_n)_{n=1}^N$ and $\boldsymbol{\Theta} = (\theta_n)_{n=1}^N$ of $N$ pixels, with $N$ potentially $O(10^4)$. The sensors are assumed to have a lower limit of sensitivity $\omega_{nl} \in \mathbb{R}$ below which
an observation is considered censored. Both an additive and multiplicative noise degrade the observations. The resulting observation model is, for $n \in [1, N]$ and $\ell \in [1, L]$,

$$y_{nl} = \max \left\{ \omega_{nl}, \varepsilon_{nl}^{(m)} \tilde{f}_{\ell}(\theta_n) + \varepsilon_{nl}^{(a)} \right\}, \quad (5.1)$$

where $\varepsilon_{nl}^{(a)} \sim N(0, \sigma_{a,nl}^2)$ is an additive Gaussian white noise, and $\varepsilon_{nl}^{(m)} \sim \text{Lognormal}(-\sigma_m^2/2, \sigma_m^2)$ is a lognormal multiplicative noise such that $E[\varepsilon_{nl}^{(m)}] = 1$. The noise terms $\varepsilon_{nl}^{(a)}$ and $\varepsilon_{nl}^{(m)}$ are assumed independent with known standard deviations $\sigma_{a,nl}$ and $\sigma_m$, respectively. For low ratios $\tilde{f}_{\ell}(\theta_n)/\sigma_{a,nl}$ low intensities, the additive noise dominates; for high ratios observations are mainly damaged by the multiplicative noise.

**Taking a mixture of noises into account** – The likelihood function associated with the observation model (Eq. 5.1) has no simple closed-form expression due to the mixture of noises. We first consider the uncensored part of the model from Eq. 5.1. In statistics, there are three main approaches to handle such a mixture of noises.

**Approach 1: neglecting one source of noise** – The first approach simply consists in neglecting one source. For instance, similar mixtures of additive and multiplicative noises also occur in medical ultrasound imaging (Krissian et al., 2007) or laser imaging and synthetic aperture radars (Durand et al., 2010). In these fields, when turning to inference, one of the two noises is generally neglected for sake of tractability (Krissian et al., 2007). As already mentioned, the Meudon PDR code spans several decades which causes the nature of the dominant noise to depend on the amplitude of $\tilde{f}_{\ell}(\theta)$. The additive noise can be neglected when $\tilde{f}_{\ell}(\theta_n) \gg \sigma_{a,nl}$, while the multiplicative noise becomes negligible when $\tilde{f}_{\ell}(\theta_n) \ll \sigma_{a,nl}$. As each of the two noise sources dominates in a physical regime, both need to be taken into account in the inversion.

**Approach 2: using a hierarchical model** – A second approach, applied e.g., in Kelly et al. (2012), relies on hierarchical models. An auxiliary variable $U \in R^{N \times L}$ such that $u_{nl} = \varepsilon_{nl}^{(m)} \tilde{f}_{\ell}(\theta_n)$ is included. The observation is then rewritten as

$$\begin{cases} \ln u_{nl} = \ln \tilde{f}_{\ell}(\theta_n) + \ln \varepsilon_{nl}^{(m)}, & \ln \varepsilon_{nl}^{(m)} \sim N(-\sigma_m^2/2, \sigma_m^2), \\ y_{nl} = u_{nl} + \varepsilon_{nl}^{(a)}, & \varepsilon_{nl}^{(a)} \sim N(0, \sigma_{a,nl}^2), \end{cases} \quad (5.2)$$

and $U$ needs to be sampled and inferred along with $\Theta$, yielding the augmented posterior

$$\pi(\Theta, U|Y) \propto \pi(Y|U)\pi(U|\Theta)\pi(\Theta). \quad (5.3)$$

The two likelihood terms $\pi(Y|U)$ and $\pi(U|\Theta)$ have simple closed-form expressions, as they are Gaussian and lognormal, respectively. This approach permits the use of the exact likelihood model. Such a model is generally sampled with a Gibbs algorithm, i.e., that alternatively samples from $\pi(\Theta|U, Y)$ and $\pi(U|\Theta, Y)$. It doubles the number of parameters to sample, which can cause memory issues with large observation maps. This divide-to-conquer approach compensates for this additional cost if sampling from the conditional distributions is simple, i.e., if the coupling between $U$ and $\Theta$ is limited.

Preliminary results on this approach with the considered case, not reported in the manuscript, suggested that sampling is hard in practice although conceptually simple. Unlike in Kelly et al. (2012), $\tilde{f}_{\ell}$ is not gradient Lipschitz continuous and thus neither is $-\ln \pi(U|\Theta)$. This absence of gradient Lipschitz regularity and the changes in dominant noise creates a high coupling between $\Theta$ and $U$, which makes the extended posterior difficult to sample. In the region of the parameter space where the additive noise dominates, exploring the distribution on $\theta_n$ is likely to be slow because the multiplicative noise is negligible. Similarly, exploring the distribution on $u_{nl}$ is likely to be slow when the additive noise is negligible. Avoiding a hierarchical model permits to always consider the dominant noise.
A possibility we did not explore yet is to resort to the ancillarity-sufficiency interweaving strategy (ASIS) sampling approach, dedicated to reducing coupling in hierarchical models. This has already been applied in dust studies, including in Kelly et al. (2012). Another option to reduce the coupling would be to resort to a splitting algorithm – see e.g., Vono et al. (2019).

**Approach 3 (proposed approach): approximating the likelihood function** – The third approach approximates the full likelihood function. For instance, Nicholson and Kaipio (2020) approximates the mixture using a purely additive Gaussian model. The approach we propose builds on this approximation. The additive noise $\epsilon_{n\ell}^{(a)}$ in Eq. 5.1 can be neglected when $\tilde{f}_{\ell}(\theta_n) \to \infty$, while the multiplicative noise $\epsilon_{n\ell}^{(m)}$ becomes negligible as $\tilde{f}_{\ell}(\theta_n) \to 0$. Therefore, for each observation $y_{n\ell}$, the true likelihood is approximated using three different regimes: low, intermediate and high values of $\tilde{f}_{\ell}(\theta_n)$. In the low value regime, the true likelihood function $\pi(y_{n\ell}|\theta_n)$ is approximated by an additive Gaussian approximation $\pi^{(a)}(y_{n\ell}|\theta_n)$ corresponding to

$$y_{n\ell} \approx \tilde{f}_{\ell}(\theta_n) + e_{n\ell}^{(a)}, \quad e_{n\ell}^{(a)} \sim \mathcal{N}(m_{a,n\ell}, s_{a,n\ell}^2), \qquad (5.4)$$

where $m_{a,n\ell}$ and $s_{a,n\ell}^2$ are obtained by matching the two first moments with the model from Eq. 5.1, which yields

$$\begin{cases} m_{a,n\ell} = 0, \\ s_{a,n\ell}^2 = \tilde{f}_{\ell}(\theta_n)^2(e^{\sigma_m^2} - 1) + \sigma_{a,n\ell}^2. \end{cases} \qquad (5.5)$$

Conversely, in the high value regime, a multiplicative lognormal approximation $\pi^{(m)}(y_{n\ell}|\theta_n)$ is used. It reads

$$y_{n\ell} \approx e_{n\ell}^{(m)} \tilde{f}_{\ell}(\theta_n), \quad e_{n\ell}^{(m)} \sim \text{Lognormal}(m_{m,n\ell}, s_{m,n\ell}^2), \qquad (5.6)$$

where moment matching with Eq. 5.1 yields:

$$\left\{ \begin{aligned} m_{m,n\ell} &= -\frac{1}{2} \left\{ \sigma_m^2 + \ln \left[ 1 + \frac{\sigma_{a,n\ell}^2}{\tilde{f}_{\ell}(\theta_n)^2 e^{\sigma_m^2}} \right] \right\}, \\ s_{m,n\ell}^2 &= -2 m_{m,n\ell} \quad \text{so that } \mathbb{E}[e_{n\ell}^{(m)}] = 1. \end{aligned} \right. \qquad (5.7)$$

For the intermediate regime, for each channel $\ell$, we introduce parameters $a_\ell = (a_{\ell,0}, a_{\ell,1}) \in \mathbb{R}^2$. $a_{\ell,0}$ pinpoints the low to intermediate value transition and $a_{\ell,1}$ the intermediate to high value transition. In this intermediate regime, i.e., $a_{\ell,0} \le \tilde{f}_{\ell}(\theta_n) \le a_{\ell,1}$, we propose to use a geometric average of the two likelihood approximations $\pi^{(a)}(y_{n\ell}|\theta_n)$ and $\pi^{(m)}(y_{n\ell}|\theta_n)$ with weights $1 - \lambda$ and $\lambda$, respectively, see the first term of Eq. 5.9 below. The weight function $\lambda$ is defined as a twice differentiable sigmoid with values in $[0, 1]$:

$$\lambda(\theta_n, a_\ell) = \begin{cases} 0 & \text{if } \tilde{f}_{\ell}(\theta_n) \le a_{\ell,0} \\ 1 & \text{if } \tilde{f}_{\ell}(\theta_n) \ge a_{\ell,1}, \\ Q\left(\frac{\ln \tilde{f}_{\ell}(\theta_n) - \ln a_{\ell,0}}{\ln a_{\ell,1} - \ln a_{\ell,0}}\right) & \text{otherwise} \end{cases} \qquad (5.8)$$

where $Q$ is a polynomial such that $Q(0) = 0$, $Q(0) = 1$ and $Q'(0) = Q'(1) = Q''(0) = Q''(1) = 0$ for $\lambda$ to be $\mathscr{C}^2$. One of the simplest such polynomials is $Q(u) = u^3(6u^2 - 15u + 10)$. Figure 5.1 illustrates the $\lambda$ function. The accuracy of this likelihood approximation clearly depends on the choice of the parameter $a_\ell$. Appendix 5.B explains how to optimize $a_\ell$ to maximize the likelihood approximation quality.

**Censorship** – To take censorship into account, let $\mathbf{C} = (c_{n\ell})_{n,\ell} \in \{0, 1\}^{NL}$ be a matrix such that $c_{n\ell} = 1$ for a censored observation, and $c_{n\ell} = 0$ otherwise. Let $F^{(a)}(\cdot|\theta_n)$ and $F^{(m)}(\cdot|\theta_n)$ be the cumulative density function (cdf) of $\pi^{(a)}(\cdot|\theta_n)$ and $\pi^{(m)}(\cdot|\theta_n)$, respectively. The likelihood of
<figure><img src="image_37.png" /><figcaption>Figure 5.1 – Illustration of the $\lambda$ function from Eq. 5.8.</figcaption></figure>

censored data involves $F^{(a)}(\omega_{n\ell}|\theta_n)$ and $F^{(m)}(\omega_{n\ell}|\theta_n)$. The proposed likelihood approximation of the model from Eq. 5.1 finally reads

$$ \tilde{\pi}(y_{n\ell}|\theta_n, \mathbf{a}_{\ell}) \propto \left[ \pi^{(a)}(y_{n\ell}|\theta_n)^{1-\lambda(\theta_n, \mathbf{a}_{\ell})} \pi^{(m)}(y_{n\ell}|\theta_n)^{\lambda(\theta_n, \mathbf{a}_{\ell})} \right]^{1-c_{n\ell}} \\ \times \left[ F^{(a)}(\omega_{n\ell}|\theta_n)^{1-\lambda(\theta_n, \mathbf{a}_{\ell})} F^{(m)}(\omega_{n\ell}|\theta_n)^{\lambda(\theta_n, \mathbf{a}_{\ell})} \right]^{c_{n\ell}}. \quad (5.9) $$

### 5.1.2 Prior distribution

In this thesis we consider maps of physical parameters $\Theta \in \mathbb{R}^{N \times D}$, with $N$ the number of pixels in the map, and $D$ the number of physical parameters per pixel. We combine two penalties to build the prior distribution on $\Theta$. The first provides validity intervals on the physical parameters, as is common in ISM studies (see Chapter 3, Section 3.1.3). The second exploits the map structure in $\Theta$ with a spatial regularization.

**Smooth approximation of uniform distribution** – The first penalty term encodes the validity of $\tilde{f}$ on a compact set $C = [l_1, u_1] \times \dots \times [l_D, u_D]$. Note that the reduced model may be defined out of $C$ but will not be considered as valid since it was not trained on such points. The most natural approach would be to use the indicator function $\iota_{C^N}$ of the set $C^N$, where $\iota_{C^N}(\Theta) = 0$ if $\Theta \in C^N$, and $+\infty$ otherwise. The PMALA kernel to be introduced in Section 5.2.1 requires twice differentiability. The indicator function being non-continuous, it can not be used as is. We use a regularized version. In optimization, one common approach is to replace the indicator function with an exterior penalty function (Nocedal and Wright, 2006). Arguably one of the most common exterior penalty is the quadratic penalty $\theta_{nd} \mapsto [\max(0, \theta_{nd} - u_d, l_d - \theta_{nd})]^2$ (Nocedal and Wright, 2006) which is differentiable and gradient Lipschitz, but not twice differentiable. To obtain a twice differentiable approximation of $\iota_{C^N}$, we use the quartic penalty:

$$ \tilde{\iota}_{C^N} : \Theta \mapsto \sum_{n=1}^{N} \sum_{d=1}^{D} [\max(0, \theta_{nd} - u_d, l_d - \theta_{nd})]^4, \quad (5.10) $$

which is not gradient Lipschitz.

**Spatial regularization** – The second penalty term favors the spatial regularity of estimations. It is based on a local regularizer $h: \mathbb{R}^N \to \mathbb{R}_+$ applied to each map $\Theta_{.d} = (\theta_{nd})_{1 \le n \le N}$, with $d \in [1, D]$. The regularizer can be the Euclidean norm of the usual gradient or Laplacian of the component map, with regularization parameter $\tau_d > 0$. In this thesis, for each of the $D$ maps, we use a $L_2$ norm on the map Laplacian $h(\Theta_{.d}) = \|\Delta\Theta_{.d}\|_2^2$. This regularization function $h$ is twice differentiable. This spatial regularization was already used in the ISM community, e.g., in the ROHSA code (Marchal et al., 2019).
Overall, the resulting prior distribution is given by

$$ \pi(\mathbf{\Theta}) \propto \exp \left( -\xi \tilde{\iota}_{C^N}(\mathbf{\Theta}) - \sum_{d=1}^{D} \tau_d h(\mathbf{\Theta}_{\cdot d}) \right), \quad (5.11) $$

where $\xi > 0$ is a penalty parameter. The higher $\xi$, the better the approximation of the indicator function and the higher the penalty out of $C$.

**Hyperparameters** – This prior relies on two hyperparameters: $\xi > 0$ weighs the smooth indicator function, and $\tau$ contains the $D$ spatial regularization weights of the $D$ maps to be inferred. To learn the spatial regularization weights from the data as in Galliano (2018) would avoid the tedious usual manual setting. In Appendix 5.A, we present a hierarchical model to learn the prior parameters from the data, along with the reasons why we did not use it. Therefore, in the remainder of this thesis, the prior parameters are set manually and iteratively to obtain a trade-off between physically and visually consistent maps. To avoid underestimating uncertainties and favoring physical consistency, the spatial regularization is set to low values, i.e., $\tau_d \in [1, 10]$ for each physical parameter $d$. By default, $\tau_d$ is set to 1.

### 5.1.3 Posterior distribution

The posterior distribution combines $NL$ independent likelihoods (Eq. 5.9) and the prior (Eq. 5.11):

$$ \pi(\mathbf{\Theta}|\mathbf{Y}) \propto \left[ \prod_{n=1}^{N} \prod_{\ell=1}^{L} \tilde{\pi}(y_{nl}|\theta_n, a_\ell) \right] \pi(\mathbf{\Theta}). \quad (5.12) $$

This posterior distribution is hard to manipulate as is, and requires sampling to derive estimators and credibility intervals. However, drawing samples from this posterior is challenging since it is non-log-concave, potentially multimodal. Besides, no gradient Lipschitz continuity is assumed for the log-posterior. In the following and as in Chapter 2, the negative log-posterior pdf $-\ln \pi(\mathbf{\Theta}|\mathbf{Y})$ is denoted $\mathcal{L}(\mathbf{\Theta})$.

## 5.2 Proposed MCMC algorithm

As the forward model spans several decades, the gradient of the negative log-likelihood $\nabla \mathcal{L}$ has a potentially very large Lipschitz constant, if any. Besides, the smooth uniform component of the prior is not gradient Lipschitz. The negative log posterior is therefore not gradient Lipschitz either. In addition, the posterior distribution (Eq. 5.12) is in general non-log-concave, potentially multimodal, which makes the sampling task challenging.

To address these two challenges, a new transition kernel is proposed as a combination of two kernels: PMALA (Xifara et al., 2014) and MTM (Liu et al., 2021). PMALA tackles the regularity issue to efficiently explore the neighborhood of a local mode, whereas MTM permits jumps between modes.

### 5.2.1 PMALA transition kernel

MALA and HMC rely on a step size inversely proportional to the Lipschitz constant of $\nabla \mathcal{L}$, if it exists. Here the forward model $\tilde{f}$ covers several decades so that this Lipschitz constant is potentially very large or even infinite. Therefore, MALA and HMC will typically fail to efficiently explore the posterior (Eq. 5.12). To accelerate the exploration, a preconditioned MALA equipped with RMSProp (Tieleman and Hinton, 2012) is introduced to perform larger steps. To simplify notation, we temporarily use the vector version of $\mathbf{\Theta}$ in lexicographic order so that $\mathbf{\Theta} \in \mathbb{R}^{ND}$.

A transition kernel that handles such situations relies on extensions of HMC and MALA to Riemannian manifolds, introduced in Chapter 2 (Section 2.2.2.5). The Riemannian manifold
MALA version was improved in Xifara et al. (2014), resulting in the so-called preconditioned Metropolis adjusted Langevin algorithm (PMALA) kernel. It permits to exploit local information geometry thanks to a position-dependent preconditioner. We propose to use the RMSProp preconditioner (Tieleman and Hinton, 2012) that was initially defined in the deep learning literature for fast neural networks training. It adaptively estimates a local variance of the gradient $\nabla \mathcal{L}$ by keeping memory of former proposals $\Theta_c^{(t)}$. At each step $t$, it updates a surrogate gradient variance vector $v^{(t)} \in \mathbb{R}^{ND}$ such that for all $i \in [1, ND]$,

$$v_i^{(t)} = av_i^{(t-1)} + (1-a) \left[ \frac{\partial \mathcal{L}}{\partial \theta_i} (\Theta_c^{(t)}) \right]^2 = (1-a) \sum_{j=1}^{t} a^{t-j} \left[ \frac{\partial \mathcal{L}}{\partial \theta_i} (\Theta_c^{(t-j)}) \right]^2, \quad (5.13)$$

where $a \in ]0, 1[$ is an exponential decay rate. Note that the variance vector $v^{(t)}$ relies on candidates $\Theta_c^{(t)}$ instead of iterates $\Theta^{(t)}$: candidates might not be kept in the Markov chain, but they still contain important information about the shape of the distribution. The RMSProp preconditioner is defined as (Tieleman and Hinton, 2012)

$$G^{(t)} = \mathrm{diag} \left( \frac{1}{\epsilon + \sqrt{v^{(t)}}} \right) \in \mathbb{R}^{ND \times ND}, \quad (5.14)$$

with $\epsilon$ a small damping parameter. This preconditioner has already been used in a MCMC context (Li et al., 2016) within an approximate sampler. The goal in Li et al. (2016) was to sample from a distribution defined over the parameters of a neural network trained on a large dataset. Accept or reject steps were omitted as they would have required expensive computations on the full dataset. Additionally, the discretization of the Langevin diffusion process equipped with a position-dependent preconditioner comes with an additional drift term (Xifara et al., 2014) that was neglected in Li et al. (2016). We correct these two approximations to sample exactly from Eq. 5.12.

Following Xifara et al. (2014), the proposal distribution corresponding to PMALA with the RMSProp preconditioner is the Gaussian distribution:

$$q(\Theta_c^{(t)} | \Theta^{(t-1)}) = \mathcal{N}(\Theta_c^{(t)} | \mu^{(t)}, \Lambda^{(t)}) \quad (5.15)$$

with

$$\left\{ \begin{array}{l} \mu^{(t)} = \Theta^{(t-1)} - \eta G^{(t-1)} \nabla \mathcal{L}(\Theta^{(t-1)}) + 2\eta \gamma^{(t-1)}, \\ \Lambda^{(t)} = 2\eta G^{(t-1)}, \end{array} \right. \quad (5.16)$$

where $\eta$ is a step size and $\gamma^{(t-1)}$ is the additional drift term due to the position-dependent preconditioner (Xifara et al., 2014). In full generality, for all $i \in [1, ND]$,

$$\gamma_i^{(t-1)} = \frac{1}{2} \sum_{j=1}^{ND} \frac{\partial G_{ij}^{(t-1)}}{\partial \theta_j^{(t-1)}}. \quad (5.17)$$

However, the RMSProp preconditioner is diagonal so that the sum in Eq. 5.17 reduces to the $j = i$ term only. Note that $\gamma^{(t-1)}$ is defined from a differentiation with respect to iterate $\Theta^{(t-1)}$ while the variance vector $v$ in Eq. 5.13 is defined from candidates. Since all iterates start as candidates, let $j^{(t)}$ be the number of iterations since last accept: $j^{(t)} = \min\{j \ge 0 | \Theta^{(t)} = \Theta_c^{(t-j)}\}$. The correction terms $\gamma_i^{(t-1)}$ are then given by

$$\gamma_i^{(t-1)} = - \frac{(1-a)a^{j^{(t-1)}}}{2\sqrt{v_i^{(t-1)}} \left(\epsilon + \sqrt{v_i^{(t-1)}}\right)^2} \left(\frac{\partial \mathcal{L}}{\partial \theta_i} \cdot \frac{\partial^2 \mathcal{L}}{\partial \theta_i^2}\right) (\Theta^{(t-1)}). \quad (5.18)$$
We now can compute the four components involved in the accept-reject probability $\rho^{(t)}$ (Eq. 2.34) except $q(\Theta^{(t-1)}|\Theta_c^{(t)}) = N(\Theta^{(t-1)}|\mu_c^{(t)}, \Lambda_c^{(t)})$. To compute this last term, one needs to update the variance $\mathbf{v}^{(t)}$ and preconditioner $\mathbf{G}^{(t)}$ and evaluate the candidate additional drift term $\gamma_c^{(t)}$. By definition, $j^{(t)} = 0$ for candidates, so for all $i \in [1, ND]$,

$$ \gamma_{c,i}^{(t)} = - \frac{(1-a)}{2\sqrt{v_i^{(t)}}(\epsilon + \sqrt{v_i^{(t)}})^2} \left( \frac{\partial \mathcal{L}}{\partial \theta_i} \cdot \frac{\partial^2 \mathcal{L}}{\partial \theta_i^2} \right) (\Theta^{(t-1)}). \quad (5.19) $$

The parameters $\mu_c^{(t)}, \Lambda_c^{(t)}$ are thus given by

$$ \left\{ \begin{array}{l} \mu_c^{(t)} = \Theta_c^{(t)} - \eta \mathbf{G}^{(t)} \nabla \mathcal{L}(\Theta_c^{(t)}) + 2\eta \gamma_c^{(t)}, \\ \Lambda_c^{(t)} = 2\eta \mathbf{G}^{(t)}. \end{array} \right. \quad (5.20) $$

Algorithm 5.1 describes the proposed PMALA kernel with RMSProp preconditioner. It relies on three scalar parameters: a damping parameter $\epsilon$, an exponential decay rate $a$ and a step size $\eta$. The first two are generally set to $\epsilon = 10^{-5}$ and $a = 0.99$ (Li et al., 2016). The step size is chosen empirically. MALA achieves optimal convergence rates with an acceptance rate equal to 0.574 when the components of $\Theta$ are independent (Robert and Casella, 2004). Despite the interdependencies in the posterior distribution, we also set $\eta$ to obtain an average acceptance rate close to 0.574, which yields good results in practice.

<caption>**Algorithm 5.1:** PMALA kernel $K_1$ at step $t$</caption>

Input: $\Theta^{(t-1)}, \mathbf{v}^{(t-1)}, j^{(t-1)}$
Output: $\Theta^{(t)}, \mathbf{v}^{(t)}, j^{(t)}$
// Propose candidate
1 $\mathbf{G}^{(t-1)}$ and $\gamma^{(t-1)}$ // using Eq. 5.14 and Eq. 5.18
2 $\mu_c^{(t)}$ and $\Lambda_c^{(t)}$ // using Eq. 5.16
3 $\Theta_c^{(t)} \sim N(\mu_c^{(t)}, \Lambda_c^{(t)})$
// Accept or reject
4 $\mathbf{v}^{(t)}, \mathbf{G}^{(t)}$ and $\gamma_c^{(t)}$ // using Eq. 5.13, Eq. 5.14, and Eq. 5.19
5 $\mu_c^{(t)}, \Lambda_c^{(t)}$ and $\rho^{(t)}$ // using Eq. 5.20, and Eq. 2.34
6 Draw $\zeta \sim \text{Unif}(0, 1)$
7 if $\zeta \le \rho^{(t)}$ then $\Theta_c^{(t)} = \Theta_c^{(t-1)}, j^{(t)} = 0$
8 else $\Theta_c^{(t)} = \Theta_c^{(t-1)}, j^{(t)} = j^{(t-1)} + 1$

Accounting for the correction term $\gamma$ is necessary in stochastic gradient MCMC (SG-MCMC) as it does not apply accept-reject steps. Alternatively, disregarding it, as in Li et al. (2016), leads to a controlled error on estimations. In our case, as we do apply a MH accept-reject step, computing the correction term is not mandatory to obtain an MCMC algorithm that converges to the correct stationary distribution $\pi(\Theta|\mathbf{Y})$. However, we keep this term to obtain a correct discretization of the associated Langevin diffusion process. As stated in Chapter 2 (Section 2.2.2.4), better discretization of the Langevin process often leads to better candidates.

## 5.2.2 MTM transition kernel

The non-log-concavity and potential multimodality of the posterior (Eq. 5.12) is the second major difficulty to be addressed. As explained in Chapter 2, samplers such as MH, MALA, HMC or even PMALA fail to explore the full distribution when modes are far away: they get stuck in one. In Section 2.A, we reviewed existing algorithms dedicated to multimodal distributions. None
of these methods can efficiently address the very high dimensionality $O(10^4)$ to be encountered in some of the applications considered in this thesis.

We define a kernel that can escape a local mode and explore other ones without any knowledge about the number, positions or variances of the modes. To do so, we propose to exploit two particularities of the considered posterior, namely its map structure and the speed of evaluation of its forward model - and thus of its full pdf. Instead of sampling the whole vector $\Theta \in \mathbb{R}^{ND}$ at once, it uses a Metropolis-within-Gibbs sampler to decompose it into $N$ individual $\theta_n$ (see Chapter 2, Section 2.2.2.6). For each conditional distribution, it harnesses an independent multiple-try Metropolis (I-MTM) approach (Liu et al., 2021; Martino, 2018) (see Chapter 2, Section 2.2.2.7). This method generates $K \ge 1$ candidates $(\theta_n^{(k)})_{k=1}^K$ independently of $\theta_n^{(t-1)}$. This divide-to-conquer approach permits to consider $N$ conditional distributions $\pi(\theta_n|\mathbf{y}_n, \Theta_n^{(t-1+(n-1)/N)})$ of small dimension, where $\Theta_n^{(t-1+(n-1)/N)} = (\theta_1^{(t)}, \dots, \theta_{n-1}^{(t)}, \theta_{n+1}^{(t-1)}, \dots, \theta_N^{(t-1)})$. Candidates are sampled from a proposal distribution $q(\theta_n|\Theta_n^{(t-1)})$ that should be permissive enough to generate candidates in all modes of $\pi(\theta_n|\mathbf{y}_n, \Theta_n^{(t-1)})$. Then, using the importance weight function $w$ from Eq. 2.42 adapted to this case,

$$w(\theta_n^{(k)}) = \frac{\pi(\theta_n^{(k)}|\mathbf{y}_n, \Theta_n^{(t-1)})}{q(\theta_n^{(k)}|\Theta_n^{(t-1)})}, \qquad (5.21)$$

one candidate is selected using a categorical distribution with selection probability $w_k$ for candidate $k$

$$w_k = \frac{w(\theta_n^{(k)})}{\sum_{j=1}^{K} w(\theta_n^{(j)})}. \qquad (5.22)$$

The MH step is then performed with the selected candidate $i$ and the following generalized acceptance probability, adapted from Eq. 2.44,

$$\hat{\rho}_n^{(t)} = \min \left( 1, \frac{w(\theta_n^{(i)}) + \sum_{j=1, j \neq i}^{K} w(\theta_n^{(j)})}{w(\theta_n^{(t-1)}) + \sum_{j=1, j \neq i}^{K} w(\theta_n^{(j)})} \right). \qquad (5.23)$$

Algorithm 5.2 summarizes the proposed Gibbs and MTM kernel. Note that due to the Gibbs approach, it updates one component at a time and returns the result of all updates. A succession of intermediate updates $(\Theta_n^{(t-1+n/N)})_{n=1}^N$ is therefore introduced (see Algorithm 2.3). This transition kernel relies on the choice of the proposal distribution $q$ and on the number of candidates $K$ generated at each step. This parameter is chosen as a trade-off between computational intensity and average acceptance probability: the higher $K$, the higher the acceptance probability, the mixing capability but also the computational cost.

In image inverse problems, many common spatial priors are based on a local operator such as the image gradient or Laplacian. In such cases, many components $\theta_n$ are conditionally independent. They can be sampled in parallel using a Chromatic Gibbs sampler (Gonzalez et al., 2011), which can significantly speed up computations. Note that this Gibbs sampling and chromatic Gibbs sampling can also be adopted in the PMALA kernel to farther accelerate local explorations.

### 5.2.3 Proposed sampler and implementation details

To combine a good local exploration of modes as well as jumps between modes, the proposed kernel mixes the PMALA and MTM transition kernels above. At every step $t$, the MTM kernel is selected with probability $p_{MTM}$, and the PMALA kernel with probability $1-p_{MTM}$. Since the MTM kernel divides the parameter space in $N$ $D$-dimensional subspaces, the PMALA global integer $j^{(t)} \ge 0$ is replaced by a vector $j^{(t)} \in \mathbb{N}^N$, where $j_n^{(t)}$ counts the number of steps since last acceptance for component $\theta_n$. When a component $\theta_n$ is accepted by the MTM kernel, the
<table><caption>**Algorithm 5.2: MTM kernel $K_2$ at step t**</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>for n = 1 to N do</td><td></td></tr><tr><td>2</td><td>// Propose candidates, select one</td><td></td></tr><tr><td>3</td><td>θ<sub>n</sub><sup>(k)</sup> ~ q(θ<sub>n</sub> | Θ<sub>n</sub><sup>(t-1+<sup>n-1</sup>/<sub>N</sub>)</sup>) for k = 1 to K</td><td></td></tr><tr><td>4</td><td>w(θ<sub>n</sub><sup>(k)</sup>) for k = 1 to K</td><td>// using Eq. 5.21</td></tr><tr><td>5</td><td>w<sub>k</sub> for k = 1 to K</td><td>// using Eq. 5.22</td></tr><tr><td>6</td><td>i ~ Cat(w<sub>1</sub>, ..., w<sub>K</sub>)</td><td></td></tr><tr><td>7</td><td>// Accept or reject</td><td></td></tr><tr><td>8</td><td>ρ̃<sub>n</sub><sup>(t)</sup></td><td>// using Eq. 5.23</td></tr><tr><td>9</td><td>Draw ζ<sub>n</sub> ~ Unif(0, 1)</td><td></td></tr><tr><td>10</td><td>if ζ<sub>n</sub> ≤ ρ̃<sub>n</sub><sup>(t)</sup> then</td><td></td></tr><tr><td></td><td>θ<sub>n</sub><sup>(t-1+<sup>n</sup>/<sub>N</sub>)</sup> = θ<sub>n</sub><sup>(i)</sup>, Θ<sub>n</sub><sup>(t-1+<sup>n</sup>/<sub>N</sub>)</sup> = Θ<sub>n</sub><sup>(t-1+<sup>n-1</sup>/<sub>N</sub>)</sup></td><td></td></tr><tr><td></td><td>else Θ<sub>n</sub><sup>(t-1+<sup>n</sup>/<sub>N</sub>)</sup> = Θ<sub>n</sub><sup>(t-1+<sup>n-1</sup>/<sub>N</sub>)</sup></td><td></td></tr></tbody></table>

counter $j_n$ is reset to 0 and the variance component $v_n \in \mathbb{R}^D$ is updated as in Eq. 5.13 with $\frac{\partial \mathcal{L}}{\partial \theta_n}(\Theta^{(t)})$.

Algorithm 5.3 reports the complete proposed sampler. Similarly to RDMC and WHMC, the proposed sampler mixes a kernel dedicated to local exploration – PMALA – and another to jump between modes – MTM. The decomposition of the parameter space into $N$ $D$-dimensional subspaces makes the sampler much simpler than previous approaches. It will perform well in structured problems that allow such decomposition, e.g., images and graphs, and poorly in high-dimensional problems that do not, e.g., Gaussian Mixtures over the full space. In particular, resorting to an Gibbs and MTM strategy leads to fast convergence to the high probability regions. Therefore, the burn-in phase is shorter than with usual MCMC algorithms. Using the Gibbs approach in the PMALA kernel also accelerates convergence and local exploration.

An optimization algorithm can be obtained from this sampler. It combines preconditioned gradient descent and simulated annealing-like component-wise updates that can escape from local minima.

Regarding theoretical properties, the PMALA kernel satisfies the detailed balance property – from Robert and Casella (2004, theorem 7.2) – and produces ergodic Markov chains – from Robert and Casella (2004, corollary 7.5). The proposed MTM kernel is a Metropolis-within-Gibbs algorithm with propositions independent to the current location and with multiple candidates $K$. In the particular case where $K = 1$, it satisfies the detailed balance property and produces uniformly ergodic Markov Chains – from Jones et al. (2014, theorem 7). Using $K > 1$ candidates in a MTM framework maintains detailed balance and ergodicity (Martino, 2018). As a mixture of kernels having the same stationary distribution, the proposed kernel also admits the posterior as a stationary distribution – from Robert and Casella (2004, chapter 10). As the MTM kernel produces uniformly ergodic Markov chains, so does the proposed mixture kernel – from Robert and Casella (2004, proposition 10.20). These results of convergence towards the posterior are mostly asymptotic and also hold for simpler algorithms such as RWMH (Roberts and Tweedie, 1996). A comparative theoretical study of non-asymptotic properties that could demonstrate a faster convergence of the proposed sampler is beyond the scope of this manuscript. However, Section 5.2.4 presents empirical results showing that the proposed sampler yields state-of-the-art performance on multimodal distributions in low-dimensional settings. As reported in Section 5.4, the proposed sampler also yields state-of-the-art performance on higher dimensional applications
<figure><img src="image_38.png" /><figcaption>**Algorithm 5.3:** Proposed sampler: PMALA and MTM</figcaption></figure>

with relevant low-dimensional conditional distributions.

### 5.2.4 Illustration: 2D Gaussian mixture model

The proposed sampler is applied to a two-dimensional Gaussian mixture model (GMM) restricted to the square $C = [-15, 15]^2$. This simple multimodal distribution, shown on Figure 5.2 (top left), is set to contain 15 modes $(\mu_i, \Sigma_i)$. It will demonstrate the ability of the proposed sampler to jump between modes. For simplicity, all the modes have an equal weight in the mixture

$$ \pi(\theta) \propto \left[ \sum_{i=1}^{15} \mathcal{N}(\theta | \mu_i, \Sigma_i) \right] \exp(-\xi \tilde{l}_C(\theta)), \quad (5.24) $$

with $\xi = 10^4$. No natural structure decomposition exists for a GMM since each observation consists of $N = 1$ point only in dimension $D = 2$. A Markov chain composed of $T_{MC} = 10000$ samples is considered, including $T_{BI} = 100$ burn-in samples. To illustrate the role of each of the two kernels in the proposed sampler, two different values are considered for the probability of selecting the MTM kernel: $p_{MTM} = 0.1$ or $p_{MTM} = 0.9$. The number of candidates of the MTM kernel is set to $K = 50$, and the proposal distribution $q$ is the smooth uniform prior on $C$ – see Appendix 5.C. The MTM candidates weights $w(\theta_n^{(k)})$ in Eq. 5.21 are then equal to the likelihood term, i.e., the sum of Gaussian pdfs. The default values $a = 0.99$ and $\epsilon = 10^{-5}$ are considered for the exponential decay and damping factor of the PMALA kernel (Li et al., 2016), and its step size is set to $\eta = 0.25$. The proposed approach is compared to the state-of-the-art wormhole Hamiltonian Monte Carlo (WHMC) (Lan et al., 2014), using the same number of samples. Note that WHMC needs the prior knowledge of mode positions $(\mu_i)_{1\le i\le 15}$, while the proposed kernel does not.

The three sampling algorithms widespread in astrophysics introduced in Chapter 2 (Appendix 2.A) are also applied. The affine-invariant sampler (Goodman and Weare, 2010) is run with 4 parallel chains, using the EMCEE implementation (Foreman-Mackey et al., 2013). The sequential MC (SMC) algorithm (Del Moral et al., 2006) is run with the state-of-the-art PYMC code<fnref n="1" /> (Abril-Pla et al., 2023). The MULTINEST algorithm (Feroz and Hobson, 2008) is run with

<footnote n="1">1. https://www.pymc.io/welcome.html</footnote>
its Python implementation PYMULTINEST (Buchner, 2016a). All three algorithms are run with their respective default parameters. Note that these three samplers are not MCMC algorithms.

Figure 5.2 shows the 2D histograms obtained with the six samplers. The affine-invariant sampler is the only sampler that fails to explore all the modes. The five others visit all the modes, approximately give them equal weights in the histogram, and their local dispersion obeys the covariance structures equally well. Table 5.1 compares their bias<fnref n="2" /> and effective sample size (ESS) for the MCMC algorithms – see Chapter 2 (Section 2.2.2) for the definition of the ESS. When $p_{MTM} = 0.9$, the proposed sampler achieves the second lowest bias and the highest ESS. In particular, it achieves better performances than WHMC, despite the absence of information about the position of the modes $\mu_i$. The high ESS values result from the 85% acceptance rate of the MTM kernel for $K = 50$. However, the MTM kernel with a fixed number of candidates $K$ would not scale up to much higher dimensions. The probability to jump between modes is proportional to the volume of the high probability regions compared to the volume of $C$, and thus decreases exponentially with the dimension of the problem. The proposed sampler would therefore fail to reach isolated modes in a high-dimensional GMM, whereas WHMC would succeed to do so by exploiting its additional information about the modes. However, the proposed approach focuses on scenarios where the parameter space can be partitioned into a collection of $N$ subspaces of limited dimension $D$, typically $D \lesssim 10$. The MTM kernel thus remains out of reach from the curse of dimension thanks to the structure of the problem. As in this simple GMM example, the proposed sampler can then outperform WHMC, even without any prior information on the modes of a multimodal distribution.

<table><caption>Table 5.1 – Samplers comparison on 2D Gaussian mixture model.</caption><thead><tr><th rowspan="2">Algorithm</th><th>Bias</th><th colspan="2">ESS</th></tr><tr><th>||<span style="text-decoration: overline;">θ̂</span><sub>MMSE</sub> - E[θ]||</th><th>θ<sub>1</sub></th><th>θ<sub>2</sub></th></tr></thead><tbody><tr><td>Affine-invariant sampler (Goodman and Weare, 2010)</td><td>4.32 · 10<sup>0</sup></td><td>—</td><td>—</td></tr><tr><td>sequential MC (Del Moral et al., 2006)</td><td>1.09 · 10<sup>-1</sup></td><td>—</td><td>—</td></tr><tr><td>MULTINEST (Feroz and Hobson, 2008)</td><td><b>3.22 · 10<sup>-2</sup></b></td><td>—</td><td>—</td></tr><tr><td>WHMC (Lan et al., 2014)</td><td>1.28 · 10<sup>-1</sup></td><td>2 753</td><td>2 993</td></tr><tr><td>Proposed, p<sub>MTM</sub> = 0.1</td><td>7.02 · 10<sup>-1</sup></td><td>395</td><td>444</td></tr><tr><td>Proposed, p<sub>MTM</sub> = 0.9</td><td>4.61 · 10<sup>-2</sup></td><td><b>6 157</b></td><td><b>5 780</b></td></tr></tbody></table>

<footnote n="2">2. To compute the bias, we consider that $E[θ] \approx \frac{1}{15} \sum_{i=1}^{15} \mu_i$, which is inexact due to the smooth indicator term. However, the probability mass of the Gaussian mixture out of the validity intervals is numerically negligible. Therefore, the error due to this approximation is negligible as well.</footnote>
<figure><img src="image_39.png" /><figcaption>Figure 5.2 – Comparison of samplers on a 2D Gaussian mixture model. The red ellipses show the probability level at $2\sigma$. All histograms are in logarithmic norm.</figcaption></figure>
## 5.3 Model checking using a predictive posterior p-value

This section presents a method that checks whether the observation model $\mathcal{M}$ can reproduce the observation $\mathbf{Y}$. As in Chapter 2 (Section 2.1), we consider that the model $\mathcal{M}$ gathers the forward model and the noise model. In astrophysics applications, the forward model is the emulator of the Meudon PDR code defined in Chapter 4. The noise model is described in Chapter 3 (Section 3.4).

To perform this model checking, we resort to the predictive posterior p-value introduced in Chapter 2 (Section 2.3.2). As we showed in Chapter 3 (Section 3.3.2), this approach was already used in astrophysics. We improve it in two ways. First, we set the discrepancy measure to the negative log-likelihood for generality. As the resulting predictive posterior p-value is intractable, we approximate it with a Monte Carlo (MC) estimator. Second, we control the associated approximation error to avoid making wrong decisions because of insufficient number of iterates in our MC estimator.

### 5.3.1 Definition of selected p-value

For inverse problems relying on a non-linear forward model or on a posterior distribution instead of a point estimate $\hat{\Theta}$, there exists no general test statistic $T$ that leads to an exact evaluation of the p-value. The choice of the test statistic should then enforce consistent behavior between the observation model and the p-value. We propose to use the negative log likelihood function as a test statistic,

$$T(\tilde{\mathbf{Y}}, \mathbf{\Theta}) = -\ln \pi(\tilde{\mathbf{Y}}|\mathbf{\Theta}, \mathcal{M}), \qquad (5.25)$$

where $\tilde{\mathbf{Y}}$ are reproductions of observations following the observation model. In addition, as the observation model handles pixels independently, we propose to evaluate one p-value per pixel. Overall, a map $\mathbf{p} = (p_n)_{n=1}^N$ of $N$ p-values is evaluated on the marginal posterior predictive distributions $\tilde{\mathbf{y}}_n|\mathbf{Y}$. The p-value $p_n$ corresponds to the measure of the set

$$I_n = \{(\tilde{\mathbf{y}}_n, \mathbf{\Theta}) | \pi(\tilde{\mathbf{y}}_n|\mathbf{\Theta}, \mathcal{M}) \le \pi(\mathbf{y}_n|\mathbf{\Theta}, \mathcal{M})\}, \qquad (5.26)$$

i.e.,

$$p_n = \int \mathbb{1}_{I_n}(\tilde{\mathbf{y}}_n, \mathbf{\Theta}) \pi(\tilde{\mathbf{y}}_n|\mathbf{\Theta}, \mathcal{M}) \pi(\mathbf{\Theta}|\mathbf{Y}, \mathcal{M}) d\mathbf{\Theta} d\tilde{\mathbf{y}}_n. \qquad (5.27)$$

Small p-values $p_n$ enable detecting anomalies in observed maps, e.g., regions that are not well modeled by a PDR. This definition is general, as it does not rely on restraining assumptions on the observation model. The only two necessary assumptions are 1) that the observation model can be used to generate observations from a given $\mathbf{\Theta}$ and 2) that the likelihood function can be evaluated.

### 5.3.2 p-value approximation and associated uncertainties

As detailed in Chapter 2 (Section 2.3.2), for general models the p-value is evaluated with a Monte Carlo (MC) estimator $\hat{p}$ (Eq. 2.57). The map of p-values $\hat{p}^{(t)} = (\hat{p}_n^{(t)})_{n=1}^N$ can be evaluated after the burn-in phase of size $T_{BI} \ge 0$:

$$\hat{p}_n^{(t)} = \frac{1}{t - T_{BI}} \sum_{\tau=T_{BI}+1}^{t} \mathbb{1}_{I_n}(\tilde{\mathbf{y}}_n^{(\tau)}, \mathbf{\Theta}^{(\tau)}). \qquad (5.28)$$

Algorithm 5.4 details how this estimator can be effortlessly included in an MCMC algorithm. In case of a point estimator $\hat{\Theta}$, the Markov chain is replaced by $\mathbf{\Theta}^{(0)} = \dots = \mathbf{\Theta}^{(T_{MC})} = \hat{\Theta}$, and $T_{BI} = 0$.
<caption>**Algorithm 5.4:** MCMC sampling with p-value computation</caption>

**Input:** Starting point $\Theta^{(T_{BI})}$, sampling kernel $\mathcal{K}$, numbers of iterations $T_{BI}$ et $T_{MC}$
**Initialization:** map of p-values $\hat{\mathbf{p}}^{(T_{BI})} = \mathbf{0}_N$, $\mathbf{u}^{(T_{BI})} = \mathbf{0}_N$
**for** $t = T_{BI} + 1, \dots, T_{MC}$
// After burn-in phase
**do**
4 &nbsp;&nbsp; $\Theta^{(t)}$, $(\pi(\mathbf{y}_n|\Theta^{(t)}, \mathcal{M})_{n=1}^N = \mathcal{K}(\Theta^{(t-1)})$
// Sample from $\mathcal{K}$
5 &nbsp;&nbsp; $\tilde{\mathbf{Y}}^{(t)} \sim \pi(\cdot|\Theta^{(t)}, \mathcal{M})$
// generate observation reproduction
6 &nbsp;&nbsp; $u_n^{(t)} = u_n^{(t-1)} + \mathbf{1}_{I_n}(\tilde{\mathbf{y}}_n^{(t)}, \Theta^{(t)})$ for $n \in [1, N]$ // update p-values
7 &nbsp;&nbsp; $\hat{\mathbf{p}}^{(T_{MC})} = \frac{1}{T_{MC}-T_{BI}}\mathbf{u}^{(T_{MC})}$
**Output:** Markov chain $\{\Theta^{(t)}\}_{t=T_{BI}+1}^{T_{MC}}$, map of p-values $\hat{\mathbf{p}}^{(T_{MC})}$

Accounting for the uncertainty inherent to a Monte Carlo estimation of the p-values makes hypothesis testing more robust. The probability of null hypothesis rejection with confidence level $\alpha$ then corresponds to the cdf of $p_n$ evaluated on $\alpha$, $\mathbb{P}[p_n \le \alpha]$. Using a threshold $\delta \in [0, 0.5]$, $\mathbb{P}[p_n \le \alpha] \in [\delta, 1 - \delta]$ implies that more samples $(\Theta^{(t)}, \tilde{\mathbf{y}}_n^{(t)})$ are necessary to make a decision with the desired confidence level. The model $\mathcal{M}$ is rejected for pixel $n$ when $\mathbb{P}[p_n^{(t)} \le \alpha] > 1 - \delta$, and is not rejected when $\mathbb{P}[p_n^{(t)} \le \alpha] < \delta$. In other words,

$$ \begin{cases} \mathbb{P}[p_n^{(t)} \le \alpha] > 1 - \delta & \Rightarrow \text{model } \mathcal{M} \text{ is rejected for pixel } n \text{ with confidence level } \alpha. \\ \mathbb{P}[p_n^{(t)} \le \alpha] < \delta & \Rightarrow \text{model } \mathcal{M} \text{ is not rejected for pixel } n \text{ with confidence level } \alpha. \quad (5.29) \\ \mathbb{P}[p_n^{(t)} \le \alpha] \in [\delta, 1 - \delta] & \Rightarrow \text{the sample size is insufficient to make a decision.} \end{cases} $$

In general, evaluating the cdf of the p-value $p_n$ would require an additional MCMC layer. We avoid this additional step and maintain a reasonable cost for this test. First, note from Eq. 5.27 that $1_{I_n}(\tilde{\mathbf{y}}_n, \Theta)$ is a binary random variable that equals 1 with probability $p_n$ and 0 with probability $1 - p_n$. In other words, it follows a Bernoulli distribution of parameter $p_n$. A natural prior distribution for $p_n$ is a uniform distribution on $[0, 1]$, as $p_n$ is a probability. This is a special case of a Beta distribution, Beta(1, 1), which is the conjugate prior of the Bernoulli distribution (Gelman et al., 2015, chapter 2, section 2.4). The uncertainty on $p_n$ thus can be reasonably modeled by an a posteriori Beta distribution with closed-form parameters

$$ p_n^{(t)} \sim \text{Beta} \left( 1 + N_n^{(t)} \hat{p}_n^{(t)}, 1 + N_n^{(t)} (1 - \hat{p}_n^{(t)}) \right), \quad (5.30) $$

for any $t > T_{BI}$, with $N_n^{(t)}$ the number of independent samples. When the posterior on $\Theta$ is reduced to a Dirac on an estimator $\hat{\Theta}$, then all samples are independent and $N_n^{(t)} = t - T_{BI}$. When a Markov chain is generated to sample from the posterior distribution $\pi(\Theta|\mathbf{Y}, \mathcal{M})$, the samples $\Theta^{(t)}$ are correlated. To compensate this correlation and avoid underestimating uncertainties, we set $N_n^{(t)}$ to the effective sample size (ESS) on $\Theta$, i.e., $N_n^{(t)} = \text{ESS}_n^{(t)}$ This approach permits accounting for errors associated with the numerical evaluation of the p-values, which makes the test more robust.

Figure 5.3 shows three cases using this definition of the p-value. In this case, In the first case, the Beta distribution has most of its weight on p-values higher than $\alpha$. The cdf at $\alpha$ is therefore very small. An associated model $\mathcal{M}$ would not be rejected. In the second case, the Beta distribution has most of its weight on p-values lower than $\alpha$. The cdf at $\alpha$ is therefore very high, and higher than the threshold $\delta$. An associated model $\mathcal{M}$ would therefore be rejected. Finally, in
the third case, the Beta distribution has significant weight both below and above the confidence
level $\alpha$. As a consequence, the cdf evaluated at $\alpha$ falls in the $[\delta, 1 - \delta]$ interval. In this case, more
samples would be necessary to make a decision.

<figure><img src="image_40.png" /><figcaption>Figure 5.3 – Illustration of the three-case test for model assessment. In this case, $\alpha = 0.1$ and $\delta = 0.1$.</figcaption></figure>

5.3.3 Illustration: 2D Gaussian distribution

As already shown in Chapter 2 (Section 2.3.2), the Bayesian *p*-value can be efficiently es-
timated in the case of a point estimate $\hat{\Theta}$ and of an additive Gaussian noise model. In this
particular case, the test statistic *T* follows a $\chi^2_{NL}$ distribution, with *N* the number of pixels and
*L* the number of observables per pixel. This subsection illustrates this equivalence and shows
that the proposed test behaves as expected. In particular, we illustrate the convergence of the
distribution on *p*(<sup>*t*</sup>) to the theoretical value *p* at speed $1/\sqrt{t}$. We consider an inverse problem
such that *N* = 1 et *L* = 2, i.e., with an observation *y* ∈ R<sup>2</sup>. The forward model is set to the
identity function *f* = *id*<sub>2</sub>. The noise is assumed Gaussian, additive and with unit variance. We
consider three estimators $\hat{\theta}^{(i)} \in R^2$. The confidence levels are set to $\alpha = 0.05$ and $\delta = 0.1$.

Figure 5.4 shows the convergence of the three *p*-values estimated with Eq. 5.28 to the theoret-
ical values. With $T_{\text{MC}} = 10^3$ independent samples drawn from the posterior predictive distribution
$\pi(\tilde{y}|\hat{\theta}^{(i)}, \mathcal{M})$, the point estimate $\hat{p}$ from Eq. 5.28 leads to a rejection of the null hypothesis from
$\hat{\theta}^{(1)}$ and from $\hat{\theta}^{(2)}$, and to a non-rejection from $\hat{\theta}^{(3)}$. These decisions are in agreement with those
corresponding to the theoretical *p*-values for $\hat{\theta}^{(1)}$ and $\hat{\theta}^{(3)}$. However, for $\hat{\theta}^{(2)}$, the rejection is
in disagreement with the theoretical *p*-values, which is higher than the confidence level $\alpha$. This
rejection is caused by the uncertainties inherent to using a Monte Carlo estimator of the *p*-value
with insufficient number of samples. Using the Beta distribution (Eq. 5.30) to quantify uncer-
tainty, one obtains $P[p_n^{(t)} \le \alpha] = 0.83$. Since $P[p_n^{(t)} \le \alpha] \in [\delta, 1 - \delta]$, the proposed three-case
test (Eq. 5.29) identifies this issue and recommends drawing more samples. With $T_{\text{MC}} = 10^4$,
$P[p_n^{(t)} \le \alpha] = 0.096 < \delta$. The proposed three-case test then considers that enough samples have
been drawn to make a confident decision, and the model associated with $\hat{\theta}^{(2)}$ is not rejected.
This decision is consistent with the corresponding theoretical *p*-value. The proposed three-case
test (Eq. 5.29) is therefore robust to using a Monte Carlo estimator of the *p*-values.
<figure><img src="image_41.png" /><figcaption>Figure 5.4 – Application of Bayesian *p*-value with three-case test on a simple case, with Gaussian noise and point estimate $\hat{\theta}$, and comparison with corresponding theoretical *p*-values. (Top left) pdf of residuals $\tilde{y} - \hat{\theta}$ (in gray levels) and $T_{MC} = 1000$ independent samples $\tilde{y} \sim \pi(\tilde{y}|\hat{\theta}^{(i)})$. The contours correspond to constant values of $||\tilde{y}-\hat{\theta}||_2^2$. (Top right) Distribution $\chi_2^2$ of $||\tilde{y}-\hat{\theta}||_2^2$. The vertical lines correspond to the contour lines in Figure 5.4a. (Bottom) Convergence of *p*-value estimator (Eq. 5.28) and of the 90% credibility intervals, obtained from Eq. 5.30.</figcaption></figure>
## 5.4 Applications

In this section, the proposed sampler is applied to two inverse problems: a sensor localization problem and a synthetic example of the general astrophysics problem addressed in this thesis.

### 5.4.1 Sensor localization

The sensor localization problem introduced in Ihler et al. (2005) is a common test case in multimodal sampling, e.g., in Ahn et al. (2013); Lan et al. (2014); Pompe et al. (2020). Three sensors have known locations and will serve as a reference to avoid ambiguities with respect to translation, rotation and negation. The goal is to estimate the unknown positions $\Theta \in \mathbb{R}^{ND}$ of $N = 8$ sensors in dimension $D = 2$. The observation matrix $\mathbf{Y} \in \mathbb{R}^{NL}$ collects noisy and partially censored pairwise distances, where $L = N + 3$ is the total number of sensors. The distance to sensor $\ell$ feeds channel $\ell$, so that the forward model is $f_{\ell}(\theta_n) = \|\theta_n - \theta_{\ell}\|$. Note that only $N + 2$ distances will really be used since $f_{\ell}(\theta_{\ell}) = 0$, and that we set $y_{n\ell} = 0$ by convention. The probability of communication from sensor $\ell \in [1, L]$ to sensor $n \in [1, N]$ is set to $\exp\{-\frac{f_{\ell}(\theta_n)^2}{2R^2}\}$ with $R = 0.3$. In absence of communication, the observation is censored, which is encoded by the binary latent variable $c_{n\ell} = 1$. Otherwise, $c_{n\ell} = 0$ when the observation occurs and is corrupted by a white Gaussian noise

$$y_{n\ell} = f_{\ell}(\theta_n) + \epsilon_{n\ell}, \quad \text{with } \epsilon_{n\ell} \sim \mathcal{N}(0, \sigma^2), \qquad (5.31)$$

with $\sigma = 0.02$, leading to

$$-\ln \pi(\mathbf{Y}|\mathbf{\Theta}) = \sum_{n=1}^{N} \sum_{\ell=1}^{L} (1 - c_{n\ell}) \left[ \frac{(f_{\ell}(\theta_n) - y_{n\ell})^2}{2\sigma^2} + \frac{f_{\ell}(\theta_n)^2}{2R^2} \right] + c_{n\ell} \ln \left[ 1 - \exp\left(-\frac{f_{\ell}(\theta_n)^2}{2R^2}\right) \right], \qquad (5.32)$$

The smoothed uniform prior on the square $C = [-0.35, 1.2]^2$ is used as a prior on the location of each sensor. The corresponding penalty parameter $\xi$ introduced in Eq. 5.11 is set to $10^4$. This prior is non-informative enough to match the results shown in Ahn et al. (2013); Lan et al. (2014); Pompe et al. (2020). The proposed sampler is compared to both regeneration darting MC (RDMC) (Ahn et al., 2013) and WHMC. A Markov chain of size 30 000 is generated by each algorithm, including 5 000 burn-in samples. The parameters of the PMALA kernel are set to $a = 0.99$, $\epsilon = 10^{-5}$ and $\eta = 1.5 \times 10^{-3}$. The MTM kernel is selected with $p_{\text{MTM}} = 0.1$ or $p_{\text{MTM}} = 0.9$. Its proposal distribution $q$ is set to the smooth uniform prior on $C$. For each sensor, the high probability regions are small compared to $C$. To obtain high acceptance rates for the MTM kernel, the number of candidates is set to $K = 1000$. Better proposal distributions can be obtained for this specific problem, which is beyond the scope of this experiment.

Figure 5.5 shows the marginal distributions of each sensor position. The four samplers identified the same modes. Table 5.2 compares the samplers in terms of ESS. With $p_{\text{MTM}} = 0.9$, the proposed sampler yields better mixing capability than WHMC and RDMC. This is due to the partition of the $ND = 16$-dimensional problem into $N = 8$ simpler $D = 2$-dimensional problems. This divide-to-conquer strategy exploits the problem structure to fight the curse of dimension.

**Bayesian model assessment** – For sake of illustration, one observed distance is divided by a factor 10 to make the observation set inconsistent. The resulting posterior is sampled using the same parameters as above and $p_{\text{MTM}} = 0.9$. Figure 5.6 indicates which observed distance was modified, and compares the marginal posterior distributions with and without this alteration. It appears that few sensors are affected, as most of the marginals are very similar in the two cases. If the ground truth was not known, one could not detect the issue from the marginal posterior distributions alone.
<figure><img src="image_42.png" /><figcaption>Figure 5.5 – Marginal distributions of the sensors positions for the considered samplers. The graph shows the true position of all sensors. The sensors with a known position are in red and those whose position is inferred are in blue. The edges of the graph indicate which pairs of sensors are observed.</figcaption></figure>

<table><caption>Table 5.2 – Effective Sample Size (ESS) on the sensor localization problem.</caption><thead><tr><th rowspan="2">MCMC sampler</th><th colspan="3">ESS</th></tr><tr><th>min</th><th>mean</th><th>max</th></tr></thead><tbody><tr><td>WHMC (Lan et al., 2014)</td><td>29</td><td>1 026</td><td>5 753</td></tr><tr><td>RDMC (Ahn et al., 2013)</td><td>168</td><td>3 354</td><td>11 192</td></tr><tr><td>Proposed, p<sub>MTM</sub> = 0.1</td><td>29</td><td>329</td><td>1 235</td></tr><tr><td>Proposed, p<sub>MTM</sub> = 0.9</td><td><strong>299</strong></td><td><strong>3 561</strong></td><td><strong>16 789</strong></td></tr></tbody></table>
The proposed extended Bayesian model assessment computes one *p*-value per sensor. In the non altered case, the observations are derived from the observation model and are thus compatible with it. The smallest estimated *p*-value is $\hat{p}_n^{(T_{MC})} = 0.58$ and the largest probability of rejection $\mathbb{P}[p^{(T_{MC})} \le \alpha] = 8 \times 10^{-99} \le \delta$. The model assessment test thus does not reject any position estimation. In the altered case, most sensors are not rejected, as they are little to not affected by this observation modification. Indeed, for these unaffected sensors, the smallest estimated *p*-value is $\hat{p}_n^{(T_{MC})} = 0.60$ and the largest probability of rejection $\mathbb{P}[p^{(T_{MC})} \le \alpha] = 5 \times 10^{-159} \le \delta$. However, the model assessment test rejects the three sensors affected by the alteration. For these three sensors, the largest estimated *p*-value is $\hat{p}_n^{(T_{MC})} = 3 \times 10^{-5}$ and the smallest probability of rejection $\mathbb{P}[p^{(T_{MC})} \le \alpha] = 1 - 3 \times 10^{-6} > 1 - \delta$. Therefore, the model assessment approach detects the generated incompatibility between the observations and the observation model. In addition, it pinpoints which sensors are problematic, which helps in identifying the origin of the incompatibility.

<figure><img src="image_43.png" /><figcaption>(a) Marginals with the true observations. The model assessment test does not reject any position estimation.</figcaption></figure>

<caption>(b) Marginals with one observation divided by 10. The model assessment test rejects the three sensors affected by this observation. The other sensors are not rejected, as their positions are little to not affected by this observation modification.</caption>

Figure 5.6 – Application of model assessment on sensor localization problem.

## 5.4.2 Realistic astrophysical data

The overall approach is now applied to a synthetic yet realistic case of the general inverse problem addressed in this thesis. The goal is to reconstruct maps of physical parameters of a molecular cloud from radio wave multispectral intensity maps. Each observation map contains $N = 4096$ pixels. Each pixel is associated with $D = 4$ physical parameters $\theta = (\kappa, P_{\text{th}}, G_0, A_V^{\text{tot}})$, so that the aim is to infer a set of parameters $\Theta = (\theta_n)_{n=1}^N$ in dimension $N \times D = 16536$. The angle $\varphi$ is set to 0 deg. As $\kappa$ is a nuisance parameter related to the conditions of observations, its ground truth value is set to 1 over the whole map. The main parameters of interest are the thermal pressure $P_{\text{th}}$, the intensity of a UV radiative field $G_0$ – measured in reference to the Habing ISRF – and the visual extinction $A_V^{\text{tot}}$, related to the cloud depth along the line of sight. The physical parameters $\Theta = (\theta_n)$ undergo the preprocessing described in Chapter 4 to have similar dynamics. As an abuse of notation, the preprocessed physical parameters are also denoted $\Theta = (\theta_n)$. In the following, all distributions are defined on the preprocessed physical parameters.

Figure 5.7 shows the ground truth parameters $\Theta^*$ in original scale. These maps are chosen
according to a plausible astrophysical scenario (Pety et al., 2017). They correspond to a PDR seen edge-on, with a UV source illuminating the right of the image. This scenario covers a wide variety of physical environments such as an atomic diffuse medium on the right of the map, an actual PDR on the right, a deep molecular cloud with a buried source on the left.

<figure><img src="image_44.png" /><figcaption>Figure 5.7 – Structure of the synthetic molecular cloud and associated true maps of physical parameters $\Theta$. The circle on the top left of the maps indicates a buried source. The surrounding region is a deep molecular cloud. The region on the right is a PDR illuminated by a UV source out of the image.</figcaption></figure>

From the true (preprocessed) maps $\Theta^*$, the emulator $\tilde{f}$ of the Meudon PDR code derived in Chapter 4 generates observation maps of $L = 10$ emission lines. These lines are $^{12}$CO lines of mid-$J$ rotational transitions, from $J = 4-3$ to $J = 13-12$. For each line $\ell$, $\tilde{f}_\ell$ ranges from $10^{-18}$ to $10^{-2}$ erg cm$^{-2}$ s$^{-1}$ sr$^{-1}$. These maps are deteriorated according to the observation model from Eq. 5.1. In other words, they are affected by an Gaussian additive uncorrelated noise, a lognormal multiplicative uncorrelated noise and censorship. Model misspecification noise is not considered in this example, as the forward model $\tilde{f}$ is used both to define the observation $\mathbf{Y}$ and to infer the physical parameters $\Theta$. The standard deviation of the multiplicative noise is set to $\sigma_m = \log(1.1)$, which roughly represents a 10% alteration in average that corresponds to a calibration error. For the additive noise, $\sigma_{a,n\ell} = \sigma_a = 1.38715 \cdot 10^{-10}$ erg cm$^{-2}$ s$^{-1}$ sr$^{-1}$ for all pixel $n$ and line $\ell$ so that the signal-to-noise ratio (SNR) varies between -81 and 79 dB. The censorship level is set to $\omega_{n\ell} = \omega = 3\sigma_a$ for all pixel $n$ and line $\ell$. Figure 5.8 shows the observation maps of two lines and the spatial distribution of censorship importance.

<figure><img src="image_45.png" /><figcaption>Figure 5.8 – Some observation maps of the astrophysical experiment. From left to right: line $\ell = 1$, line $\ell = 10$, proportion of censored lines per pixel.</figcaption></figure>

The likelihood approximation is obtained as indicated in Section 5.1.1, and its parameters $a_\ell$ are adjusted as described in Appendix 5.B. The validity set $\mathcal{C}$ of physical parameters is set as in Wu et al. (2018), and the penalty parameter $\xi$ of the smooth uniform prior is set to $10^4$. Given
<figure><img src="image_46.png" /><figcaption>Figure 5.9 – Inference results: (left) ground truth $\Theta^*$; (middle) MMSE estimate from the proposed transition kernel; (right) size of the 95% credibility interval (CI), in % of the size of the validity intervals. All values of $\Theta$ are displayed in the preprocessed space.</figcaption></figure>
the smoothness of the true maps, for each d the chosen spatial regularizer h is taken as

$$h(\Theta_{\cdot d}) = \|\Delta\Theta_{\cdot d}\|_2^2 = \sum_{n=1}^{N} \sum_{i \in V_n} (\theta_{nd} - \theta_{id})^2, \qquad (5.33)$$

where Δ is the discrete 5-point based 2D Laplacian operator, and $V_n$ is the set of neighbors of pixel n induced by Δ. The hyperparameter $\tau$ from Eq. 5.11 is fixed to $\tau = (10, 2, 3, 4)$.

Inference is carried out using 10000 iterations of a Markov chain including 1500 burn-in samples. The parameters of the proposed sampler are set to $a = 0.99$, $\epsilon = 10^{-5}$ and $\eta = 5 \times 10^{-7}$ for PMALA, and to $p_{MTM} = 0.5$ and $K = 50$ for MTM. Since the operator $\Delta$ only compares a pixel to its four neighbors and since the indicator prior and likelihood are pixel-wise, the set of pixels can be partitioned into two conditionally independent subsets of pixels. A two sites Chromatic Gibbs sampling (Gonzalez et al., 2011) is therefore performed in the MTM kernel to speed up computations. Note that using the smooth uniform prior as a proposal distribution in MTM is inefficient due to the small size of high probability regions compared to the volume of $C$. The proposal distribution $q$ is based on the spatial prior (Eq. 5.33) instead. For any pixel, one can show that the conditional spatial prior is a Gaussian distribution centered on the mean of the set of neighboring pixels $V_n$. Since maps are assumed to be smooth, the likelihood functions for a pixel n and its neighbors should correspond to similar modes in the parameters' domain. If the neighbors are not all in the same mode, the mean of the neighbors will in general not fall in a high probability region. Therefore, for a pixel n, the proposal distribution is defined as a Gaussian mixture whose modes are all the means of non-empty subsets $V \in \mathcal{P}(V_n)$ of $V_n$:

$$q(\theta_n | \Theta_{\setminus n}) \propto \prod_{d=1}^{D} \sum_{V \in \mathcal{P}(V_n)} \exp \left[ -2\tau_d \sum_{i \in V} (\theta_{nd} - \theta_{id})^2 \right] \qquad (5.34)$$

$$\propto \prod_{d=1}^{D} \sum_{V \in \mathcal{P}(V_n)} \exp \left[ -2\tau_d |V| \left( \theta_{nd} - \frac{1}{|V|} \sum_{i \in V} \theta_{id} \right)^2 \right]. \qquad (5.35)$$

Performance is assessed for the minimum mean square error (MMSE) estimate $\hat{\Theta}$. Recall that the inferred parameters $\Theta$ correspond to normalized logarithms of physical parameters $\Theta$. Therefore, prediction errors on the $D$ parameter maps $\Theta_{\cdot d}$ are comparable. The quality of the reconstruction is quantified with the mean squared error (MSE) $\|\hat{\Theta} - \Theta^*\|_2^2$ and the reconstruction signal-to-noise ratio (R-SNR) $20 \log_{10} \left( \frac{\|\Theta^*\|}{\|\hat{\Theta} - \Theta^*\|} \right)$.

Figure 5.9 shows the estimation results. The MMSE estimate $\hat{\Theta}$ (middle) is very close to the ground truth $\Theta^*$ (left). The reconstructions are qualitatively very consistent with the underlying physics. The parameter $\Theta_{.4}$, corresponding to $A_V^{\text{tot}}$, is known by astrophysicists to be the most difficult to retrieve from $^{12}$CO low and mid-J lines. Indeed, these lines get optically thick, which means that observed photons come from the surface of the cloud, and that the integrated intensities saturate with $A_V^{\text{tot}} \gtrsim 7$ mag. Such pixels appear in the top left corner of the ground truth map.

Table 5.3 shows the MSE and the R-SNR for each parameter $\Theta_{\cdot d}$, and the relative size of the credibility intervals with respect to the associated (normalized) validity interval $C$. As expected, the MSE is larger for $\Theta_{.4}$ ($\leftrightarrow A_V^{\text{tot}}$), and the relative size of its credibility intervals are overall the largest, about 16.2%. The problem is also very ill-posed for all parameters in pixels with very low SNR, where most of the lines are censored, see Figure 5.8 (right). To interpret the results from an astrophysical viewpoint, performances are computed over two subsets of pixels with either less or more than 50% of censored lines. As expected, the credibility intervals of the latter are about twice as large as the former. Finally, all the parameters but $A_V^{\text{tot}}$ are well constrained for pixels with less than 50% of censored lines. The inference remains challenging since the posterior contains many local modes with high $L$ values, but the proposal distribution $q$ permits the Markov
chain to successfully reach the mode of interest. The relative quadratic error results in an R-SNR between 15.5 dB and 23.4 dB. Credibility intervals at 95% level remain small, ranging from 5.7% to 9.3% of the admissible interval C.

<table><caption>Table 5.3 – Reconstruction metrics and relative size of credible intervals for the astrophysics experiment. The R-SNR is not defined for $\Theta_{.1}$, as its ground truth is 0 everywhere.</caption><thead><tr><th rowspan="2"></th><th colspan="2">MMSE</th><th colspan="3">Mean 95% credibility intervals size</th></tr><tr><th>MSE</th><th>R-SNR (dB)</th><th colspan="2">censorship</th><th>overall</th></tr><tr><th></th><th></th><th></th><th>≤ 50%</th><th>> 50%</th><th></th></tr></thead><tbody><tr><th>Θ<sub>.1</sub></th><td>0.017</td><td>–</td><td>6.1 %</td><td>11.9 %</td><td>6.8 %</td></tr><tr><th>Θ<sub>.2</sub></th><td>0.019</td><td>16.8</td><td>9.3 %</td><td>20.6 %</td><td>9.9 %</td></tr><tr><th>Θ<sub>.3</sub></th><td>0.009</td><td>23.4</td><td>5.7 %</td><td>19.8 %</td><td>6.5 %</td></tr><tr><th>Θ<sub>.4</sub></th><td>0.034</td><td>15.5</td><td>16.3 %</td><td>14.5 %</td><td>16.2 %</td></tr></tbody></table>

Combining all the difficulties addressed in the general inverse problem addressed in this thesis, this synthetic inverse problem illustrates the good performances of the proposed approach in a challenging scenario. The proposed likelihood approximation enabled handling the censorship and mixture of noises present in the observation model. Dealing with a multimodal posterior distribution, the MTM kernel allows the different modes to be visited, while the PMALA kernel permits to explore them efficiently. The proposed sampler provides high quality estimates and informative credibility intervals.

## 5.5 Conclusion

In this chapter, we addressed a family of inverse problems that combine several difficulties: a non-linear black-box forward model, potentially non-injective, that covers multiple decades; observations damaged by both censorship and a mixture of additive and multiplicative noises. The proposed approach takes into account as many sources of uncertainty as possible. The likelihood is intractable and leads to a potentially multimodal posterior distribution. An approximation of the likelihood was proposed, based on a model reduction and an approximate parametric noise mixture model with controlled error. The prior distribution combined a spatial regularization and a smoothed uniform distribution encoding validity constraints on the physical parameters.

To efficiently sample from the resulting multimodal posterior, an original MCMC algorithm combining two kernels was proposed. The Gibbs-like MTM kernel permits jumps between modes, while the PMALA kernel efficiently explores the local geometry of each mode. The proposed sampler was shown to be competitive with state-of-the-art multimodal sampling methods on a Gaussian mixture model and a sensor localization problem. A realistic application to a challenging inverse problem on a large observation map has shown the interest and the good performances of the proposed approach. Estimation errors remain small and uncertainties are quantified.

In addition to the MCMC algorithm, the Bayesian hypothesis testing procedure was extended. Uncertainties due to the Monte Carlo evaluation of the p-value are accounted for in the rejection or non-rejection decision. This approach permits in real applications to assess the compatibility between the observations and the observation model.

In the next chapter, the proposed approach is applied to real observations of photodissociation regions.
## Appendix 5.A Tuning automatically the prior hyperparameters

The considered prior from Eq. 5.11 is reminded here:

$$ \pi(\Theta) \propto \exp \left( -\xi \tilde{\tau}_{CN}(\Theta) - \sum_{d=1}^{D} \tau_d h(\Theta_{.d}) \right). $$

This prior relies on two hyperparameters: $\xi > 0$ weighs the smooth indicator function, and $\tau$ contains the $D$ spatial regularization weights of the $D$ maps to be inferred. The indicator prior parameter $\xi$ has little impact on estimations. However, the values of the spatial regularization weights $\tau$ have a dramatic impact on the posterior distribution. In some cases, this parameter can be learned from the data using a hierarchical model, as in Galliano (2018). For instance, Pereyra et al. (2015) proposes a general approach to estimate the best regularization parameter along with the physical parameters $\Theta$. However, this approach requires the prior to be k-homogeneous. In our case, the box constraints encoding the validity intervals violate this condition. Using hierarchical model and sampling from the extended posterior distribution is thus much harder.

We tried applying the maximum marginal likelihood estimator presented in Vidal et al. (2020) to automatically tune $\tau$. However, the additive noise variance is sometimes overestimated in astrophysics applications, as we will show in Chapter 6. In such cases, the spatial regularization weights $\tau_d$ diverge, leading to constant estimated maps. Such maps are unrealistic for the environments considered in Chapter 6. To avoid such divergences, we chose to set manually the model hyperparameters, depending on the expected smoothness in the final estimated maps and with inversions with different values.

## Appendix 5.B Optimization of the approximation parameter

For each channel $\ell$, the parameter $a_\ell = (a_{\ell,0}, a_{\ell,1})$ locates the frontiers between low, intermediate and high values regimes of $\ln \tilde{f}_\ell$ in the definition of $\lambda$ (Eq. 5.8). It has a critical influence on the approximation quality. It should be adjusted to $\ln \tilde{f}_\ell$, $\sigma_a$ and $\sigma_m$. For simplicity, in this subsection, likelihood functions are conditioned with respect to $z = \ln \tilde{f}_\ell(\theta) \in \mathbb{R}$ instead of $\theta \in \mathbb{R}^D$. The true likelihood is not explicit, but the model in Eq. 5.1 can be easily sampled from, and the approximation (Eq. 5.9) is known.

The parameter $a_\ell$ is set to obtain an approximation as close as possible to the true likelihood, with respect to some divergence criterion. The Kullback-Leibler (KL) divergence would be a natural choice. However, due to the number of decades spanned, the standard deviation of KL estimators is in practice larger than the quantity of interest (Kraskov et al., 2004), which prevents from performing optimization. The Kolmogorov-Smirnov (KS) distance is not affected by this property: for a given $z$, it only requires ordered samples $(y^{(i)})_{i=1}^M$. It reads

$$ \hat{D}_{\mathrm{KS}}(z, \mathbf{a}_{\ell}) = \sup_{y \in \mathbb{R}} |\hat{F}_{M}(y|z) - \tilde{F}(y|z, \mathbf{a}_{\ell})|, \qquad (5.36) $$

where $\hat{F}_M(\cdot|z)$ is the empirical cdf of the true likelihood $\pi(\cdot|z)$ estimated from $M$ samples $y^{(i)}$, and $\tilde{F}(\cdot|z, \mathbf{a}_\ell)$ is the cdf of the proposed approximation (Eq. 5.9). Assuming that $\theta$ follows a uniform distribution on $C$ yields a distribution on $z$ with pdf $\pi(z)$ which can be estimated by kernel density estimation (KDE). The function to minimize is

$$ \varphi(\mathbf{a}_\ell) = \mathbb{E}_z [\hat{D}_{\mathrm{KS}}(z, \mathbf{a}_\ell)] = \int \hat{D}_{\mathrm{KS}}(z, \mathbf{a}_\ell)\pi(z)dz. \qquad (5.37) $$

An estimator $\hat{\varphi}$ can be obtained using numerical integration on $z$ over $S$ bins. The higher $M$ and $S$, the better the estimation accuracy. Minimizing $\hat{\varphi}$ can be performed using a grid search, which is quite computationally intensive. A cheaper alternative is to use a Bayesian optimization (BO) procedure (Shahriari et al., 2016). This optimization was applied for each channel in the
astrophysical application described in Section 5.4.2. Both grid search and BO approaches were used. The KDE of $\pi(z)$ was performed from 810 000 samples. The BO procedure was run with $S = 100$ and $M = 250 000$ using Nogueira (2014–) with default parameters. Figure 5.10 shows the results for one channel. The proposed approximation with adjusted $a_\ell$ is closer to the true likelihood than a purely additive Gaussian approximation, i.e., $a_{\ell,0} > \max_j z^{(j)}$, or a purely multiplicative lognormal approximations, i.e., $a_{\ell,1} < \min_j z^{(j)}$.

<figure><img src="image_47.png" /><figcaption>Figure 5.10 – Maximization of $-\log_{10} \hat{\varphi}$ using both Bayesian optimization (BO) and grid search for one channel of the astrophysical case detailed in 5.4.2. In BO, a Gaussian process (GP) replaces the function to optimize (left column). The red dashed vertical bar represents the value of $a_0+a_1/2$ for which the additive and multiplicative noises have equal variances, i.e. $\sigma_a^2 = \tilde{f}_\ell(\theta)^2 \text{Var}[\epsilon_{n,\ell}^{(m)}]$, at $\lambda = \frac{1}{2}$. For clarity, all scales are displayed in $\log_{10}$ scale, while computations are done in $\ln$ scale.</figcaption></figure>

## Appendix 5.C Sampling from the smoothed indicator distribution

This section describes the algorithm to draw samples from the real-valued probability distribution with density $\pi(\theta) \propto \exp(-\xi \tilde{\iota}_{[l,u]}(\theta))$, with $l < u$ and $\tilde{\iota}_{[l,u]}$ introduced in Eq. 5.10. To this aim, consider the generalized normal distribution $GN(0, 1/\xi^4, 4)$ of pdf (Nadarajah, 2005)

$$p_{GN}(\theta) = \frac{2\xi^{\frac{1}{4}}}{\Gamma(1/4)} \exp(-\xi \theta^4). \qquad (5.38)$$
Note that $\pi(\theta)$ is a continuous extension of a uniform distribution and of this generalized normal distribution at 0.

$$ \pi(\theta) \propto \begin{cases} p_{GN}(\theta - l) & \text{if } \theta < l, \\ p_{GN}(0) & \text{if } \theta \in [l, u], \\ p_{GN}(u - \theta) & \text{if } \theta > u. \end{cases} \qquad (5.39) $$

The normalizing constant of $\pi(\theta)$ is $1 + p_{GN}(0)(u - l)$. The weight of the uniform section in the combination is therefore

$$ w_{\text{Unif}} = \frac{1}{1 + \frac{\Gamma(1/4)}{2} \frac{1}{\xi^{\frac{1}{4}}(u-l)}} \qquad (5.40) $$

Algorithm 5.5 summarizes the procedure to sample from $\pi(\theta)$.

**Algorithm 5.5:** Sampling from the smooth distribution in Eq. 5.39

**Input:** scale factor $\xi$, bounds $l, u \in \mathbb{R}$ such that $l < u$

**Output:** sample $\theta$

1 $w_{\text{Unif}}$ // using Eq. 5.40

2 $z \sim B(w_{\text{Unif}})$

3 **if** $z = 1$ **then** $\theta \sim \text{Unif}(l, u)$

4 **else**
5 $\theta \sim GN(0, 1/\xi^4, 4)$ // using Nardon and Pianca (2009)

6 $\left\lfloor \text{if } \theta < 0 \text{ then } \theta = \theta + l \text{ else } \theta = \theta + u \right\rfloor$
# Chapter 6

# Application to real data

> "In theory, there is no difference between theory and practice. In practice there is."

> Benjamin Brewster, "The Yale Literary Magazine", February 1882

## Contents

<table><tr><td>6.1 Summary of the inversion procedure</td><td>142</td></tr><tr><td>6.2 NGC 7023</td><td>144</td></tr><tr><td>6.3 The Carina nebula</td><td>148</td></tr><tr><td>6.4 Orion molecular cloud 1 (OMC-1)</td><td>156</td></tr><tr><td>6.5 Conclusion</td><td>166</td></tr><tr><td>Appendix 6.A Orion Bar</td><td>167</td></tr></table>

In Chapter 4, we derived an approximation of the Meudon PDR code, that we now use to model photodissociation regions (PDRs). Chapter 5 presented the full observation model, prior and posterior distributions, and the proposed method to sample from the posterior. It also introduced the Bayesian test of hypothesis that permits to assess the compatibility of the model with the observations. In this Chapter, we apply the full inference workflow to real observations of photodissociation regions in star forming regions in order to determine the physical conditions in these environments. A journal article is currently in preparation.

Section 6.1 summarizes the full inversion procedure. In Section 6.2, we apply our method to observations of NGC 7023 ($N = 1$ and $L = 17$) studied in Joblin et al. (2018). Appendix 6.A describes a similar analysis on the Orion bar, that was also studied in Joblin et al. (2018). Section 6.3 presents a first application to real-life multi-pixel observation maps. It studies the Carina nebula maps analyzed in Wu et al. (2018) ($N = 176$ pixels and $L = 12$ lines for inversion). These first two analyses demonstrate that our results are consistent with those of Joblin et al. (2018) and Wu et al. (2018). They also show that our approach is richer as it provides complete uncertainty quantification in addition to the point estimates. Section 6.4 studies the OMC-1 observations introduced in Goicoechea et al. (2019) ($N = 2475$ pixels and $L = 4$ lines for inversion). The proposed analysis is the first for these OMC-1 observations.
## 6.1 Summary of the inversion procedure

This section summarizes the inversion procedure presented in Chapters 4 and 5. Figure 6.1 illustrates the full procedure.

<figure><img src="image_48.png" /><figcaption>Figure 6.1 – Full inference workflow. Observations and noise standard deviations are assumed known. The indicated model and sampler parameters need to be adjusted for each application.</figcaption></figure>

**Setting the posterior distribution** – In this chapter, the posterior distribution is defined on maps of physical parameters $\Theta = (\theta_n)_{n=1}^N \in \mathbb{R}^{N \times D}$. The physical parameter vectors contain $D=4$ elements: $\theta = (\kappa, P_{\text{th}}, G_0, A_V^{\text{tot}})$, with $\kappa$ a scaling parameter, $P_{\text{th}}$ the cloud thermal pressure, $G_0$ the intensity of the incident radiative field and $A_V^{\text{tot}}$ the cloud total visual extinction – see Chapter 1 (Section 1.1.4) for a description of these parameters. The posterior distribution relies on a likelihood function and on a prior distribution.

The likelihood function relies on a forward model and a noise model. The noise model described in Chapter 3 (Section 3.4) combines a Gaussian additive noise and a lognormal multiplicative noise. The corresponding standard deviations, $(\sigma_{a,n\ell})$ and $\sigma_m$, are assumed known. The forward model $\mathbf{f}$ is the Meudon PDR code, which requires a few hours for each evaluation. As inference requires many evaluations, we replace it with a fast and accurate neural network-based emulator $\tilde{\mathbf{f}}$ (Chapter 4). The resulting likelihood function has no simple closed-form expression. Chapter 5 (Section 5.1.1) proposed a simple approximation with controlled error. This approximation relies on parameters $(a_\ell)$ to be adjusted to minimize the error – see Chapter 5 (Appendix 5.B).

The prior distribution combines an indicator term that restricts values to an acceptable range for each parameter, and a spatial regularization term that favors smooth reconstructed maps for
multi-pixel observations. In this chapter, the indicator prior is set on the validity intervals used to train the emulator $\tilde{f}$ in Chapter 4, i.e., $P_{\text{th}} \in [10^5, 10^9]$ K cm$^{-3}$, $G_0 \in [1, 10^5]$ $^1$, and $A_V^{\text{tot}} \in [1, 40]$ mag. The scaling parameter $\kappa$ is limited to $[0.1, 10]$. The smoothing of the indicator prior at the edges of the domain is performed using $\xi = 10^4$. The spatial regularization term relies on a weight parameter $\tau \in \mathbb{R}^4$, that will be adjusted for multi-pixel observations.

**Performing inference with the proposed sampler** – Markov chain Monte Carlo (MCMC) algorithms produce sets of $T_{\text{MC}} \ge 1$ correlated samples $(\Theta^{(t)})_{t=1}^{T_{\text{MC}}}$ from the posterior distribution. In each application, the number of iterates $T_{\text{MC}}$ is chosen so that the posterior distribution is well sampled. In particular, the duration of the burn-in phase $T_{\text{BI}}$ is set so that estimators are only evaluated from iterates in the Markov chain stationary phase.

The MCMC proposed in Chapter 5 (Section 5.2) combines two update steps. The first update step is called PMALA. It is an MCMC variant of a preconditioned gradient descent that performs efficient local exploration of the posterior distribution. We resort to the RMSProp diagonal preconditioner, which offers a good trade-off between computational costs and preconditioning efficiency – see Chapter 2 (Section 2.2.2.5). It relies on three parameters: a step size $\eta$, a damping parameter $\epsilon$ and an exponential decay rate $a$. In this chapter, the damping parameter and exponential decay rate are set to their default values, i.e., $\epsilon = 10^{-5}$ and $a = 0.99$. In each application, the step size is chosen so that the average acceptance probability is close to the MALA optimum acceptance probability, i.e., around 50%-60%.

The second update step is called MTM. It permits to globally explore the posterior and to escape local minima for multimodal distributions. It handles pixels individually. For each pixel, it generates $K \ge 1$ candidates from a proposal distribution $q$, and selects one to perform an accept-reject step. At each iteration $t$, one of the two update steps is selected with probability $p_{\text{MTM}} \in [0, 1]$ for MTM and $1-p_{\text{MTM}}$ for PMALA. For each application, the number of candidates $K$ and the selection probability $p_{\text{MTM}}$ are set to explore efficiently the posterior distribution, i.e., to obtain high effective sample sizes (ESSs) – see Chapter 2 (Section 2.2.2.2).

The output correlated set $(\Theta^{(t)})_{t=1}^{T_{\text{MC}}}$ can be used to perform checks including on the chain convergence – did the sampling reach high probability region of the posterior distribution? – and on the mixing performance – did the sampling explore well the posterior distribution? These tests were performed for all the applications presented in this chapter, but are not shown for conciseness.

**Bayesian hypothesis testing** – The ability of the Meudon PDR code and of the noise model to explain the observations is checked for each application. We resort to the model assessment method presented in Chapter 5 (Section 5.3). This method evaluates a $p$-value for each pixel $n$. A $p$-value below a confidence level $\alpha$ pinpoints an incompatibility between the model and the observation. The proposed sampling algorithm evaluates Monte Carlo (MC) estimators $\hat{p}_n^{(T_{\text{MC}})}$ of the $p$-values – see Chapter 2 (Section 2.2.2) for a description on MC estimators. The uncertainty inherent to the MC estimator is accounted for, allowing us to evaluate a probability of rejection $\mathbb{P}[p_n \le \alpha]$. We proposed a three-case rule by introducing a probability threshold $\delta$ to make the test robust to the MC estimator error. In this chapter, we choose $\alpha = 0.05$ and $\delta = 0.1$.

**MMSE, credibility intervals and uncertainty factor** – The set of correlated samples $(\Theta^{(t)})_{t=1}^{T_{\text{MC}}}$ are used to evaluate MC estimators. In this chapter, two estimators are considered: the minimum mean square error (MMSE) and credibility intervals (CIs). The MMSE $\hat{\Theta}_{\text{MMSE}}$ is the sample mean and approximates the posterior expectation. Credibility intervals enable to quantify the uncertainty associated with each inferred physical parameter. See Chapter 2 (Section 2.1.2) for a description of these estimators.

The parameters inferred in this chapter cover multiple decades. Uncertainties on these parameters are thus better described with a multiplicative error. We define the uncertainty factor (UF)

<footnote n="1">1. In this thesis, $G_0$ is defined in reference to the Habing ISRF. See Chapter 1 (Section 1.1.4) for more details.</footnote>
to quantify this multiplicative error. For a confidence level $\alpha$ and a credibility interval $[l_\alpha, u_\alpha]$,

$$UF_{\alpha} = \sqrt{\frac{u_{\alpha}}{l_{\alpha}}}, \qquad (6.1)$$

so that the true value is expected to be within a factor $UF_\alpha$ below or above the estimated value. For a lognormal posterior distribution on a physical parameter, the $UF_{68\%}$ corresponds to a $1\sigma$ error. Similarly, $UF_{95\%}$ corresponds to a $\simeq 2\sigma$ error. By default, we denote $UF = UF_{68\%}$.

## 6.2 NGC 7023

Joblin et al. (2018) analyzed two one-pixel PDR observations: one of NGC 7023, and one of the Orion Bar. We applied the proposed inversion procedure to both. In this section, we only present results on NGC 7023 for conciseness. We demonstrate that the proposed inversion procedure yields consistent results with Joblin et al. (2018) and provides complete uncertainty quantification in addition to point estimates. Appendix 6.A presents the results for the Orion Bar.

Located in the Cepheus constellation, NGC 7023 is a PDR illuminated by the spectroscopic binary system HD 200775 [RA(2000) = 21h01m36.9s; Dec(2000) = +68 09047.800]. Its distance from the Sun was estimated at 320 ± 51 pc in Benisty et al. (2013). NGC 7023 is extensively studied to understand the effect of radiative feedback because of its edge-on geometry, brightness and proximity. The two stars in the binary system are classified as B3Ve and B5 (Alebian et al., 2008). Chokshi et al. (1988) estimated an intensity of $G_0 = 2600$ and a proton density of $n_H \simeq 4 \times 10^3 \text{ cm}^{-3}$ from the C$^+$ 158 µm and O 63 µm lines.

Joblin et al. (2018) studied the impact of radiative feedback in NGC 7023 from Herschel observations combined with some Spitzer-IRS observations for rotational lines of H$_2$. The resulting observations gathered $L = 17$ lines on $N = 1$ pixel, including $^{12}$CO lines (from $J = 11 \rightarrow 10$ to $J = 19 \rightarrow 18$), rotational H$_2$ lines (from S(0) to S(5)) and low level CH$^+$ rotational lines (from $J = 1 \rightarrow 0$ to $J = 3 \rightarrow 2$). These observations were modeled with a Gaussian additive and uncorrelated noise model. The standard deviations $(\sigma_{a,\ell})_{\ell=1}^L$ were assumed known and included calibration error. The analysis was performed with the version 1.5.4 of the Meudon PDR code with a fixed value of $G_0 = 2600$. As the observed lines only trace the warm molecular layer of the PDR, the visual extinction $A_V^{\text{tot}}$ was fixed to 10 mag. The observation angle $\varphi$ was set to 60 deg to account for the edge-on geometry of the cloud – the Meudon PDR code cannot run for $\varphi = 90$ deg due to its infinite slab with finite thickness geometry, and $\varphi = 60$ deg is the closest to an edge-on inclination. The inference was performed on $\theta = (\kappa, P_{\text{th}})$ using a grid search on $P_{\text{th}}$ and a simple continuous optimization for $\kappa$. They obtained $P_{\text{th}} = 1 \times 10^8 \text{ K cm}^{-3}$ and $\kappa = 0.7$. This high $P_{\text{th}}$ value indicates that the radiative feedback compresses the PDR.

**Inversion setup** – The posterior distribution is defined on the physical parameter vector $\theta = (\kappa, P_{\text{th}}, G_0, A_V^{\text{tot}}) \in \mathbb{R}^4$. As in Joblin et al. (2018), we set the observation angle to $\varphi = 60$ deg. The calibration error is included in the additive error. The multiplicative noise source accounts only for the model misspecification. We set $\sigma_m = \ln 1.3$, where a $1\sigma$ error for this multiplicative term corresponds to a factor 1.3. The resulting observation model is

$$\forall \ell \in [1, L], \quad y_{\ell} = \varepsilon_{\ell}^{(m)} \tilde{f}_{\ell}(\theta_{\ell}) + \varepsilon_{\ell}^{(a)}, \qquad (6.2)$$

with $\varepsilon_{\ell}^{(m)} \sim \text{Lognormal}(-\sigma_m^2/2, \sigma_m^2)$ and $\varepsilon_{\ell}^{(a)} \sim \mathcal{N}(0, \sigma_{a,\ell}^2)$ for all lines $\ell$. The likelihood approximation parameters $a_\ell$ are optimized using the Bayesian optimization procedure described in Chapter 5 (Appendix 5.B). The proposed sampler is run for $T_{\text{MC}} = 20000$ iterations including $T_{\text{BI}} = 500$ of burn-in. The MTM kernel has a selection probability of $p_{\text{MTM}} = 0.5$. As the observation only contains one pixel, the proposal $q$ is set to the smoothed uniform prior distribution and the number of candidates to $K = 2000$. The PMALA step size is set to $\eta = 0.05$.
**Model assessment and Bayesian p-value** – The proposed sampling algorithm led to an estimated p-value of $\hat{p}^{(T_{\text{MC}})} = 0.56$ and to a rejection probability of $\mathbb{P}[p^{(t)} \le \alpha] < 10^{-200} < \delta$, which does not lead to a rejection. This indicates that the Meudon PDR code and the noise model are compatible with the observations.

Figure 6.2 shows how the observations y compare to the marginal posterior predictive dis- tributions $\pi(\tilde{y}_{nl}|\mathbf{y})$ for lines used and not used for the inversion. It permits to visualize the compatibility between the observation y and the posterior distribution on the physical parameters $\theta$. However, these plots are imperfect as the Gaussian additive noise can lead to negative predicted observations $\tilde{y}_\ell$ in the low SNR regime. These negative values are not displayed. Underestimating the integrated intensity $y_\ell$ may thus not imply incompatibility between an observation $y_\ell$ and the observation model.

Figure 6.2a shows that for the L lines used for the inversion, the observations $y_\ell$ fall in high probability regions of the marginal predictive distributions $\pi(\tilde{y}_{nl}|\mathbf{y})$, in agreement with the results of Joblin et al. (2018). In particular, the Meudon PDR prediction from the MMSE $\hat{f}(\hat{\theta}_{\text{MMSE}})$ successfully reproduces all the lines at once. For lines that were not used in the inversion procedure, Figure 6.2b shows that this inversion reproduces well the low J lines of $^{12}\text{CO}$, but underestimates the $^{13}\text{CO}$ and $\text{C}^{18}\text{O}$. This discrepancy was already noted in Joblin et al. (2018). It might be due to the fact that the considered grid of Meudon PDR code simulations does not include mutual radiative shielding between $^{12}\text{CO}$ and its isotopologues. Figure 6.2c shows satisfying compatibility for the marginals except for the $O\,3p\,J = 1-2$ at 63 µm, which is overestimated by a factor $\approx 10$. Joblin et al. (2018) showed the same overestimation for the $O\,63\,\mu\text{m}$ line and underestimations for $^{13}\text{CO}$ and $\text{C}^{18}\text{O}$ lines.

**Inference results** – Figure 6.3 shows the pairwise histograms of the posterior samples obtained with the proposed MCMC algorithm. We first note that posterior distribution indicates high uncertainties as it covers a large portion of the physical parameter space. The visual extinction $A_V^{\text{tot}}$ is not constrained at all as the full interval $[1, 40]$ mag is covered. This is due to the fact that the observed lines are warm gas tracers and are thus emitted in the warm external layer of the PDR. Therefore, they do not trace the total column density of matter. However, the triplet $(\kappa, P_{\text{th}}, G_0)$ is well constrained, although the individual credibility intervals are large. In particular, $\kappa$ and $G_0$ are strongly anti-correlated. The $G_0$ thus has a large uncertainty. However, its MMSE and posterior mode are of a few $10^3$, which is compatible with past estimations. The thermal pressure $P_{\text{th}}$ is well constrained at $10^8 \, \text{K cm}^{-3}$.

Joblin et al. (2018) identified a positive correlation among multiple sources between estimated $P_{\text{th}}$ and $G_0$. This correlation could have been caused by a positive correlation in uncertainties on individual ($P_{\text{th}}, G_0$) estimations. Figure 6.3 shows a negative slope in the ($P_{\text{th}}, G_0$) joint histogram, which indicates that the correlation in uncertainties is negative. Figure 6.21 in Appendix 6.A shows that a similar inference on Orion bar observations leads to the same negative correlation in uncertainties. Therefore, the positive correlation among multiple sources should have a physical origin, as hypothesized in Joblin et al. (2018).

Table 6.1 shows the physical parameters estimated values and their credibility intervals. For $P_{\text{th}}, G_0$ and $A_V^{\text{tot}}$, the values from Joblin et al. (2018) fall in the credibility intervals, which means that the two estimations are consistent. In addition, the MMSE is close to the values from Joblin et al. (2018). However, this is mostly a coincidence considering how spread the posterior distribution is over the parameter space. The scaling parameter $\kappa$ is the only parameter for which the two estimations are incompatible. In principle, $\kappa$ corresponds to observational effects such as the beam filling factor and the inclination of the PDR with respect to the line of sight. However, the forward models differ in the two estimations. Joblin et al. (2018) relies on the version 1.5.4 of the Meudon PDR code, while we use an emulator of version 1.7 which contains better implementation of physical processes and updated atomic and molecular data. Besides, Joblin et al. (2018) exploited dust grain properties corresponding to dense gas, while we used average galactic grain properties. These differences appear to result in a multiplicative
<figure><img src="image_49.png" /><figcaption>(a) Lines used for the inversion.</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>(b) Additional lines of CO and isotopologues, unseen during the inversion.</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>(c) Additional lines of other molecules, atoms and ions, unseen during the inversion.</figcaption></figure>

Figure 6.2 – Posterior predictive assessment for NGC 7023. Comparison of observations **Y** and associated noise model with posterior predictive distributions on $\tilde{f}(\Theta)$, with $\Theta \sim \pi(\Theta|\mathbf{Y})$.
factor for the lines used in the inversion that was compensated for by $\kappa$.

<figure><img src="image_52.png" /><figcaption>Figure 6.3 – Inference results for NGC 7023. Two-dimensional marginal histograms in the physical parameters $\Theta$ space. All histograms are in logarithmic norm.</figcaption></figure>
<table><caption>Table 6.1 – Inference results on NGC 7023. Joblin et al. (2018) only infers $\kappa$ and $P_{th}$.</caption><thead><tr><th></th><th></th><th>κ<br>–</th><th>P<sub>th</sub><br>(K cm<sup>-3</sup>)</th><th>G<sub>0</sub><br>–</th><th>A<sup>tot</sup><sub>V</sub><br>(mag)</th></tr></thead><tbody><tr><td rowspan="2">(Joblin et al., 2018)<br>MMSE</td><td></td><td>0.7</td><td>1 × 10<sup>8</sup></td><td>2.6 × 10<sup>3</sup></td><td>10</td></tr><tr><td></td><td>2.9</td><td>9.7 × 10<sup>7</sup></td><td>4.0 × 10<sup>3</sup></td><td>7.3</td></tr><tr><td rowspan="3">68% credibility interval</td><td>lower bound l<sub>68%</sub></td><td>1.2</td><td>8.1 × 10<sup>7</sup></td><td>1.8 × 10<sup>3</sup></td><td>2.2</td></tr><tr><td>upper bound u<sub>68%</sub></td><td>6.0</td><td>1.2 × 10<sup>8</sup></td><td>1.0 × 10<sup>4</sup></td><td>24.9</td></tr><tr><td>UF<sub>68%</sub></td><td>2.2</td><td>1.2</td><td>2.4</td><td>3.4</td></tr><tr><td rowspan="3">95% credibility interval</td><td>lower bound l<sub>95%</sub></td><td>0.8</td><td>4.7 × 10<sup>7</sup></td><td>1.3 × 10<sup>3</sup></td><td>1.1</td></tr><tr><td>upper bound u<sub>95%</sub></td><td>9.0</td><td>1.5 × 10<sup>8</sup></td><td>8.9 × 10<sup>4</sup></td><td>41.2</td></tr><tr><td>UF<sub>95%</sub></td><td>3.3</td><td>1.8</td><td>8.2</td><td>6.1</td></tr></tbody></table>

## 6.3 The Carina nebula

Having tested the method on known single-pixel observations, we now move to multi-pixel maps. In this section, we perform inference on an already studied observation of the Carina nebula. This analysis aims at completing the validation of the proposed method, at showing the relevance of the spatial regularization and of the sampling approach compared to optimization.

The Carina nebula is located in the Milky Way, at an estimated distance of 2.4 kpc from the Sun (Smith, 2006). Illuminated by the massive star clusters Trumpler 14 and Trumpler 16, it is the brightest nebula in the sky of the Southern hemisphere. It includes two major components, Car I and Car II, separated by ≈ 7 pc. The structure of Car I is farther divided into three bright regions: Car I-W, Car I-E, and Car I-S, located in the west, east and south of Car I, respectively (Whiteoak, 1994). Wu et al. (2018) studied the CO and C emission observed by Herschel SPIRE/FTS. Figure 6.4 shows the field of view, which includes Car I-E, Car I-S, and a region at the intersection of Car I and Car II (Car I/II).

<figure><img src="image_53.png" /><figcaption>(a) Localization of the field of view with respect to the Trumpler 14 star cluster.</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>(b) Effective field of view (gray pixels). The background image is taken in the DSS2-red band. The magenta crosses indicate the pointings used during observation. The blue and yellow crosses pinpoint two massive stars. Taken from Wu et al. (2018).</figcaption></figure>

<caption>Figure 6.4 – Observation of the Carina nebula.</caption>
The $G_0$ in the Car I-E region was estimated of the order of $10^4$ from stellar composition and FIR observations (Brooks et al., 2003; Mizutani et al., 2004). However, estimations with PDR models led to lower values: $G_0 = 1390$ according to Oberst et al. (2011) and $G_0 = 3200$ according to Kramer et al. (2008).

In Wu et al. (2018), the authors inferred multiple parameters – including $G_0$ – from $L = 12$ lines, including $10^{12}$CO lines (from $J = 4 \rightarrow 3$ to $J = 13 \rightarrow 12$) and 2 atomic C lines (3p, $J = 1 \rightarrow 0$ and $J = 2 \rightarrow 1$). They used a RBF interpolation of a grid of the version 1.5.4 of the Meudon PDR model with PAHs in the dust population. The noise model involved an additive, uncorrelated Gaussian noise. The prior distribution relied on an indicator term with wide validity intervals on the physical parameters and on an ad hoc physics-informed regularization term. This ad hoc regularization term enforced the reconstructed integrated intensities to be decreasing for the CO lines $J = 11 \rightarrow 10$, $J = 12 \rightarrow 11$, and $J = 13 \rightarrow 12$, as this condition seems always satisfied in observations. This ad hoc constraint is not generalizable to observations of other lines that the high $J$ rotational lines of $^{12}$CO. The inversion estimated the maximum a posteriori (MAP) using a preconditioned gradient descent algorithm with the limited memory BFGS preconditioner. In the Car I-E region, they obtained $G_0 \approx 2 \times 10^4$.

**Inversion setup** – The inference is performed with the same $L = 12$ lines as in Wu et al. (2018). Figure 6.5 shows the integrated intensity maps $\mathbf{Y}$. These maps contain $N = 176$ pixels.

<figure><img src="image_55.png" /><figcaption>Figure 6.5 – Observations of the Carina nebula. The Car I-E is brightest for all lines, as it is directly illuminated by the Trumpler 14 star cluster.</figcaption></figure>

The inference is performed with the data used in Wu et al. (2018). The standard deviations $\sigma_{a,n\ell}$ of the additive noise already include multiplicative errors. As the signal-to-noise ratio is low in this map, we decided not to include additional multiplicative noise in the observation model. Besides, some standard deviations $\sigma_{a,n\ell}$ are missing. To avoid biasing the reconstructed maps, the associated observations $y_{n\ell}$ are disregarded for the inversion $^2$. The resulting observation model is

$$y_{n\ell} = \tilde{f}_{\ell}(\theta_n) + \varepsilon_{n\ell}^{(a)}, \quad \varepsilon_{n\ell}^{(a)} \sim \mathcal{N}(0, \sigma_{a,n\ell}^2). \qquad (6.3)$$

<footnote n="2">2. The associated observations are disregarded by setting $\sigma_{a,n\ell} = 1 \text{ erg cm}^{-2} \text{ s}^{-1} \text{ sr}^{-1}$, an arbitrarily high value.</footnote>
The spatial regularization parameter is set to $\tau_d = 1$ for each of the $D = 4$ physical parameters. We tried multiple values between 0.1 and 10 for each parameter, and selected the value that led to physically consistent maps that reproduce the structures in the observations. Note that we do not resort to the same prior as in Wu et al. (2018): the article used a line-specific pixel-wise ad hoc constraint on the predicted integrated intensities, while we rely on a more generally applicable spatial regularization term.

The observation angle $\varphi$ varies in the field of view, as the Car I-E is considered to be edge-on while the Car I-S is face-on. For optically thin lines such as excited $^{12}$CO lines, $\varphi$ is degenerated with the scaling parameter $\kappa$. Therefore, we chose not to infer it, and set it to $\varphi = 0$ deg.

As this use case is the first map of real observations, we compare the maximum likelihood estimator (MLE) and the MMSE. The MLE is evaluated for comparison, as a naive estimator. It is the estimator most compatible with the observations for each pixel, and is thus very sensitive to noise. In particular, it does not exploit the spatial regularization and is not accompanied by an uncertainty quantification. The goal is to demonstrate the interest of the spatial regularization.

The MLE is evaluated using the optimization algorithm adapted from the proposed sampler mentioned in Chapter 5 (Section 5.2.3), run for 500 iterations. This variant combines two update steps. The first update step is a preconditioned gradient descent step with the RMSProp preconditioner. It was run with the step size $\eta = 0.03$. At each step $t$, it was selected with probability $1 - p_{MTM} = 0.99$. The second update step, selected with probability $p_{MTM} = 0.01$, is a simulated annealing-like variant of the MTM kernel. Unlike gradient descent algorithms, this second update step enables to escape from local minima. As the MLE does not include the spatial regularization, the proposal $q$ is set to the smoothed uniform prior. It generates $K = 200$ candidates.

The MMSE and the credibility intervals are obtained with the approach we proposed in this thesis. Both are evaluated from a Markov chain with $T_{MC} = 10000$ iterates, including 500 of burn-in. The PMALA kernel relies on a step size $\eta = 0.003$. The MTM kernel is selected with probability $p_{MTM} = 0.5$. It generates $K = 20$ candidates from the proposal $q$ described in Chapter 5 (Eq. 5.35).

**Model assessment and Bayesian p-value** – Figure 6.6 shows the maps of estimated p-values $\hat{p}_n^{(T_{MC})}$, of rejection probabilities $\mathbb{P}[p_n \le \alpha]$ and of the resulting rejection decision. All the estimated p-values are above 0.5. This indicates a good fit between the observations and the observation model. When accounting for the uncertainty on the p-values estimation, the highest model rejection probability is $4 \times 10^{-81}$. None of the pixels was rejected.

<figure><img src="image_56.png" /><figcaption>Figure 6.6 – Model assessment on the Carina observations from Wu et al. (2018).</figcaption></figure>

Figure 6.7 compares the observations with the marginal posterior predictive distributions for the two highlighted pixels. For the Car I-E pixel, the bright $^{12}$CO lines are well reproduced. However, the predicted integrated intensities for the atomic C lines differ from the true observations by up to a factor 5. Yet, these values for the C lines are considered compatible as the standard deviations of the additive noise are large. In the Car I-S pixel, the $^{12}$CO lines are not as bright
as in the Car I-E, and are roughly as bright as atomic C lines. The errors between predicted integrated intensities and observations are thus more uniformly distributed among $^{12}$CO lines and atomic C lines.

<figure><img src="image_57.png" /><figcaption>Figure 6.7 – Posterior predictive assessment two pixels of the Carina nebula on the lines used for the inversion. Comparison of observations *Y* and associated noise model with posterior predictive distributions on $\tilde{f}(\Theta)$, with $\Theta \sim \pi(\Theta|Y)$. Note that these plots are imperfect as the Gaussian additive noise can lead to negative predicted observations $\tilde{y}_\ell$ in the low SNR regime. These negative values are not displayed. For the lower intensities, underestimating the integrated intensity $y_\ell$ may thus be compatible with the additive Gaussian noise model.</figcaption></figure>

**Inference results: MLE, MMSE and UF** – Figure 6.8 shows the obtained MLE, MMSE and UF. It also shows the estimated thermal pressure $P_{\text{th}}$ and incident radiative field intensity $G_0$ maps from Wu et al. (2018). The maps of scaling factor $\kappa$ and of visual extinction $A_V^{\text{tot}}$ were not displayed in the article.

The MLE is visually not good. It presents unrealistic values and variations especially in $P_{\text{th}}$ and $A_V^{\text{tot}}$ because of its sensitivity to noise. Therefore, it cannot be exploited to extract meaningful information on the physical parameters in the observed region.
Conversely, the MMSE is smooth and spatially consistent. The maps of $P_{\text{th}}$ and $G_0$ show a clear frontier between the Car I-E and Car I-S regions, with an increase by a factor ~3–5 from Car I-S to Car I-E. They also show a blurry frontier between the Car I-S and Car I/II regions, with higher $P_{\text{th}}$ and $G_0$ in the Car I/II region. These lower values in the Car I-S region are due to the fact that this region is farther in the back than the other two, and thus farther away from the Trumpler 14 star cluster. The scaling parameter $\kappa$ ranges from 0.1 to 1.2, with its highest values attained in the Car I-E region. The lower values may correspond to smaller and spatially unresolved PDR edges. The higher values may be due to the more edge-on geometry and to a higher beam filling factor. Finally, the observed lines do not trace the total visual extinction for deep clouds as they are emitted mostly from the warm surface of the cloud. It is therefore expected not to recover well the spatial structure in the $A_V^{\text{tot}}$ map. Figures 6.10 and 6.9 show pairwise two-dimensional histograms of the posterior distribution samples. Unlike in NGC 7023, there is no strong degeneracy between $\kappa$ and $G_0$, although there is a negative correlation between the two parameters. In both pixels, there is also a negative correlation between $P_{\text{th}}$ and $G_0$. In the bright Car I-E region, this correlation is higher, which might be due to the higher SNR and constraining power of the observed lines.

The UF maps show that the uncertainty is on average lower in the bright Car I-E region. This is due to the larger signal-to-noise ratio in this region. It also illustrates the constraining power of the spatial regularization. Indeed, on each map, the pixels with highest UF are on the border of the image and with only one neighbor.

The inferred maps from Wu et al. (2018) are spatially very similar to the MMSE. However, the $G_0$ values differ. In the bright Car I-E region, the article found $G_0 \simeq 2 \times 10^4$, while in the MMSE $G_0 \simeq 5 \times 10^3$. Figure 6.10 shows that the point from the article falls in a low posterior density region, even though it is close to the posterior mode. This difference is due to the Meudon PDR code version change – Wu et al. (2018) relied on version 1.5.4 with polycyclic aromatic hydrocarbons (PAH) while we resort to version 1.7 with its standard dust population, which does not include PAHs. To check whether the differences came from the difference in the model used and not from the inversion procedure, we trained a neural network emulator of the grid that was used in the article and performed the inversion again. The $G_0$ values in the obtained MMSE maps were similar to those of Wu et al. (2018).
<figure><img src="image_58.png" /><figcaption>Figure 6.8 – Inference results for the Carina nebula.</figcaption></figure>
<figure><img src="image_59.png" /><figcaption>Figure 6.9 – Inference results for the Car I-S pixel of the Carina nebula. Two-dimensional marginal histograms in the physical parameters $\Theta$ space. All histograms are in logarithmic norm.</figcaption></figure>
<figure><img src="image_60.png" /><figcaption>Figure 6.10 – Inference results for the Car I-E pixel of the Carina nebula. Two-dimensional marginal histograms in the physical parameters $\Theta$ space. All histograms are in logarithmic norm.</figcaption></figure>
## 6.4 Orion molecular cloud 1 (OMC-1)

The first two applications considered observations that had already been previously studied with the Meudon PDR code. In this section, we apply our inference procedure on the Orion molecular cloud 1 (OMC-1) far infrared and millimetric maps described in Goicoechea et al. (2019). This inference is the first performed on the full OMC-1 map. Unlike in the two previous applications, it is not clear a priori whether the considered lines can constrain the physical parameters.

The OMC-1 is a bright region of the Orion A cloud. It is located at about 414 pc (Menten et al., 2007) from the Sun, making it the closest region of star formation for stars of intermediate and high masses. Figure 6.11 shows its position in Orion A and its general structure. OMC-1 is illuminated by the Trapezium star cluster, which contains many heavy stars. The cluster lies within the HII bubble. It heats and photodissociates the Orion bar and the East PDR. The two major star-forming sites are the Becklin-Neugebauer/Kleinmann-Low region (BN/KL) and Orion South (Orion S). The Orion bar is one of the most famous and studied PDRs.

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

(a) The Orion A star-formation cloud seen by ESA's Herschel space observatory. The black rectangle outlines the observed region of OMC-1, centered on the Trapezium star cluster. The image is a composite of the wavelengths of 70 µm (blue), 160 µm (green) and 250 µm (red), i.e., dust emission. It spans about 1.3 × 1.6 deg². The moon is shown for scale. Credits: ESA/Herschel/Ph. André, D. Polychroni, A. Roy, V. Könyves, N. Schneider for the Gould Belt survey Key Programme.

(b) OMC-1 composite image covering ~85 arcmin² at ~12" with three emission lines: the C⁺ 158 µm line which traces the FUV-illuminated surface of the molecular cloud (red), the C¹⁸O, J = 2 → 1 line which traces cold dense gas (blue) and the HCO⁺, J = 3 → 2 line (green). Taken from Goicoechea et al. (2019).

<caption>Figure 6.11 – Structure of the OMC-1 cloud.</caption>

**Inversion setup** – The considered observation Y contains L = 4 molecular emission lines: CH⁺ (J = 1 → 0), CO (J = 2 → 1), CO (J = 10 → 9), and HCO⁺ (J = 3 → 2). Other lines such as C⁺ 158 µm and ¹³CO J = 2 → 1 were measured. We tried performing inversion with
these lines, which resulted in an incompatibility between the models and observations. Therefore,
these lines are removed for this analysis.

Figure 6.12 shows the observation maps associated with the *L* = 4 lines used in the inversion. The BN/KL region is masked (in white, in the center of the maps) as it is known to be dominated by shocks. As such events are not included in PDRs models, this region is removed for inversion to avoid biasing the neighboring pixels with the spatial regularization. Dedicated models such as the Paris-Durham code (Godard et al., 2019) would better simulate this region. Note that the model assessment approach might have rejected the model for these pixels. Four pixels are highlighted in this figure: in the Orion bar (lowest square), in the East PDR (mid-height square), in the North-East border of the map (top left square) and in the North-West ridge edge (top right square). Results on these pixels will be detailed.

<figure><img src="image_62.png" /><figcaption>Figure 6.12 – Observations of OMC-1 used for inversion. The white region is the middle of the maps is the BN/KL region, dominated by shocks, and not considered in the inversion. The remaining white pixels correspond to negative integrated intensities due to the additive Gaussian noise. The inference results will be detailed for the four pixels highlighted with a black square.</figcaption></figure>

The observation model is identical to Eq. 6.2. The integrated intensity maps **Y** and the associated additive noise standard deviations ($\sigma_{a,n\ell}$)$_{n\ell}$ are computed directly from the hyperspectral cubes. The multiplicative noise is modeled as a lognormal distribution that combines two terms. The model misspecification is set so that a 3σ error corresponds to an error of a factor 3, i.e., $\sigma_{\text{mod}} = \frac{1}{3}\ln 3 \approx \ln 1.44$. An estimated 8% calibration error is considered for the whole map, i.e., $\sigma_c = \ln 1.08$. The total multiplicative standard deviation is $\sigma_m = \sqrt{\sigma_{\text{mod}}^2 + \sigma_c^2} \approx 0.373 \approx$
ln 1.452. The parameters $a_\ell$ of the likelihood approximation should be optimized for each pixel, which is very computationally demanding for $N > 10^3$ pixels. For simplicity, we set $a_\ell$ to high values so that only the additive Gaussian approximation is used. After testing different values, the spatial regularization parameter $\tau$ is set to 1 for the four physical parameters to obtain smooth and physically consistent maps.

The observation angle $\varphi$ between the surface of the cloud and the line of sight varies in the observed map. For instance, the Orion South (Orion S) region is considered to be mostly face-on ($\varphi \sim 0$ deg) while the Orion bar is close to edge-on ($\varphi \approx 90$ deg). To avoid inferring more parameters than observed lines, we set $\varphi = 0$ deg. Like in the Carina nebula inversion, we assume that the inclination effect will be approximately captured by the scaling factor $\kappa$.

Like for the Carina nebula, we compare the MLE, as a naive estimator yet widespread in astrophysics – see Chapter 3 –, and the MMSE from the method proposed in this thesis. The MLE does not exploit spatial regularization and handles all pixels independently. It is obtained with the optimization algorithm used for the Carina nebula run for 500 iterations. At each step, the optimization variant of the MTM kernel is used with probability $p_{\text{MTM}} = 0.02$. The proposal distribution is set to the smooth uniform prior, and the number of candidates to $K = 50$. The preconditioned gradient descent update step size is set to $\eta = 0.03$.

The MMSE is evaluated from a Markov chain of $T_{\text{MC}} = 10\,000$ iterates, including $T_{\text{BI}} = 500$ of burn-in, generated with the proposed MCMC algorithm. At each iteration $t$, the MTM kernel is selected with probability $p_{\text{MTM}} = 0.5$. It generates $K = 50$ candidates for each pixel $\theta_n$. As in the Carina inversion, the proposal $q$ is set to the Gaussian mixture defined with the neighboring pixels described in Chapter 5 (Eq. 5.35). The PMALA step size is set to $\eta = 0.1$.

**Model assessment and Bayesian p-value** – Figure 6.13 shows the map of estimated p-values $(\hat{p}_n^{(T_{\text{MC}})})_{n=1}^N$. As computing the effective sample size (ESS) is costly for large observation maps, the uncertainties due to the Monte Carlo estimation of the p-values is not accounted for. However, Figure 6.13a shows that the estimated p-values is larger than 0.5 on the regions of interest of the maps. As a result, as shown in Figure 6.13b, only 7 pixels out of 2475 lead to a model rejection, i.e., less than 0.3%. In particular, none of the 4 highlighted pixels lead to a model rejection.

Figure 6.14 shows that the predicted integrated intensities are very close to the actual observations. The only exception is the CH$^+$ $J = 1 \rightarrow 0$ emission line for the pixel in the North-West ridge edge, where the MMSE is a factor $\sim 3-5$ too high. Note that the two pixels from PDRs contain brighter lines than the two other pixels, which is also visible in Figure 6.12.

<figure><img src="image_63.png" /><figcaption>Figure 6.13 – Model assessment for OMC-1.</figcaption></figure>
<figure><img src="image_64.png" /><figcaption>Figure 6.14 – Posterior predictive assessment for four pixels of OMC-1 on the lines used for the inversion. Comparison of observations **Y** and associated noise model with posterior predictive distributions on $\hat{f}(\Theta)$, with $\Theta \sim \pi(\Theta|\mathbf{Y})$.</figcaption></figure>
**Inference results: MLE, MMSE, and UF maps** – Figure 6.15 shows inference results for both likelihood minimization and posterior sampling. The spatial structures in the observations are recovered with the MMSE when including the spatial smoothness prior, but hardly visible when using a simple MLE. Due to the low number of lines *L*, the spatial regularization plays a key role in recovering the correct spatial structure. As there are few lines and thus few constraints for the models, the likelihood is very sensitive to noise and the MLE presents unphysical variations. Conversely, the spatial regularization permits each pixel to access the information contained in its neighbors, which share similar physical conditions.

We now focus on the MMSE and UF maps. The MMSE permits to recover the overall structure of OMC-1. The Orion bar, the HII region, the Orion S region, the East PDR and the North-West ridge are visible on each of the *D* = 4 reconstructed maps. The scaling factor *κ* is close to 1 in the bright regions. In the Orion S region, in the East PDR and in the Orion bar, the thermal pressure *P*<sub>th</sub> is estimated at ≈ 4 × 10<sup>8</sup> K cm<sup>-3</sup>, and the radiative intensity factor *G*<sub>0</sub> at ≈ 5 × 10<sup>4</sup>. In these regions, the recovered visual extinction *A*<sub>V</sub><sup>tot</sup> is roughly 10 mag. These values are consistent with the literature and with the analysis of the one-pixel observation of the Orion bar presented in Appendix 6.A. Our results show for the first time that the high thermal pressures previously found in the Orion Bar appear to be widespread in the PDRs of OMC-1.

Over the whole map, the UF is lower than 1.5 for the scaling factor $\kappa$, than 2 for the visual extinction $A_V^{\text{tot}}$. This means that the typical error is at most 50% for $\kappa$, and of a factor 2 for $A_V^{\text{tot}}$. Despite using only $L=4$ lines to infer $D=4$ parameters for each pixel, these two parameters are thus well constrained. Similarly, on most of the map, the UF is lower than 3 for the thermal pressure $P_{\text{th}}$ and than 4 for the intensity of the incident radiative field $G_0$. These two parameters are thus relatively well constrained considering their large validity intervals and the low number of lines used for inversion ($L=4$). The UF maps of $P_{\text{th}}$ and $G_0$ show that the darker regions yield higher associated uncertainties. Besides, the UF attains a maximum in the North-West ridge edge, around one of the highlighted pixels.

**Inference results: analysis of highlighted pixels** – Figure 6.16 shows the pairwise histograms for the Orion bar pixel (lowest of the four squares). Figure 6.17 shows the pairwise histograms for the East PDR pixel (mid-height square). They show that the marginal posterior distributions for these two pixels are simple, as they could be well approximated by a Gaussian distribution with a diagonal covariance matrix. Uncertainties are large but do not show degeneracy between parameters. Figure 6.18 shows the North-East pixel (top left square). This pixel and its surrounding region yields a more moderate $P_{\text{th}}$ (≈ 7 × 10<sup>7</sup> K cm<sup>-3</sup>) while the $G_0$ remains similar to what is found in the Orion Bar (≈ 5 × 10<sup>4</sup>). This seems to contradict the correlation between $P_{\text{th}}$ and $G_0$ described in Joblin et al. (2018) and Wu et al. (2018). Further investigation is thus required to better understand this specific region.

Finally, Figure 6.19 shows that the UF peak around the North-West ridge edge pixel (top right pixel) corresponds to a bimodality in the joint distribution ($P_{\text{th}}, G_0$). The sampler successfully identified the two modes in the ($P_{\text{th}}, G_0$) joint distributions, which is challenging for standard MCMC algorithms. The two modes correspond to two different environments. One has a very high $G_0$ ≈ 3 × 10<sup>4</sup> and a moderate $P_{\text{th}}$ ≈ 3 × 10<sup>7</sup> K cm<sup>-3</sup>. The other, which is more likely, has a moderate $G_0$ ≈ 5 × 10<sup>3</sup> and a very high $P_{\text{th}}$ ≈ 3 × 10<sup>8</sup> K cm<sup>-3</sup>. Different hypotheses may explain this bimodality. For instance, the existence of two solutions may be due to the low number of lines $L=4$. In this case, adding other lines for the inversion may remove one mode. Another possibility is that this pixel mixes emissions from distinct environments at different positions along the line of sight. This is also possible as this pixel is on the edge of a ridge. Note that although the MMSE values represented on the reconstructed maps fall in a high probability region, it does not fall in a mode. In such a case, the MAP would be a more representative estimator of maps of $P_{\text{th}}$ and $G_0$, but would require running an additional optimization algorithm.
<figure><img src="image_65.png" /><figcaption>Figure 6.15 – Inference results for OMC-1.</figcaption></figure>
<figure><img src="image_66.png" /><figcaption>Figure 6.16 – Inference results for OMC-1: Orion bar pixel. Two-dimensional marginal histograms in the physical parameters $\Theta$ space. All histograms are in logarithmic norm.</figcaption></figure>
<figure><img src="image_67.png" /><figcaption>Figure 6.17 – Inference results for OMC-1: East PDR pixel. Two-dimensional marginal histograms in the physical parameters $\Theta$ space. All histograms are in logarithmic norm.</figcaption></figure>
<figure><img src="image_68.png" /><figcaption>Figure 6.18 – Inference results for OMC-1: North-East edge pixel. Two-dimensional marginal histograms in the physical parameters $\Theta$ space. All histograms are in logarithmic norm.</figcaption></figure>
<figure><img src="image_69.png" /><figcaption>Figure 6.19 – Inference results for OMC-1: North-West ridge edge pixel. Two-dimensional marginal histograms in the physical parameters $\Theta$ space. All histograms are in logarithmic norm.</figcaption></figure>
## 6.5 Conclusion

In this chapter, we applied the full inference procedure to three real observations – NGC 7023, the Carina nebula and OMC-1. The preliminary analysis of the presented inference results already provides interesting astrophysical insights. For instance, Joblin et al. (2018) and Wu et al. (2018) describe a positive correlation between point estimates of $P_{\text{th}}$ and $G_0$ among multiple sources. In contrast, our inference results on NGC 7023, the Carina nebula and the Orion bar revealed a negative correlation between these two parameters in their joint distribution. Putting things together, this result indicates that the positive correlation among sources from the literature was not due to inference uncertainties, and is likely to have a physical origin. However, in view of the large estimation uncertainties, deriving a numerical relation between $P_{\text{th}}$ and $G_0$ from point estimates seems inaccurate. In particular, the power law between proposed in Joblin et al. (2018) and Wu et al. (2018) may be improved by accounting for these inference uncertainties.

We emphasize that this is the first study of maps of ionic, atomic or molecular line observations that permits to infer maps of physical parameters and to quantify the associated uncertainties at once. In particular, the proposed method proved to be able to detect multiple modes and to compare their relative weights in the posterior distribution. It also proved to be able to detect correlations between inferred physical parameters. These properties lead to finer and more trustful results and interpretations.
# Appendix 6.A Orion Bar

The Orion Bar PDR lies ~ 2' south-east of the Trapezium star cluster (Allers et al., 2005). This cluster of massive stars creates a HII region that penetrates into the parent molecular cloud. The distance between the Orion nebula and the Sun has been measured at 414 ± 7 pc (Menten et al., 2007). Because of its proximity and edge-on geometry, the Bar is one of the most studied PDRs.

Joblin et al. (2018) analyzed the Orion bar with an observation of $L = 24$ lines in $N = 1$ pixel. The observed lines used in the inversion include $^{12}$CO lines (from $J = 11 \rightarrow 10$ to $J = 23 \rightarrow 22$), rotational H$_2$ lines (from S(0) to S(5)) and low level CH$^+$ rotational lines (from $J = 1 \rightarrow 0$ to $J = 6 \rightarrow 5$). The noise on the observation y was assumed additive, Gaussian and uncorrelated with known standard deviations ($\sigma_{a,l}$)$_{l=1}^L$. Like NGC 7023, the authors inferred $\theta = (\kappa, P_{\text{th}})$, while fixing the observation angle to $\phi = 60$ deg and the total extinction to $A_V^{\text{tot}} = 10$ mag – as the observed lines did not provide sufficient constraint on $A_V^{\text{tot}}$. Based on Tielens and Hollenbach (1985) and Marconi et al. (1998), they set $G_0 = 2 \times 10^4$. Like NGC 7023, the fit was performed with a grid search on $P_{\text{th}}$ and a simple continuous optimization on $\kappa$.

They obtained $P_{\text{th}} = 2.8 \times 10^8 \text{ K cm}^{-3}$ and $\kappa = 1.3$.

**Inversion setup** – In this analysis, we use the same observation model as for NGC 7023 (Eq. 6.2). The prior hyperparameters and validity intervals are identical. The parameters of the proposed sampler are the same as in the previous section, except for the step size of the PMALA kernel. We empirically set it to $\eta = 0.03$ to improve the acceptance probability.

**Model assessment and Bayesian p-value** – In Joblin et al. (2018), the authors stated that the fit did not reproduce the observations as well as for NGC 7023. Applying our inference procedure, we obtain an estimated p-value is $\hat{p}^{(T_{\text{MC}})} = 0.52$, and a rejection probability of $\mathbb{P}[p^{(t)} \le \alpha] < 10^{-200} < \delta$, which does not lead to a rejection. Therefore, the considered observation model is compatible with the observations. Figure 6.20 compares the marginal posterior predictive distributions with the true observations. Most H$_2$ lines and the CH$^+$ $J = 1-0$ line are relatively poorly reconstructed, with an error of a factor ~3 for the MMSE on these lines. However, this factor 3 remains compatible with the errors $\sigma_{a,l}$ and $\sigma_m$ in the observation model.

**Inference results** – Figure 6.21 shows our estimation results. The anticorrelation between $P_{\text{th}}$ and $G_0$ observed in NGC 7023 is also visible here. This is expected as we are using almost the same set of lines as observational constraints. As in NGC 7023, the results for $P_{\text{th}}$, $G_0$ and $A_V^{\text{tot}}$ are overall consistent with the values from Joblin et al. (2018). Table 6.2 details the obtained values. The obtained $G_0 = 3.0 \times 10^4$ – with $G_0 \in [4.8 \times 10^3, 1.3 \times 10^5]$ with probability 95% – is consistent with the literature. The inferred thermal pressure $P_{\text{th}} = 1.1 \times 10^8 \text{ K cm}^{-3}$ – with $P_{\text{th}} \in [4.7 \times 10^7, 2.8 \times 10^8] \text{ K cm}^{-3}$ with probability 95% – is slightly lower than the $2.8 \times 10^8 \text{ K cm}^{-3}$ from Joblin et al. (2018). Conversely, the scaling factor $\kappa = 5.2$ – with $\kappa \in [2.7, 10]$ with probability 95% – is much larger than 1.3. Like in NGC 7023, this might be due to the change of Meudon PDR version – 1.5.4 in Joblin et al. (2018), and 1.7 in our work. It may also be due to the difference in some secondary parameters.
<figure><img src="image_70.png" /><figcaption>Figure 6.20 – Posterior predictive assessment for the Orion Bar on the lines used for the inversion. Comparison of observations **Y** and associated noise model with posterior predictive distributions on $\hat{f}(\Theta)$, with $\Theta \sim \pi(\Theta|\mathbf{Y})$.</figcaption></figure>

<table><caption>Table 6.2 – Inference results on Orion Bar. In Joblin et al. (2018), only the scaling parameter $\kappa$ and the thermal pressure $P_{\text{th}}$ are inferred.</caption><thead><tr><th></th><th></th><th>κ<br>-</th><th>P<sub>th</sub><br>(K cm<sup>-3</sup>)</th><th>G<sub>0</sub><br>-</th><th>A<sub>V</sub><sup>tot</sup><br>(mag)</th></tr></thead><tbody><tr><td rowspan="2">(Joblin et al., 2018)</td><td rowspan="2">MMSE</td><td>1.3</td><td>2.8 × 10<sup>8</sup></td><td>2 × 10<sup>4</sup></td><td>10</td></tr><tr><td>5.2</td><td>1.1 × 10<sup>8</sup></td><td>3.0 × 10<sup>4</sup></td><td>11.6</td></tr><tr><td rowspan="3">68% credibility interval</td><td>lower bound l<sub>68%</sub></td><td>3.7</td><td>8.0 × 10<sup>7</sup></td><td>9.5 × 10<sup>3</sup></td><td>4.1</td></tr><tr><td>upper bound u<sub>68%</sub></td><td>7.5</td><td>1.6 × 10<sup>8</sup></td><td>8.5 × 10<sup>4</sup></td><td>31.5</td></tr><tr><td>UF<sub>68%</sub></td><td>1.4</td><td>1.4</td><td>3.0</td><td>2.8</td></tr><tr><td rowspan="3">95% credibility interval</td><td>lower bound l<sub>95%</sub></td><td>2.7</td><td>4.7 × 10<sup>7</sup></td><td>4.8 × 10<sup>3</sup></td><td>1.9</td></tr><tr><td>upper bound u<sub>95%</sub></td><td>10.0</td><td>2.8 × 10<sup>8</sup></td><td>1.3 × 10<sup>5</sup></td><td>42.4</td></tr><tr><td>UF<sub>95%</sub></td><td>1.9</td><td>2.5</td><td>5.2</td><td>4.7</td></tr></tbody></table>
<figure><img src="image_71.png" /><figcaption>Figure 6.21 – Inference results for the Orion Bar. Two-dimensional marginal histograms in the physical parameters $\Theta$ space. All histograms are in logarithmic norm.</figcaption></figure>
# Conclusions and perspectives

"People from different backgrounds approach a subject in different ways and ask different questions."

Jocelyn Bell Burnell

## Contents

<table><tr><td>Conclusions</td><td>171</td></tr><tr><td>Contributions in statistics</td><td>172</td></tr><tr><td>Contributions in astrophysics</td><td>173</td></tr><tr><td>Perspectives</td><td>174</td></tr></table>

## Conclusions

The study of the ISM carries fundamental questions such as the regulation of star formation or the development of molecular complexity, possibly leading to the formation of prebiotic molecules. In this thesis, we studied the impact of the radiative feedback of newborn massive stars on their parent molecular clouds. New facilities at IRAM, ALMA and the JWST might lead to breakthroughs in the coming years, thanks to the very rich hyperspectral data they provide. For instance, the IRAM-30m Large Program "Orion B" observed the Orion B cloud at dense core resolution, resulting in a million-pixel map, with 240 000 spectral bands containing emissions of dozens of ionic, atomic or molecular tracers (Pety et al., 2017). In this thesis, we addressed five problems to extract as much information as possible from such observations combined with an ISM numerical model such as the Meudon PDR code:

- Deriving fast, accurate and light emulators of the Meudon PDR code, an ISM numerical model, to be able to make inference scalable. This task was addressed in Chapter 4.

- Approximating the likelihood function to obtain a simple closed-form expression without neglecting a source of noise. The obtained approximation simplifies the sampling task while introducing a small and controlled error. This task was addressed in Chapter 5 (Section 5.1.1).

- Proposing a spatial regularization prior to improve the quality of estimations in low signal-to-noise ratio regions. This task was addressed in Chapter 5 (Section 5.1).

- Estimating physical parameter maps such as the thermal pressure and the total visual extinction from observations, to bring new insights on the interstellar medium (ISM) and star formation. The proposed Bayesian inference method accounts for as much physics as possible – in the forward model, in the observation model, and in the prior distribution – and scales well enough to be applied to observation maps of up to $O(10^4)$ pixels. It
provides uncertainty quantification through credibility intervals. This task was addressed in Chapter 5 (Section 5.2).

— Assessing the compatibility between the observation model and the observations, to potentially provide feedback to the astrophysicists who build these models. This task was addressed in Chapter 5 (Section 5.3).

An additional task was also addressed. It consists in selecting the most informative observables to minimize uncertainty on inferred physical parameters. The resulting variable selection method is based on the conditional differential entropy – or, equivalently, on the mutual information. Although we already have preliminary results, this work is still in progress in collaboration with other members of the ORION-B consortium. This task was not presented in the body of this manuscript. The interested reader can find a description of the preliminary results on our side, i.e., on the numerical model side, in Appendix A.

As part of an interdisciplinary project, this thesis yielded contributions in both statistics and astrophysics for each of the considered problems.

## Contributions in statistics

The main methodological contributions of this thesis are the proposition of a general method to build emulators, of a likelihood approximation, of a new Markov chain Monte Carlo (MCMC) algorithm, and of an extension of the Bayesian test of hypothesis.

**Derivation of fast and accurate approximations of ISM numerical models** – State-of-the-art ISM numerical models have prohibitively long evaluation times. They are thus often replaced with approximations for inference, which induces an error. Interpolation methods are often used in the ISM community (Wu et al., 2018; Ramambason et al., 2022), and the associated error is not always quantified. In Chapter 4, we proposed to apply general strategies to emulate a numerical model that predicts many observables at once and punctually produces outliers. We applied these strategies to emulate the Meudon PDR code, and compared the obtained artificial neural networks (ANNs) with interpolation methods with respect to memory requirements, speed, and accuracy. The proposed ANNs significantly outperformed all the considered interpolation methods.

**Likelihood approximation** – The observation model introduced in Chapter 3 (Section 3.4), which involves a multiplicative lognormal noise and an additive Gaussian noise, does not lead to a simple closed-form likelihood function. In Section 5.1.1, instead of neglecting one source of noise, we proposed a parametric closed-form approximation of the likelihood function. This approximation builds on that from Nicholson and Kaipio (2020), which is purely additive and Gaussian. We combined a similar additive and Gaussian approximation with a multiplicative and lognormal one using a weighted geometric mean. The weight of this geometric mean is a function of the physical parameters such that the multiplicative approximation dominates in the high SNR regime, and the additive approximation in the low SNR regime. In Section 5.B, we proposed a grid search and a Bayesian optimization approach to tune the weight parameters so that the approximation error is minimized.

**A new MCMC algorithm** – The posterior distribution presented in Chapter 5 (Section 5.1.3) is challenging to sample from. The log-posterior is considered to be non log-concave and non gradient Lipschitz continuous, or with an unusable Lipschitz constant. We proposed in Section 5.2 a new MCMC algorithm. This algorithm combines two sampling kernels, each addressing one of these two difficulties.

The first kernel is based on a preconditioned MALA algorithm (Xifara et al., 2014) that exploits information on the local geometry of the log-posterior to propose relevant candidates. We resorted to the RMSProp preconditioner (Tieleman and Hinton, 2012), which was already
applied in an MCMC algorithm (Li et al., 2016). We improved the algorithm from Li et al. (2016) by performing accept-reject steps and by accounting for the correction term $\gamma$ that emerges from a position-dependent preconditioner.

The second kernel combines three algorithms. A Gibbs sampling (Geman and Geman, 1984) divides the physical parameter of dimension $N \times D$ into $N$ conditional distributions of dimension $D$ that are easier to explore. A multiple-try Metropolis algorithm (Martino, 2018) increases the probability of acceptance, and uses a proposal distribution that allows escapes from local modes. A chromatic Gibbs sampling algorithm (Gonzalez et al., 2011) permits to perform the sampling of many pixels in parallel, which greatly accelerates the sampling.

**Extension of Bayesian hypothesis testing for model assessment** – Solving an inverse problem relies on the hypothesis of a possibly misspecified observation model. The Bayesian test described in Gelman et al. (1996) permits to assess the model compatibility with the observations. It relies on a test statistic that permits to compare true observations with observations reproduced from the observation model. For simple cases such as a point estimate, a Gaussian additive uncorrelated noise and the $L_2$ norm as the test statistic, the $p$-value associated with the test can be computed exactly. Otherwise, it is approximated by a Monte Carlo (MC) estimator, which induces an error that can lead to a wrong decision on the compatibility or incompatibility of the model with the observations. In Chapter 5 (Section 5.3), we extended the test to account for the uncertainty due to the MC approximation error. We considered the $p$-value as a random variable. Using a simple likelihood model and a conjugate prior, we obtained a simple uncertainty description of the $p$-value. This description represents a negligible additional computational cost with respect to the standard MC estimator. The proposed test can detect cases where the uncertainty on the MC estimator is too large to reject or not the numerical model. Too large uncertainties indicate that longer Markov chains should be used.

## Contributions in astrophysics

Chapter 3 reviewed applications of statistical inference in ISM studies, and the position of our work with respect to the literature. We now review our contributions for the ISM community, namely the derivation of a fast, accurate and light surrogate for numerical models, the use of a spatial regularization prior, the use of Bayesian inference for high dimensional observation maps of integrated intensities, the application of the model assessment test, and the insights obtained from the application of our methods to real data.

**Informative spatial regularization prior** – To exploit the map structure of the physical parameters, we resorted to a spatial regularizing prior – introduced in Chapter 5 (Section 5.1.2). Like Marchal et al. (2019), we use a $L_2$ norm of the map Laplacian. However, this article aimed at fitting one mixture of Gaussian profiles per spectrum on a hyperspectral map, while we directly infer more complex physical parameters. A common approach in the ISM community for low SNR regions is to stack pixels to increase the SNR, thus reducing spatial resolution. This spatial regularization approach permits to better exploit low SNR observations, as it enables pixels to access the information contained in their neighbors. In addition, the spatial regularization is negligible when the likelihood function constrains well the physical parameters, i.e., in high SNR observations and good tracers of the physical parameters. It can thus be seen as a form of adaptive stacking.

**Bayesian inference for high dimensional maps** – When the observed lines constrain poorly the physical conditions or when the Signal-to-Noise Ratio is low, multiple solutions might reconstruct observations equivalently well. Moreover, the non-linearity of astrochemistry models leads to non-convex and even multi-modal problems. Many methods used in interstellar astrophysics do not consider these degeneracies and only return one estimated map (Joblin et al., 2018; Wu et al., 2018). Most methods get trapped in local minima and estimations have no optimality guarantees.
Heuristics can overcome this issue for small maps but would be unrealistically slow for larger ones. In Chapter 6, using a Bayesian approach, we obtained meaningful uncertainty quantification on inferred physical parameters. In the ISM community, dust studies such as Galliano (2018); Galliano et al. (2021) had already applied Bayesian sampling methods in high dimensional settings. However, sampling methods had never been used – to the best of our knowledge – to reconstruct maps of physical parameters from observations of integrated intensities of ionic, atomic or molecular emission lines.

**Model assessment with Bayesian hypothesis test** – In ISM studies, the compatibility of the model with the observations is rarely assessed. However, recently, Lebouteiller and Ramambason (2022) used multiple information criteria, the Bayesian evidence and a model assessment to estimate the number of distinct environments that are necessary to explain observations. Similarly, Galliano (2022) used the Bayesian hypothesis test to assess the quality of the model. In Chapter 6, we applied our extended test of hypothesis which includes the error that comes from the MC approximation. We also apply it to each pixel of our reconstructed maps, which permits to identify pixels corresponding to environments that are poorly modeled by the Meudon PDR code. Such pixels may not correspond to PDR, or may pinpoint an insufficiency from the Meudon PDR code. This feedback is therefore valuable for astrophysicists working on numerical models, as it can indicate potentially necessary future development.

**New insights** – Applications of the proposed method to real observations already brought new astrophysical insights. Here, we highlight two of them. First, in Joblin et al. (2018) and Wu et al. (2018), the authors outline a positive correlation between point estimates of the thermal pressure $P_{\text{th}}$ and the intensity of the incident UV radiative field $G_0$ among multiple sources. In contrast, our inference results on NGC 7023, the Carina nebula and the Orion bar revealed a negative correlation between these two parameters in their joint distribution. This result indicates that the positive correlation among sources from the literature is not due to inference uncertainties, and is likely to have a physical origin. However, the large estimation uncertainties indicate that the numerical power law relations between $P_{\text{th}}$ and $G_0$ proposed in Joblin et al. (2018) and Wu et al. (2018) may be improved by accounting for these inference uncertainties.

Second, the presented analysis on the OMC-1 cloud is the first inference performed on the full observation map. The results are consistent with the literature, with e.g., $P_{\text{th}} \simeq 3-5 \times 10^8 \text{ Kcm}^{-3}$ and $G_0 \simeq 3-5 \times 10^4$ in the Orion bar. New information about the region are provided, such as the fact that the pressure is high in all the PDRs around the Trapezium star cluster. Additional analyses are still required for a full physical interpretation of the inference results.

## Perspectives

This manuscript represents a milestone in our research on the ISM and star formation with the ORION-B consortium. This work could now be extended in multiple directions. In the following, we list some promising directions.

**Derivation of better emulators for ISM numerical models** – The proposed neural network-based emulation strategy is quite general and could be applied to other ISM numerical models. In general, numerical models are based on a system of coupled equations, potentially partial differential equations. Deriving an emulator from a set of input-output pairs seems suboptimal as the emulator could in principle be trained to satisfy the system of equations. The physics-informed neural networks (Raissi et al., 2019) form a potentially promising approach to derive surrogate models for the ISM community. Such ANNs are trained directly with the system of physical differential equations. Therefore, they have access to more information during training, and thus generalize better with limited datasets. This seems particularly promising for chemistry models, for which the differential equations are easy to write. This lead was considered for the Meudon
PDR code. However, like many complex ISM codes, it includes many microphysical processes and thus many equations. The physics-informed ANNs were thus badly suited to this task.

A second potentially promising lead regarding the emulator is the management of outliers. We proposed a three-step approach that includes a semi-manual check. Preventing a numerical model from producing such outliers would be statistically simpler, and would improve the confidence in the model predictions. This is very challenging for a code as complex as the Meudon PDR code, in which bistability – the existence of two solutions for an equation system – can occur at different levels. Then, deriving a deterministic binary test that detects outliers would prevent from using any in the training of the emulator. Again, this is very challenging for the Meudon PDR code, as the number of profiles to check is high. Finally, if a statistical approach needs be used for outliers, a potentially better approach would be to learn the surrogate model and fit the binary mask on outliers at once, as detailed in Chapter 4 (Appendix 4.A). However, this approach would require careful hyperparameter tuning to avoid the classification of physically meaningful points that are difficult to reproduce as outliers.

**Observation model** – Several improvements of the observation model could be explored. First, using a hierarchical model instead of the proposed approximation of the likelihood function would permit to use the exact observation model in the inversion. We did try to adapt the proposed sampler to the hierarchical model, but using the approximation proved to considerably simplify the sampling. A possibility we did not explore yet is to resort to the ancillarity-sufficiency interweaving strategy (ASIS) sampling approach, already applied in dust studies, for instance in (Kelly, 2011; Kelly et al., 2012; Galliano, 2018).

Second, the observation model specified in Chapter 3 (Section 3.4) neglected all the correlations in the noise. This choice greatly simplified the definition of the likelihood function and the computation of its first and second order derivative, as it enabled to consider the observation elements $y_{nl}$ individually. However, it is known that the noise realizations in a telescope are not independent, neither spatially nor spectrally. For instance, in Einig et al. (2023), the authors describe the spatial and spectral noise structure in the Orion-B observations. Similarly, as we showed in Chapter 3 (Section 3.1.2), dust studies now frequently include non diagonal terms in the noise covariance. Accounting for these noise correlations would provide additional information in the observation model, which is likely to result in lower uncertainties on estimates.

Finally, the noise model presented in Section 3.4 and exploited for OMC-1 in Section 6.4 slightly overestimates the additive noise standard deviation. Indeed, we set a constant interval for the spectral integration of the intensities, and used large intervals to avoid cutting any signal. Using a line profile fitting algorithm such as ROHSA (Marchal et al., 2019) or CUBEFIT Paumard et al. (2022) may lead to more accurate integrated intensity estimates. Similarly, the denoiser proposed in (Einig et al., 2023) improves the SNR without assuming a line profile. However, these approaches may introduce a bias in the estimated integrated intensities. Besides, it is not simple how to describe the associated uncertainties, which is crucial for the inversion.

To improve the observation model with respect to the last two points, one could infer at once both 1) the multiline maps of integrated intensities from the hyperspectral maps with the complete noise model and 2) the physical parameters from the multiline maps of integrated intensities. The resulting complex posterior distribution would be written with a hierarchical model and sampled from using a Gibbs sampler. Such an approach might define a heavier inverse problem. However, the effective noise on the multiline maps would be reduced, which would lead to smaller uncertainties in inferred physical parameters.

**Choice of a preconditioner** – We chose the RMSProp preconditioner because it is a diagonal preconditioner that proved efficient in high dimensions both in theory and in practice (Dauphin et al., 2015). As it is diagonal and only requires first order derivatives, the correction term $\gamma$ in PMALA that comes from the position-dependent preconditioner can be written and evaluated with acceptable costs – see Chapter 5 (Section 5.2.1). This correction term ensures a good dis-
cretization of a Langevin diffusion process, which favors good candidates. However, accounting
for this term is not necessary to ensure ergodicity or that the posterior is the invariant distribution
of the generated Markov chain, even if the correction term is neglected. Other preconditioners
that take into account the spatial regularization may lead to better mixing properties. As the
spatial regularization prior is equivalent to a Gaussian distribution with a band precision matrix,
implementing PMALA with such a preconditioner could exploit classic methods to sample from
a Gaussian distribution (Vono et al., 2021). Such a preconditioner may perform better than RM-
SProp, e.g., for large spatial regularization parameter $\tau$.

**Analysis of larger observation maps** – The proposed sampler was applied to cases with dimensions up to $ND = 4 \times 10^4$. Although this is already large, new and future maps are even larger. For instance, the Orion-B map (Pety et al., 2017) contains $10^6$ pixels. There are two necessary and independent approaches to scale to these maps.

First, the Orion-B map contains multiple environments. Indeed, the Horsehead nebula is better modeled with a PDR model, but the dense cores in the cloud require a dark cloud time-dependent model. Similarly, the original observations of OMC-1 studied in Chapter 6 (Section 6.4) contain the BN/KL region, which is dominated by shocks. To be able to perform inference on these observations using the Meudon PDR code, some pixels were masked. In addition, observations cannot be spatially well resolved, which means that some pixels may contain multiple environments that require different models. Inference on such large maps would therefore require simultaneously unmixing the environments on each pixel and estimating the physical parameters of each environment. This is part of the PhD project of Antoine Zakardjian started in 2022, member of the ORION-B consortium, supervised by Jérôme Pety and Annie Hughes.

Second, even neglecting these mixtures of environments, the proposed algorithm is unlikely
to scale to these dimensions without additional effort. Distributed algorithms are a common
solution to divide one computationally very expensive task into smaller tasks addressed in parallel
by different workers. Thanks to the simple regularization prior and diagonal noise correlation,
the proposed MCMC algorithm could certainly be distributed, using for instance the framework
presented in Thouvenin et al. (2023). This is also in part the PhD project of Maxime Bouton
started in 2022, supervised by Pierre Chainais and Pierre-Antoine Thouvenin.
# Appendix A

# Entropy of probability distributions and line selection

"You're going to have to make a choice, Mr. Anderson."

The Matrix (1999), Lana and Lilly
Wachowski

## Contents

<table><tr><td>A.1 A line selection method that compares posterior distributions</td><td>178</td></tr><tr><td>A.1.1 Using the MSE as a quantitative criterion to rank sets of lines</td><td>178</td></tr><tr><td>A.1.2 Differential entropy and the Fano Bound</td><td>179</td></tr><tr><td>A.1.3 Selection of the best set of lines as a discrete optimization problem</td><td>180</td></tr><tr><td>A.2 Illustration: a Gaussian and linear case</td><td>180</td></tr><tr><td>A.3 Conclusion</td><td>182</td></tr><tr><td>Appendix A.A Derivation of the multivariate Fano Bound</td><td>184</td></tr></table>

In Chapter 4, we derived a surrogate model of the full Meudon PDR code. This surrogate model emulates all the 5375 emission lines predicted by the Meudon PDR code. Chapter 5 introduced the posterior distribution involved in the inverse problem considered in this thesis. In actual ISM observations, due to observational constraints, only $L \simeq 5$ to 30 lines are observed. Observations contain different lines depending on the telescope or observed wavelength range. In particular, observing lines in the infrared (IR) domain requires space telescopes, as these wavelengths are absorbed by the Earth atmosphere. In general, observing maps of many lines is extremely expensive. Besides, observing all lines at once is unfeasible due to the large differences in wavelengths, from near infrared to millimetric domains. It is therefore crucial to prepare observations to measure lines that permit to infer physical parameters of interest with low uncertainty.

In this chapter, we present a statistical method that determines, for a fixed number of emission lines, the set of lines that leads to the lowest uncertainty on inferred physical parameters. This uncertainty is quantified with the conditional differential entropy – or equivalently, the mutual information. We evaluate the best set of lines for a variety of true physical conditions $\theta \in \mathbb{R}^D$ and for different sizes $L$ of sets.

This project is a collaboration with Lucas Einig and Antoine Roueff, also members of the ORION-B consortium. The overall goal is to apply this variable selection strategy both on nu-
merical models and the Orion-B dataset (Pety et al., 2017), and to compare the results. In this appendix, we only present the theoretical developments of a line selection method. This is still ongoing work and an application to the Meudon PDR code with the noise model presented in Chapter 3 is currently under study.

Section A.1.2 introduces the mean squared error (MSE) considered in this context, the conditional differential entropy and the Fano bound that links the two. Section A.2 applies the Fano bound to a simple case with an affine forward model and an additive uncorrelated Gaussian noise to provide intuitive insights.

# A.1 A line selection method that compares posterior distributions

In this chapter, we aim at determining the set of $K \in [1, L]$ lines that leads to the lowest uncertainty on inferred physical parameters, i.e., that is the most informative. This uncertainty is quantified with a general mean squared error (MSE). As this MSE is in general not accessible, we resort to a lower bound on the MSE called the Fano bound. The Fano bound relies on the conditional differential entropy, that can be evaluated with a Monte Carlo (MC) estimator from samples. The set of K lines that minimizes the conditional differential entropy also minimizes the Fano bound.

## A.1.1 Using the MSE as a quantitative criterion to rank sets of lines

**Distribution on the observation $\pi(y|\theta^*)$** – The results of such a search depend on the type of observed environment. We account for this dependence and determine the most informative set of K lines for different physical conditions $\theta^* \in \mathbb{R}^D$. Applying the observation model presented in Chapter 3 (Eq. 3.23) with the “true” physical parameters $\theta^*$, $\omega = -\infty$ – no censorship –, and a set of lines $S \in 2^{[1,L]}$ yields

$$ \forall \ell \in S, \quad y_{\ell} = \varepsilon_{n\ell}^{(m)} \tilde{f}_{\ell}(\theta^*) + \varepsilon_{n\ell}^{(a)}, \quad \varepsilon_{n\ell}^{(a)} \sim \mathcal{N}(0, \sigma_a^2), \quad \varepsilon_{n\ell}^{(m)} \sim \text{Lognormal}\left(-\frac{\sigma_m^2}{2}, \sigma_m^2\right). \quad (A.1) $$

This observation model defines a distribution on observations denoted $\pi(y|\theta^*)$.

**Joint distribution $\pi(\theta, y)$** – From each observation y, the posterior distribution $\pi(\theta|y)$ is defined as in Chapter 5 (Eq. 5.12). In this case, as we consider a single pixel, there is no spatial regularization. The combination of the distribution on observations (Eq. A.1) and the posterior distribution defines a joint distribution on $(\theta, y)$

$$ \pi(\theta, y) = \pi(\theta|y)\pi(y|\theta^*). \quad (A.2) $$

Sampling from this joint distribution is equivalent to solving one inverse problems per sample $y \sim \pi(y|\theta^*)$, i.e., per realization of noise on an observation of the environment characterized by $\theta^*$.

**A general mean squared error (MSE) as quantitative criterion** – To select a set of lines $s \subset [1, L]$, we need to compare joint distributions on $(\theta, y)$ with respect to a quantitative criterion. As already mentioned, the goal of this chapter is to determine the set of lines that best constrains the physical parameters. The MSE quantifies this notion of well constrained parameters. Here, the MSE is integrated with respect to both observation y and physical parameters $\theta$,

$$ \text{MSE} = \mathbb{E}_{(\theta, y)} \left[ \|\theta - \hat{\theta}(y)\|_2^2 \right] = \int \|\theta - \hat{\theta}(y)\|_2^2 \pi(\theta, y) d\theta dy \quad (A.3) $$

The MSE is minimized with one estimator $\hat{\theta}$: $y \mapsto \hat{\theta}(y)$, the posterior expectation $\hat{\theta}(y) = \mathbb{E}[\theta|y]$. Chapters 2 and 5 covered how to evaluate the posterior expectation for a single value of
y using an Monte Carlo (MC) estimator. When y follows a distribution, evaluating the posterior expectation would require running an MCMC algorithm for each y value, which is highly inefficient.

An alternative to using the MSE directly to rank sets of lines is to compute a lower bound on the MSE. In the following, we consider an entropy-based lower bound, called the Fano bound.

### A.1.2 Differential entropy and the Fano Bound

Selecting the set of lines that best constrains the physical parameters is equivalent, in a sense, to selecting the set that brings most information on the physical parameters, i.e., that minimizes the uncertainties in inference. This notion of uncertainty can be quantified using the entropy of a distribution.

**Differential entropy** – The differential entropy is similar to a Shannon entropy counterpart for continuous random variables. For $\theta \in \mathbb{R}^D$,

$$h(\theta) = E_{\theta}[-\ln \pi(\theta)] = -\int \ln \pi(\theta) \pi(\theta) d\theta, \quad (A.4)$$

with $\pi(\theta)$ any distribution on $\theta$. Table A.1 provides the entropy formula for a few common parametric distributions. The larger the entropy of a distribution, the larger the uncertainties it describes. For instance, the entropy of a univariate Gaussian distribution increases with the log-variance and is independent to its mean. Similarly, the entropy of a uniform distribution on a compact set is the log-volume of this set. The smaller the volume, the smaller the uncertainties on the associated random variable. These examples also show that the differential entropy is not always positive, unlike the Shannon entropy that applies to discrete random variables. For a univariate Gaussian distribution, when the variance tends to 0, the differential entropy tends to $-\infty$. For a uniform distribution on a compact set, the differential entropy tends to $-\infty$ when the volume of the set tends to 0.

<table><caption>Table A.1 – Differential entropy for a few common distributions. The Vol set operator corresponds to the volume of a set. The det function corresponds to the determinant of a matrix.</caption><thead><tr><th></th><th>distribution</th><th>differential entropy h</th></tr></thead><tbody><tr><td rowspan="3">univariate<br>θ ∈ ℝ</td><td>θ ~ N(μ, σ²), μ ∈ ℝ, σ > 0</td><td>1/2 ln(2πe) + 1/2 ln σ²</td></tr><tr><td>θ ~ Unif(a, b), a < b</td><td>ln(b - a)</td></tr><tr><td>θ ~ Lognormal(μ, σ²), μ ∈ ℝ, σ > 0</td><td>μ + 1/2 ln(2πe σ²)</td></tr><tr><td rowspan="2">multivariate<br>θ ∈ ℝ<sup>D</sup>, D > 1</td><td>θ ~ N(μ, Σ)</td><td>D/2 ln(2πe) + 1/2 ln det(Σ)</td></tr><tr><td>θ ~ Unif(C), C ⊂ ℝ<sup>D</sup></td><td>ln Vol C</td></tr></tbody></table>

**Conditional differential entropy** – The conditional differential entropy is defined for an observation variable $y \in \mathbb{R}^L$ and a parameter $\theta \in \mathbb{R}^D$,

$$h(\theta|y) = E_{(\theta,y)}[-\ln \pi(\theta|y)] = -\int \ln \pi(\theta|y) \pi(\theta, y) d\theta dy. \quad (A.5)$$

One essential property of the conditional differential entropy is that

$$h(\theta|y) = - \int \ln \frac{\pi(\theta, y)}{\pi(y|\theta^*)} \pi(\theta, y) d\theta dy = h(\theta, y) - h(y). \quad (A.6)$$
This means that evaluating the conditional entropy is equivalent to computing the differential entropy of two linked distributions. Note that differential entropy and conditional differential entropy are linked to the mutual information $MI(\theta, y)$, as

$$MI(\theta, y) = h(\theta) - h(\theta|y). \qquad (A.7)$$

Unlike differential entropy and conditional differential entropy, the mutual information is always positive, as $h(\theta) \ge h(\theta|y)$.

**Fano bound** – The conditional differential entropy provides a lower bound on the general MSE from Eq. A.3. This lower bound is called the Fano bound, and reads

$$\mathbb{E}_{(\theta, y)} [\|\theta - \hat{\theta}(y)\|_2^2] \ge \frac{D}{2\pi e} \exp \left( \frac{2}{D} h(\theta|y) \right). \qquad (A.8)$$

This entropy-based lower bound on the MSE does not require the choice of an estimator nor its computation for all $y$ drawn from the joint distribution $\pi(\theta, y)$. It is only based on a summary statistic of the joint distribution $\pi(\theta, y)$, the conditional differential entropy. Besides, this lower bound is quite tight. The inequality turns into an equality in the case where the posterior distribution $\theta$ is Gaussian, and there exists $\lambda > 0$ such that $C = \lambda I_D$. Appendix A.A derives this lower bound. In essence, this derivation relies on the fact the Gaussian distribution maximizes the entropy, and on the arithmetic mean geometric mean inequality

$$[\det C]^{\frac{1}{D}} \le \frac{\mathrm{Tr} C}{D} \qquad (A.9)$$

Finally, this bound can be shown to be tighter than other lower bounds such as the Bayesian Cramér-Rao bound (Aras et al., 2019).

### A.1.3 Selection of the best set of lines as a discrete optimization problem

Unlike the simple cases listed in Table A.1, the distributions on $(\theta, y)$ and $y$ are not parametric. Therefore, their differential entropies do not have simple closed-form expressions. Both $K, D \lesssim 10$, thus a grid-based quadrature of the differential entropies in Eq. A.6 would be possible. However, for simplicity, we resort to MC estimators $\hat{h}^{(T_{MC})}(\theta^*, S)$

$$\hat{h}^{(T_{MC})}(\theta^*, S) = -\frac{1}{T_{MC}} \sum_{t=1}^{T_{MC}} \ln \frac{\pi(\theta^{(t)}, \mathbf{y}^{(t)})}{\pi(\mathbf{y}^{(t)})} \qquad (A.10)$$

$$= \frac{1}{T_{MC}} \sum_{t=1}^{T_{MC}} \ln \pi(\mathbf{y}^{(t)}) - \frac{1}{T_{MC}} \sum_{t=1}^{T_{MC}} \ln \pi(\theta^{(t)}, \mathbf{y}^{(t)}), \qquad (A.11)$$

with $\mathbf{y}^{(t)} \sim \pi(\mathbf{y}|\theta^*)$ and $\theta^{(t)} \sim \pi(\theta|\mathbf{y}^{(t)})$ for all $t$. The best set of $K$ lines for a given $\theta^*$ is the solution of the discrete combinatorial optimization problem

$$S_K^*(\theta^*) \in \arg \min_{S \in S_K} \hat{h}^{(T_{MC})}(\theta^*, S), \qquad (A.12)$$

with $S_K \subset 2^{[1,L]}$ the set of sets of $K$ lines.

## A.2 Illustration: a Gaussian and linear case

In this section, we consider a multivariate case in which the noise model is Gaussian, the forward model is affine and the prior on $\theta$ is also Gaussian. In this simple case, the MSE can be written in closed form. The MSE is compared with the Fano bound to visualize how tight this
bound is.

The considered simple inverse problem consists in estimating a physical parameter vector $\theta \in \mathbb{R}^D$ from an observation vector $y \in \mathbb{R}^L$. The forward model is set to an affine function $\theta \mapsto A\theta + b$, with $A \in \mathbb{R}^{L \times D}$ and $b \in \mathbb{R}^L$. The considered noise $\varepsilon$ is Gaussian additive and uncorrelated. We consider the following observation model, that involves

$$y = A\theta + b + \varepsilon, \text{ with } \varepsilon \sim \mathcal{N}(0, \sigma^2 I_L). \qquad (A.13)$$

The prior is set to a zero-mean Gaussian distribution

$$\theta \sim \mathcal{N}(0, \tau^2 I_D). \qquad (A.14)$$

In this setting, this prior is conjugate. The posterior distribution can thus be written simply by combining Eq. A.13 and Eq. A.14. It is given by

$$\theta|y \sim \mathcal{N}(\mu, C), \qquad (A.15)$$

with

$$\left\{ \begin{aligned} C &= \operatorname{Cov}(\theta|y) = \left[ \frac{1}{\sigma^2} A^T A + \frac{1}{\tau^2} I_D \right]^{-1} = \tau^2 \left[ I_D + \frac{\tau^2}{\sigma^2} A^T A \right]^{-1}, \\ \mu &= \mathbb{E}[\theta|y] = \frac{1}{\sigma^2} C A^T (y - b). \end{aligned} \right. \qquad (A.16)$$

The posterior expectation $\mu = \mathbb{E}[\theta|y]$ is the estimator $\hat{\theta}(y)$ that minimizes the MSE. With this estimator, the MSE reduces to

$$\mathbb{E}_{(\theta, y)} \left[ \| \theta - \hat{\theta}(y) \|^2 \right] = \operatorname{Tr} C = \tau^2 \operatorname{Tr} \left( \left[ I_D + \frac{\tau^2}{\sigma^2} A^T A \right]^{-1} \right). \qquad (A.17)$$

In the following, we derive the Fano lower bounds on the MSE to compare how tight it is.

**Fano Bound** – The computation of the Fano bound relies on $h(\theta|y) = h(\theta, y) - h(y)$. As $(\theta, y)$ and $y$ are Gaussian random variables,

$$\left\{ \begin{aligned} h(\mathbf{y}) &= \frac{1}{2} \ln \left[ (2\pi e)^L \det(\sigma^2 \mathbf{I}_L + \tau^2 \mathbf{A} \mathbf{A}^T) \right], \\ h(\theta, \mathbf{y}) &= \frac{1}{2} \ln \left( (2\pi e)^{D+L} \tau^{2D} \sigma^{2L} \right), \end{aligned} \right. \qquad (A.18)$$

and

$$h(\theta|y) = \frac{1}{2} \ln \left[ (2\pi e)^D \frac{\tau^{2D} \sigma^{2L}}{\det(\sigma^2 \mathbf{I}_L + \tau^2 \mathbf{A} \mathbf{A}^T)} \right]. \qquad (A.19)$$

The Fano Bound reads

$$\mathbb{E}_{(\theta, y)} \left[ \| \theta - \hat{\theta}(y) \|^2 \right] \geq \frac{D}{2\pi e} \exp \left( \frac{2}{D} h(\theta|y) \right) = D \frac{\tau^2 \sigma^{2(L/D)}}{\det(\sigma^2 \mathbf{I}_L + \tau^2 \mathbf{A} \mathbf{A}^T)^{(1/D)}}. \qquad (A.20)$$

The Weinstein-Aronszajn identity permits to simplify the comparison with the true MSE. It reads

$$\det(\sigma^2 \mathbf{I}_L + \tau^2 \mathbf{A} \mathbf{A}^T) = \sigma^{2(L/D)} \det \left( \mathbf{I}_D + \frac{\tau^2}{\sigma^2} \mathbf{A}^T \mathbf{A} \right). \qquad (A.21)$$
One can thus rewrite Eq. A.20 to

$$ \mathbb{E}_{(\theta, y)} [\|\theta - \hat{\theta}(y)\|^2] \geq \tau^2 D \det \left( \left[ \mathbf{I}_D + \frac{\tau^2}{\sigma^2} \mathbf{A}^T \mathbf{A} \right]^{-1} \right)^{(1/D)} . \quad (A.22) $$

This bound is very similar to the true MSE. The main difference comes from the arithmetic mean geometric mean inequality.

The Fano bound attains the true MSE if and only if the posterior on $\theta$ is Gaussian, the estimator is the posterior expectation, and the covariance matrix $C = \lambda \mathbf{I}_D$ for some $\lambda \ge 0$. The first two conditions are verified in this linear and Gaussian simple illustration. For the last one, considering $C$ from Eq. A.16,

$$ \begin{aligned} \exists \lambda \ge 0, \quad & \mathbf{C} = \lambda \mathbf{I}_D \Leftrightarrow \exists \lambda_2 \ge 0, \quad \mathbf{A}^T \mathbf{A} = \lambda_2 \mathbf{I}_D \\ & \Leftrightarrow \mathbf{A} = 0 \text{ or } \exists \lambda_2 > 0, \quad \frac{\mathbf{A}^T}{\sqrt{\lambda_2}} \frac{\mathbf{A}}{\sqrt{\lambda_2}} = \mathbf{I}_D \end{aligned} $$

The first case, $\mathbf{A} = 0$, implies that the observation model does not provide any information on $\theta$, and that the posterior distribution on $\theta$ is the same as the prior distribution. The second case can only be satisfied if $D \le L$. Indeed, if $D > L$, it is impossible for $\mathbf{A}^T \mathbf{A}$ to be full rank.

**Special cases** – Table A.2 summarizes the formulae obtained in this case for the true MSE and the Fano bound. The first two formulae illustrate the role of the arithmetic mean geometric mean inequality (Eq. A.9) in the relevance of the Fano bound. The Fano bound equals exactly the true MSE when the prior is dominant, when the matrix $\mathbf{A}$ is square with constant diagonal, or when both the observation and physical parameter vectors are reduced to scalars. The last two cases show that the linear map $\mathbf{A}$ minimizing the true MSE is the largest gradient $a$. In the general case, the most informative linear map $\mathbf{A}$ has the largest singular values. To obtain large singular values, the lines of $\mathbf{A}$ should have large norms and be as orthogonal as possible.

<table><caption>Table A.2 – Comparison of true MSE and Fano bound in Gaussian linear inverse problem</caption><thead><tr><th>Case</th><th>True MSE<br>(Eq. A.17)</th><th>Fano bound<br>(Eq. A.22)</th></tr></thead><tbody><tr><td>General</td><td>τ<sup>2</sup> Tr (<span style="text-decoration: overline">[<b>I</b><sub>D</sub> + <span style="text-decoration: overline">τ<sup>2</sup> <b>A</b><sup>T</sup><b>A</b></span>]<sup>-1</sup>)</td><td>τ<sup>2</sup> D det (<span style="text-decoration: overline">[<b>I</b><sub>D</sub> + <span style="text-decoration: overline">τ<sup>2</sup> <b>A</b><sup>T</sup><b>A</b></span>]<sup>-1</sup>)<sup>(1/D)</sup>)</td></tr><tr><td>σ → 0</td><td>σ<sup>2</sup> Tr (<span style="text-decoration: overline">[<b>A</b><sup>T</sup><b>A</b></span>]<sup>-1</sup>)</td><td>σ<sup>2</sup> D det (<span style="text-decoration: overline">[<b>A</b><sup>T</sup><b>A</b></span>]<sup>-1</sup>)<sup>(1/D)</sup>)</td></tr><tr><td>||<b>A</b>||<sub>F</sub> → 0<br>or τ → 0</td><td>Dτ<sup>2</sup></td><td>Dτ<sup>2</sup></td></tr><tr><td>D = L<br>and <b>A</b> = a<b>I</b><sub>D</sub></td><td><div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">τ<sup>2</sup>σ<sup>2</sup></div><div style="padding: 0 5px;">a<sup>2</sup>τ<sup>2</sup> + σ<sup>2</sup></div></div></td><td><div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">τ<sup>2</sup>σ<sup>2</sup></div><div style="padding: 0 5px;">a<sup>2</sup>τ<sup>2</sup> + σ<sup>2</sup></div></div></td></tr><tr><td>D = L = 1</td><td><div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">τ<sup>2</sup>σ<sup>2</sup></div><div style="padding: 0 5px;">a<sup>2</sup>τ<sup>2</sup> + σ<sup>2</sup></div></div></td><td><div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">τ<sup>2</sup>σ<sup>2</sup></div><div style="padding: 0 5px;">a<sup>2</sup>τ<sup>2</sup> + σ<sup>2</sup></div></div></td></tr></tbody></table>

## A.3 Conclusion

In this chapter, we presented a variable selection algorithm that permits identifying the lines that lead to the lowest uncertainties in inference. The MSE is used as a quantitative criterion to
quantify this uncertainty. As it cannot be evaluated, we resort to comparing lower bounds on the
MSE. The Fano bound is a tight lower bound. It is based on the conditional differential entropy
of a joint distribution $\pi(\theta, y)$, which can be evaluated efficiently with an MC estimator.

Future work involves adapting the code of the sampling algorithm proposed in Chapter 5 to
include a distribution on the observation $\pi(y|\theta^*)$. Then, it will be applied to compare sets of lines
predicted by the Meudon PDR code. The full set of $L = 5375$ lines can be considerably reduced to
a few hundreds by only considering lines that are observable, i.e., that reach intensities larger than
$\sigma_a$. However, even with this smaller number of lines, the number of sets to evaluate becomes very
large for $K \ge 2$. A heuristic approach may be needed to perform selection. Such methods include
e.g., greedy algorithms, stepwise forward selection or backward elimination methods (Shalev-
Shwartz and Ben-David, 2014, chapter 25, section 25.1), and meta-heuristics (Dreo et al., 2006).
# Appendix A.A Derivation of the multivariate Fano Bound

## The scalar case

This subsection's goal is to serve as a simple case to the next one. It is essentially a more detailed version (particularly about hypothesis) than Cover and Thomas (2006, theorem 8.6.6).

**Theorem A.A.1** (Estimator MSE and differential entropy in 1D case). For any scalar random variable $\theta \in \mathbb{R}$ and any estimator $\hat{\theta}$,

$$ \mathbb{E}[(\theta - \hat{\theta})^2] \geq \frac{1}{2\pi e} e^{2h(\theta)}, \qquad (A.23) $$

with equality if and only if $\theta$ is Gaussian and $\hat{\theta} = \mathbb{E}[\theta]$

*Proof.* Let $\hat{\theta}$ be any estimator of $\theta$. Then

$$ \mathbb{E}[(\theta - \hat{\theta})^2] \geq \min_{\hat{\theta}} \mathbb{E}[(\theta - \hat{\theta})^2] = \mathbb{E}[(\theta - \mathbb{E}[\theta])^2] = \mathrm{Var}(\theta). \qquad (A.24) $$

Then, to get a relation between the variance of $\theta$ and its entropy, one uses the fact that for scalar random variables with given variance, the Gaussian distribution has maximum entropy. Let $\sigma^2 = \mathrm{Var}(\theta)$ and $Z \sim \mathcal{N}(\mu, \sigma^2)$, Then,

$$ h(\theta) \leq h(Z) = \frac{1}{2} \ln(2\pi e \sigma^2), \qquad (A.25) $$

which gives

$$ \mathrm{Var}(\theta) = \sigma^2 \geq \frac{1}{2\pi e} e^{2h(\theta)}, \qquad (A.26) $$

with equality if and only if $\theta$ is Gaussian.

□

**Corollary A.A.1.1.** Given side information $y \in \mathbb{R}^L$ and estimator $\hat{\theta}(y)$, it follows that

$$ \mathbb{E}_{(\theta, y)}[(\theta - \hat{\theta}(y))^2] \geq \frac{1}{2\pi e} e^{2h(\theta|y)}, \qquad (A.27) $$

with equality if and only if $\theta$ is Gaussian and the estimator $\hat{\theta}(y)$ is the posterior expectation, i.e., $\hat{\theta}(y) = \mathbb{E}[\theta|y]$.

*Proof.* Let $\hat{\theta}(y)$ be any estimator of $\theta$ from observation $y$. Then

$$ \mathbb{E}_{(\theta, y)}[(\theta - \hat{\theta}(y))^2] \geq \min_{\hat{\theta}} \mathbb{E}_{(\theta, y)}[(\theta - \hat{\theta}(y))^2] = \mathbb{E}_{(\theta, y)}[(\theta - \mathbb{E}[\theta|y))^2] \qquad (A.28) $$

$$ = \mathbb{E}_y[\mathrm{Var}(\theta|y)] \qquad (A.29) $$

Note that $\mathbb{E}_y[\mathrm{Var}(\theta|y)]$ is the average variance of $\theta|y$ (when averaging over $y$). Likewise, by definition, $h(\theta|y)$ is the average entropy of $\theta|y$ (averaged over $y$). Therefore, by applying the theorem above, one gets the expected inequality. Like in the theorem above, equality is attained when $\theta$ is Gaussian and the estimator $\hat{\theta}(y)$ is the posterior expectation, ie $\hat{\theta}(y) = \mathbb{E}[\theta|y]$. □

No assumption over the distribution or dimensionality of $y$ is necessary. This implies that this bound works whether the relation between $\theta$ and $y_{nl}$ is linear or not.

In the general case, we can't compute $h(\theta|y)$ analytically, and we need to estimate it with an MC estimator. When the pdfs can be evaluated – which is the case in this chapter – then the standard estimator from Eq. A.10 can be used. When the pdfs cannot be evaluated, the entropy can be evaluated with the estimator from Kraskov et al. (2004).