# Sampling methods for statistical inference of non-linear inverse problems: spatial distribution of physico-chemical properties of the interstellar medium

Pierre Palud

► To cite this version:

Pierre Palud. Sampling methods for statistical inference of non-linear inverse problems: spatial distribution of physico-chemical properties of the interstellar medium. Signal and Image Processing. Centrale Lille, 2023. English. NNT: tel-04424965v2

HAL Id: tel-04424965

https://theses.hal.science/tel-04424965v2

Submitted on 1 Mar 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# THÈSE

présentée en vue d'obtenir le grade de

DOCTEUR

dans la spécialité

« AUTOMATIQUE, GÉNIE INFORMATIQUE, TRAITEMENT DU SIGNAL ET DES IMAGES »

par

Pierre PALUD

DOCTORAT DÉLIVRÉ PAR CENTRALE LILLE

## Sampling methods for statistical inference of non-linear inverse problems : spatial distribution of physico-chemical properties of the interstellar medium

Méthodes d'échantillonnage pour l'inférence statistique de problèmes inverses non linéaires : distribution spatiale des propriétés physico-chimiques du milieu interstellaire

Soutenue le 7 décembre 2023 devant le jury d'examen :

<table><tr><td>Président, Rapporteur</td><td>JÉRÔME IDIER</td><td>DR CNRS, LS2N (Nantes)</td></tr><tr><td>Rapporteur</td><td>OLIVIER BERNÉ</td><td>DR CNRS, RAP (Toulouse)</td></tr><tr><td>Examinatrice</td><td>FLORENCE FORBES</td><td>Directrice de recherche, INRIA Grenoble Rhone-Alpes</td></tr><tr><td>Examinatrice</td><td>ÉMILIE HABART</td><td>Maîtresse de conférences, Univertité Paris-Saclay - IAS</td></tr><tr><td>Examinateur</td><td>EMERIC BRON</td><td>Astronome adjoint, Observatoire de Paris</td></tr><tr><td>Examinateur</td><td>PIERRE-ANTOINE THOUVENIN</td><td>Maître de conférences, Centrale Lille Institut</td></tr><tr><td>Directeur de Thèse</td><td>PIERRE CHAINAIS</td><td>Professeur des universités, Centrale Lille Institut</td></tr><tr><td>Directeur de Thèse</td><td>FRANCK LE PETIT</td><td>Astronome, Observatoire de Paris</td></tr></table>

THÈSE PRÉPARÉE DANS LE **Laboratoire CRISTAL ET L'Observatoire de Paris**
ECOLE DOCTORALE MADIS 631
# Résumé long

## Méthodes d'échantillonnage pour l'inférence statistique de problèmes inverses non linéaires : distribution spatiale des propriétés physico-chimiques du milieu interstellaire

Cette thèse est une collaboration entre statisticiens et astrophysiciens. L'objectif est de développer des outils statistiques pour l'étude de la physique du milieu interstellaire (MIS). Le milieu interstellaire est un milieu très diffus qui remplit l'immense volume entre les objets célestes tels que les étoiles et les trous noirs au sein d'une galaxie. Il rassemble une grande variété d'environnements. La plupart du volume du MIS est chaude, ionisée et diffuse. À l'inverse, la plupart de la masse du MIS est froide, neutre et suffisamment dense pour que l'hydrogène — qui représente 90% du MIS — soit sous forme moléculaire. Ces derniers environnements forment des régions appelées nuages moléculaires.

**Pourquoi étudier le milieu interstellaire ?** – L'étude du MIS soulève des questions fondamentales comme la formation d'étoiles et le développement de la complexité moléculaire menant aux molécules pré-biotiques. Cette thèse se concentre sur le processus de formation d'étoiles, en particulier sur les effets de rétroaction d'une étoile sur les nuages moléculaires. Les étoiles naissent de l'effondrement gravitationnel de parties de nuages moléculaires. Une nouvelle étoile influence son nuage parent via son champ radiatif ultraviolet (UV) et via des vents solaires. Une étoile massive mourante influence également son nuage moléculaire parent à travers une explosion en supernova. L'impact global de ces effets de rétroaction sur ce qu'il reste du nuage parent est partiellement compris à ce jour. Ces effets pourraient dissiper le nuage, ce qui diminuerait sa capacité à former de nouvelles étoiles. Ils pourraient également comprimer localement les régions du nuage qui lui sont proches, et donc favoriser l'apparition de nouvelles étoiles.

Dans cette thèse, l'attention est portée aux nuages moléculaires illuminés et chauffés par des étoiles massives voisines émettant des photons UV. La couche de surface de tels nuages, où le champ radiatif UV chauffe et dissocie le gaz moléculaire, est appelée région de photodissociation (PDR). L'étude des conditions physiques dans une PDR permet de mieux comprendre l'effet de rétroaction radiative d'une étoile massive sur son nuage parent. Les ions, atomes et molécules présents dans une PDR se refroidissent en partie via des émissions de photons observables, à des fréquences correspondant à des transitions quantiques.

**Des observations à l'estimation des paramètres physiques** – Cette thèse analyse des cartes hyperspectrales de PDRs dans les domaines infrarouge lointain et radio, obtenues par des télescopes spatiaux ou terrestres. Ces cartes hyperspectrales sont réduites à des cartes multispectrales que l'on note $\mathbf{Y} = (y_n) \in \mathbb{R}^{N \times L}$, où chaque pixel $y_n$ contient l'intensité intégrée de $L \sim 5-30$ raies d'émission. Les cartes multispectrales considérées dans cette thèse contiennent de $N = 1$ à $N = 10\,000$ pixels. Ce type de carte est voué à devenir de plus en plus commun dans les années à venir, grâce à l'arrivée de nouveaux instruments comme le télescope spatial James Webb (JWST).

Ces cartes d'intensités $\mathbf{Y}$ peuvent être comparées avec les prédictions d'un modèle numérique du MIS. De nombreux modèles existent, pour les PDRs et pour d'autres types d'environnements. Dans cette thèse, nous utilisons le code PDR de Meudon (Le Petit et al., 2006), que l'on note $f$ et qui est état de l'art. Ce code modélise le nuage comme un bloc de gaz plan-parallèle à l'état stationnaire, et y résout les équations du transfert radiatif, de la chimie et de l'équilibre
thermique. Un modèle numérique de ce type peut calculer des cartes d'intensité intégrée de raies à partir de cartes de paramètres physiques. Nous notons ces cartes de paramètres physiques $\Theta = (\theta_n) \in \mathbb{R}^{N \times D}$, où chaque vecteur $\theta_n$ contient $D \lesssim 10$ paramètres comme la densité du gaz, la pression thermique, l'intensité du champ radiatif incident ou l'épaisseur totale du nuage le long de la ligne de visée. Les paramètres $\Theta$ vivent donc en haute dimension, mais présentent une structure simple en pixels. Cette thèse vise à estimer des cartes de paramètres physiques $\Theta$ à partir d'une carte d'observation $\mathcal{Y}$ et du code PDR de Meudon $f$.

**Modèle statistique** – Ce problème est une instance d'une classe générale de problèmes statistiques, appelés problèmes inverses. Avec une approche bayésienne, poser un tel problème nécessite de définir deux distributions : une vraisemblance et une distribution a priori.

La vraisemblance encode un modèle d'observation, qui combine un modèle direct et un modèle de bruit. Une évaluation du code PDR de Meudon dure de 6 heures à quelques jours, ce qui empêche son utilisation directe pour de l'inférence. Une procédure générale est proposée pour construire une approximation de modèles numériques. Elle exploite un réseau de neurones spécifique et surpasse les méthodes d'interpolation – communes dans les études du MIS – en termes de précision, de poids mémoire et de temps d'évaluation. Elle est appliquée au code PDR de Meudon pour 5375 raies ioniques, atomiques et moléculaires. Le réseau de neurones obtenu est donc utilisable dans une procédure d'inférence quelle que soit la combinaison de raies observée.

Le modèle de bruit tient compte d'autant de sources d'incertitudes que possible. Il fait intervenir un bruit multiplicatif lognormal et un bruit additif gaussien. Le bruit multiplicatif encode une erreur de calibration et une incertitude sur la validité du code PDR de Meudon. Le bruit additif inclut du bruit thermique ou une erreur sur le comptage de photons, selon la gamme de fréquences. La vraisemblance associée à ce mélange de bruit n'a pas de forme simple. Nous en proposons une approximation avec contrôle de l'erreur. Le problème inverse défini avec la vraisemblance seule est mal posé : plusieurs solutions sont compatibles avec les observations.

La distribution a priori mélange deux termes afin de lever partiellement ces dégénérescences sur les solutions. Des intervalles de validité sont considérés sur chaque paramètre physique. Une régularisation spatiale améliore les estimations en favorisant des cartes plus lisses. Cette régularisation intervient principalement dans les régions à faibles rapport signal sur bruit, qui correspondent à du gaz diffus où les paramètres évoluent lentement sur une carte. De fait, ces régions sont les plus affectées par le caractère mal posé du problème inverse. Dans les régions à haut rapport signal sur bruit, cette régularisation spatiale sera dominée par la vraisemblance et y jouera donc un rôle moindre.

La distribution a posteriori combine vraisemblance et distribution a priori. C'est la distribution exploitée pour inférer les paramètres physiques. Cette distribution présente trois difficultés principales. Elle est non log-concave et potentiellement multimodale, puisque le code PDR de Meudon est non linéaire. Explorer cette distribution requiert donc une capacité d'échapper à un mode local mais non global. Le log-posterior est non gradient Lipschitz, puisque le code PDR de Meudon couvre de nombreux ordres de grandeur. Cette condition de régularité est requise pour beaucoup d'algorithmes classiques. Enfin, les paramètres physiques $\Theta$ vivent en haute dimension. Exploiter cette distribution requiert donc des algorithmes pouvant passer à l'échelle. Peu de cas dans la littérature présentent ces trois difficultés simultanément. Le problème inverse considéré est donc général et difficile.

**Inférence bayésienne** – Une nouvelle méthode d'inférence traitant ces trois difficultés est proposée. Elle s'appuie sur une approche par échantillonnage pour fournir des quantifications d'incertitudes en plus d'estimateurs ponctuels de cartes de paramètres physiques. Cette quantification d'incertitude est critique afin de compenser l'absence de vérité terrain, inhérente à l'astrophysique.

L'algorithme proposé de Monte Carlo par chaîne de Markov (MCMC) (Robert et Casella, 2004, chapitres 6 et 7) combine deux noyaux d'échantillonnage. Le premier noyau identifie les minima locaux dans l'espace des paramètres. Il combine un algorithme dit "Metropolis à essais multiples"
(MTM) (Martino, 2018) et un algorithme de "Gibbs chromatique" (Gonzalez et al., 2011), qui exploite la structure de $\Theta$. Le second noyau explore efficacement les minima locaux en exploitant des informations sur leur géométrie locale. Il étend l'algorithme de Langevin avec ajustement de Metropolis (MALA) (Roberts et Stramer, 2002) avec l'utilisation d'un préconditionneur (Girolami et Calderhead, 2011). Le préconditionneur adopté est RMSProp (Tieleman et Hinton, 2012), très utilisé pour l'entraînement de réseaux de neurones. Chaque noyau est efficace en grande dimension grâce à l'exploitation de la structure de $\Theta$ et au préconditionnement.

Le processus d'inférence proposé se base sur l'hypothèse de validité du modèle. Cette hypothèse est testée pour s'assurer de la pertinence des résultats. Nous avons étendu un test bayésien par valeur-$p$. Ce test porte sur la capacité du modèle d'observation considéré à reproduire les observations $Y$.

**Applications à des données synthétiques** – La méthode proposée est validée à l'aide de trois jeux de données synthétiques. Un mélange de gaussiennes en deux dimensions illustre les propriétés des deux noyaux proposés. Les gaussiennes étant isolées, cette distribution est difficile à explorer pour un algorithme MCMC classique. L'algorithme proposé parvient à visiter correctement chaque mode.

Le problème de localisation de détecteurs (Ihler et al., 2005) montre la capacité de l'algorithme à exploiter efficacement la structure des paramètres $\Theta$. Ce problème à 16 dimensions est lui aussi multimodal. L'algorithme proposé identifie tous les modes de la distribution. Il se compare favorablement à d'autres algorithmes état de l'art dédiés aux distributions multimodales.

Une carte multispectrale synthétique mais réaliste démontre la capacité de l'algorithme à résoudre des problèmes inverses sur de grandes cartes. Un nuage fictif est défini pour couvrir plusieurs types d'environnements physiques. Ce nuage correspond à une carte de paramètres théoriques $\Theta^*$ de $N = 4096$ pixels et $D = 4$ paramètres. Les observations sont générées à partir de cette carte théorique, du code PDR de Meudon et de réalisations du modèle de bruit. Elles contiennent $L = 10$ raies de $^{12}$CO. L'échantillonneur proposé mène à des estimations $\hat{\Theta}$ très proches des valeurs théoriques, avec une quantification des incertitudes associées. La cohérence physique des résultats est vérifiée.

**Applications à des données réelles** – Trois applications à des observations réelles de nuages moléculaires sont présentées. Les résultats de l'inférence sont analysés pour fournir des interprétations astrophysiques.

L'observation de NGC 7023 (Joblin et al., 2018) contient $N = 1$ pixel et $L = 17$ raies, dont des raies de $^{12}$CO (de $J = 11 \rightarrow 10$ à $J = 19 \rightarrow 18$), de $H_2$ (de S(0) à S(5)) et de $CH^+$ (de $J = 1 \rightarrow 0$ à $J = 3 \rightarrow 2$). Le test du modèle montre que le code PDR de Meudon parvient à reproduire simultanément toutes les raies observées. L'extinction visuelle n'est pas contrainte puisque les raies observées sont principalement émises depuis la surface du nuage, et n'en sondent donc pas la profondeur. A l'inverse, la pression thermique et l'intensité du champ radiatif sont bien contraints et les valeurs obtenues sont cohérentes avec des estimations précédentes.

L'observation de la nébuleuse de la Carène (Wu et al., 2018) contient $N = 176$ pixels et $L = 12$ raies, dont des raies de $^{12}$CO (de $J = 4 \rightarrow 3$ à $J = 13 \rightarrow 12$) et de C. Pour chaque pixel, le code PDR de Meudon reproduit simultanément toutes les raies observées. Dans Wu et al. (2018), les cartes de paramètres physiques estimées présentaient des variations spatiales brutales. Grâce à la régularisation spatiale et à l'approche par échantillonnage, nos reconstructions sont plus cohérentes physiquement. Les 3 régions présentes sur la carte y sont identifiables. Les valeurs estimées de pression de Wu et al. (2018) sont compatibles avec les incertitudes obtenues. En revanche, les intensités du champ radiatif obtenues avec notre approche sont plus faibles d'un facteur 4. Cette différence est justifiée par un changement de version du code PDR de Meudon et donc de la microphysique considérée.

Enfin, l'observation d'OMC-1 (Goicoechea et al., 2019) contient $N = 2475$ pixels et $L = 4$ raies, dont deux de $^{12}$CO, une de $CH^+$ et une de $HCO^+$. Ce jeu de données n'avait jusqu'à présent
jamais été étudié dans une procédure d'estimation, contrairement aux deux précédents. Cette carte inclut de nombreux environnements : plusieurs PDRs dont la barre d'Orion, un filament, Orion S, et une région HII. L'ensemble de la carte est illuminé par le cluster d'étoiles du Trapèze, qui émet un très fort rayonnement UV. Le test du modèle montre que le code PDR de Meudon parvient à reproduire simultanément toutes les raies observées pour tous les pixels, sauf quelques pixels isolés aux bords de la carte. Ces pixels mal reproduits correspondent à des régions de gaz diffus, et sont donc non essentiels pour l'analyse de l'effet de rétroaction radiative sur le gaz moléculaire. Les estimations obtenues sont lisses et permettent d'identifier les différentes régions. Les incertitudes obtenues sont restreintes, ce qui permet de tirer des conclusions sur les environnements observés. En particulier, nos inversions montrent pour la première fois que les pressions thermiques élevées, connues pour la barre d'Orion, se retrouvent dans les autres PDRs d'OMC-1.

Dans les applications présentées, au sein d'un même pixel, les incertitudes sur la pression et sur l'intensité du champ radiatif incident présentent une anti-corrélation. À l'inverse, entre différents pixels, la pression et l'intensité du champ radiatif incident estimés présentent une corrélation positive. Cette corrélation positive a déjà été décrite dans la littérature (Joblin et al., 2018; Wu et al., 2018). Elle semble indiquer que la rétroaction d'une étoile massive favorise l'apparition de nouvelles étoiles à partir du nuage parent. L'anti-corrélation dans les incertitudes démontre que cette corrélation positive entre différentes sources n'est pas un artefact, et a probablement une origine physique.
Acknowledgements

Today, a few weeks after the defense, I feel very proud of these three years of PhD. I feel
proud about the work I accomplished, the difficulties I overcame, and the expertise I developed.
I did not do all this by myself: many people contributed to this success. I want to acknowledge
and thank them here.

First and foremost, I want to thank my two supervisors, Pierre Chainais and Franck Le Petit, and my two co-advisors, Emeric Bron and Pierre-Antoine Thouvenin. Thank you for entrusting me with this project in the first place, and for all the time you dedicated to discuss, orient and review my work. Every comment you made pushed me towards better science and better structured, clearer and more rigorous papers and presentations.

Thank you in particular to Pierre-Antoine for the many (maybe too many?) times I came to your office asking "Do you have 5 minutes?" and ended up staying 45. Thank you as well for giving me the opportunity to teach some of your students. Thank you, Emeric and Franck, for the good moments of conviviality in and out of the office. By the way, your precious recommendations gave me back the pleasure of reading fiction!

Thank you to the other members of the jury, namely Jérôme Idier, Olivier Berné, Florence Forbes and Émilie Habart, for accepting to review my work. Thank you in particular, Jérôme Idier, for the rich discussion we had about multiple technical points. Your many comments contributed to improve this manuscript to its current version.

I divided my research time in two laboratories, the SIGMA team from the CRISTAL laboratory in Lille, and the LERMA laboratory in the Observatory of Paris, in the Meudon site. I am very grateful to the interns, PhD students, PostDocs and to the permanent staff that made these two working environments very enjoyable. For the time I spent in Lille, thank you Clémence, Étienne, Quentin, Jimmy, Maxime, Rony, Martin, Solène, Yoann, Jan, Barbara, Jérémie and Patrick for the coffee breaks, the PhD coffees, and the afterwork sessions. Thank you Markus, my PhD twin, for sharing this PhD journey with me. For the time I spent in Meudon, thank you, Raphaël, Vincent, Japhar, David and Nicolas. Thank you also to the administrative staff that helped me throughout these years, Sonia Waskowiak and Aurélien Denis.

In addition to these two laboratories, I had the luck to integrate the ORION-B consortium, led by Jérôme Pety and Maryvonne Gérin. Thank you, Jérôme and Maryvonne, for the effort you put on the coordination of this consortium and on the organization of the biannual Orion weeks. The collaborative atmosphere and the rich discussions that arose from our meetings were inspiring and led me to explore new ideas and methods, despite our different backgrounds and vocabularies. Besides, I believe the now traditional discretely-say-a-given-word-during-a-presentation game from the Orion weeks is a good indicator of the friendly mood in this consortium.

I thank all the members for the insightful exchanges, in particular Miriam and Javier for the OMC-1 data and our many discussions on this region, and Antoine R. and Pierre G. for your enthusiasm and feedback. And to Antoine Z.: I hope you manage to use BEETROOTS for your PhD, best of luck! Finally, a huge thank you to Lucas Einig! Thank you for onboarding with me on the model approximation and entropy adventures. Your expertise in PYTORCH and your implication turned these questions into fascinating research projects. I greatly enjoyed collaborating with you.

On a more personal side, thank you to my parents Pierre-Yves and Isabelle, my brother Florent, my sister Romane, and my fiancée Isabela. Thank you for your unconditional support, long before and during these three years, in the day-to-day life and in the big moments. I love you.
# Contents

Résumé long iii

Acknowledgements vii

Contents ix

List of Figures xiii

List of Tables xv

List of Algorithms xvii

Notation xix

Introduction 1

Scientific context 1

PhD project summary 2

Structure of the manuscript 6

List of publications 7

I Backgrounds: interstellar medium, statistical inference and interactions 9

1 Background on the interstellar medium (ISM) 11

1.1 The interstellar medium (ISM) 13

1.1.1 A short overview of the ISM 13

1.1.2 Star formation, stellar feedback and photodissociation regions (PDRs) 15

1.1.3 Chemical complexity and planet formation 17

1.1.4 PDRs: structure and physical parameters 18

1.2 Linking physical conditions and observations: astrophysical numerical models 20

1.2.1 Radiative transfer models 20

1.2.2 Astrochemical models 21

1.2.3 Dust models 22

1.2.4 Holistic models 22

1.3 The Meudon PDR code 23

1.3.1 Physics and chemistry taken into account 23

1.3.2 Input parameters and considered environments 25

1.3.3 Limitations 26

1.4 Conclusion 26

2 Background on Bayesian statistical modeling and inference 29

2.1 Bayesian statistical modeling 30

2.1.1 Main random variables and probability distributions 30

2.1.2 Estimators 32

2.2 Statistical inference 33

2.2.1 Evaluating estimators defined as solutions of optimization problems 34

2.2.2 Approximating integrals with Monte Carlo estimators 39
<table>
  <tbody>
    <tr>
      <td>2.3 Comparing and checking observation models</td>
      <td>49</td>
    </tr>
    <tr>
      <td>2.3.1 Model selection and evaluation</td>
      <td>50</td>
    </tr>
    <tr>
      <td>2.3.2 Model checking using Bayesian hypothesis testing</td>
      <td>52</td>
    </tr>
    <tr>
      <td>2.4 Conclusion</td>
      <td>54</td>
    </tr>
    <tr>
      <td>Appendix 2.A Samplers widespread in astrophysics dedicated to multimodal distributions</td>
      <td>56</td>
    </tr>
    <tr>
      <td>2.A.1 Adapting meta-heuristics to MCMC</td>
      <td>56</td>
    </tr>
    <tr>
      <td>2.A.2 Prior or parallel identification of the modes</td>
      <td>58</td>
    </tr>
    <tr>
      <td>2.A.3 Nested sampling</td>
      <td>59</td>
    </tr>
    <tr>
      <td><b>3 Review of statistical inference in ISM studies and position of our problem</b></td>
      <td><b>61</b></td>
    </tr>
    <tr>
      <td>3.1 Statistical modeling</td>
      <td>62</td>
    </tr>
    <tr>
      <td>3.1.1 Forward model: handling a numerical model</td>
      <td>62</td>
    </tr>
    <tr>
      <td>3.1.2 Noise model and forward model misspecification</td>
      <td>65</td>
    </tr>
    <tr>
      <td>3.1.3 Prior distributions and regularization functions</td>
      <td>70</td>
    </tr>
    <tr>
      <td>3.2 Statistical inference</td>
      <td>72</td>
    </tr>
    <tr>
      <td>3.2.1 Optimization-based inference</td>
      <td>72</td>
    </tr>
    <tr>
      <td>3.2.2 Sampling-based inference</td>
      <td>72</td>
    </tr>
    <tr>
      <td>3.3 Comparing and checking observation models</td>
      <td>74</td>
    </tr>
    <tr>
      <td>3.3.1 Model comparison</td>
      <td>74</td>
    </tr>
    <tr>
      <td>3.3.2 Model checking with posterior predictive assessment</td>
      <td>75</td>
    </tr>
    <tr>
      <td>3.4 Our observation model</td>
      <td>76</td>
    </tr>
    <tr>
      <td>3.4.1 Noise model on the original hyperspectral cube</td>
      <td>76</td>
    </tr>
    <tr>
      <td>3.4.2 From the hyperspectral cube to maps of integrated intensities</td>
      <td>79</td>
    </tr>
    <tr>
      <td>3.4.3 Forward model, observational effects and censorship</td>
      <td>80</td>
    </tr>
    <tr>
      <td>3.5 Modeling, inference and model assessment choices</td>
      <td>81</td>
    </tr>
    <tr>
      <td>3.6 Conclusion</td>
      <td>82</td>
    </tr>
    <tr>
      <td><b>II Solving inverse problems on ISM multiline maps</b></td>
      <td><b>83</b></td>
    </tr>
    <tr>
      <td><b>4 Fast simulations of photodissociation region models</b></td>
      <td><b>85</b></td>
    </tr>
    <tr>
      <td>4.1 Deriving emulators with interpolation or regression</td>
      <td>86</td>
    </tr>
    <tr>
      <td>4.1.1 Interpolation methods</td>
      <td>87</td>
    </tr>
    <tr>
      <td>4.1.2 Regression methods</td>
      <td>87</td>
    </tr>
    <tr>
      <td>4.2 Neural networks for regression</td>
      <td>88</td>
    </tr>
    <tr>
      <td>4.2.1 Generalities on neural networks</td>
      <td>88</td>
    </tr>
    <tr>
      <td>4.2.2 Fitting a neural network to a dataset</td>
      <td>89</td>
    </tr>
    <tr>
      <td>4.3 Dataset structure</td>
      <td>90</td>
    </tr>
    <tr>
      <td>4.4 Illustration: the log Rosenbrock function</td>
      <td>92</td>
    </tr>
    <tr>
      <td>4.5 Towards an emulator of the Meudon PDR code</td>
      <td>96</td>
    </tr>
    <tr>
      <td>4.5.1 Datasets</td>
      <td>96</td>
    </tr>
    <tr>
      <td>4.5.2 Comparison metrics</td>
      <td>98</td>
    </tr>
    <tr>
      <td>4.6 Designing and training adapted ANNs</td>
      <td>100</td>
    </tr>
    <tr>
      <td>4.6.1 Removing outliers from the training set</td>
      <td>100</td>
    </tr>
    <tr>
      <td>4.6.2 Exploiting correlations between line intensities</td>
      <td>102</td>
    </tr>
    <tr>
      <td>4.6.3 A polynomial transform to learn nonlinearities</td>
      <td>104</td>
    </tr>
    <tr>
      <td>4.6.4 Dense networks to reuse intermediate computations</td>
      <td>106</td>
    </tr>
    <tr>
      <td>4.7 Experiments: application to the Meudon PDR code</td>
      <td>107</td>
    </tr>
    <tr>
      <td>4.7.1 Performance analysis</td>
      <td>107</td>
    </tr>
    <tr>
      <td>4.7.2 Removing outliers is crucial</td>
      <td>109</td>
    </tr>
    <tr>
      <td>4.7.3 The importance of the polynomial feature augmentation</td>
      <td>109</td>
    </tr>
    <tr>
      <td>4.8 Conclusion</td>
      <td>110</td>
    </tr>
    <tr>
      <td>Appendix 4.A Automatic outlier detection procedure</td>
      <td>111</td>
    </tr>
  </tbody>
</table>
<table><tr><td>Appendix 4.B Content of clusters of lines</td><td>112</td></tr><tr><td><b>5 An MCMC algorithm for efficient inversion with quantified uncertainty</b></td><td><b>113</b></td></tr><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><tr><td><b>6 Application to real data</b></td><td><b>141</b></td></tr><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><tr><td><b>Conclusions and perspectives</b></td><td><b>171</b></td></tr><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><tr><td><b>A Entropy of probability distributions and line selection</b></td><td><b>177</b></td></tr><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><tr><td><b>Acronyms</b></td><td><b>187</b></td></tr><tr><td><b>References</b></td><td><b>189</b></td></tr><tr><td><b>Résumé court</b></td><td><b>205</b></td></tr><tr><td><b>Abstract</b></td><td><b>206</b></td></tr></table>
# List of Figures

<table>
  <tbody>
    <tr>
      <td>1</td>
      <td>Orion B giant molecular cloud</td>
      <td>2</td>
    </tr>
    <tr>
      <td>2</td>
      <td>General principle of an inference procedure</td>
      <td>3</td>
    </tr>
    <tr>
      <td>3</td>
      <td>Illustration of two types of solution multiplicity in an optimization problem</td>
      <td>4</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Replacing a loss function by a posterior distribution</td>
      <td>4</td>
    </tr>
    <tr>
      <td>1.1</td>
      <td>Illustration of the parsec definition</td>
      <td>12</td>
    </tr>
    <tr>
      <td>1.2</td>
      <td>Scales in astrophysics</td>
      <td>13</td>
    </tr>
    <tr>
      <td>1.3</td>
      <td>Hubble and JWST observation of the pillars of creation</td>
      <td>14</td>
    </tr>
    <tr>
      <td>1.4</td>
      <td>Observed flow of baryons in the Milky Way</td>
      <td>15</td>
    </tr>
    <tr>
      <td>1.5</td>
      <td>Key steps of the star formation process</td>
      <td>16</td>
    </tr>
    <tr>
      <td>1.6</td>
      <td>ALMA observations of two protoplanetary disks</td>
      <td>17</td>
    </tr>
    <tr>
      <td>1.7</td>
      <td>Structure of a PDR region</td>
      <td>18</td>
    </tr>
    <tr>
      <td>1.8</td>
      <td>Illustration of a spectral energy density for dust</td>
      <td>23</td>
    </tr>
    <tr>
      <td>1.9</td>
      <td>Summary of Meudon PDR code modeling components</td>
      <td>24</td>
    </tr>
    <tr>
      <td>2.1</td>
      <td>Illustration of the condition number of the Hessian matrix</td>
      <td>36</td>
    </tr>
    <tr>
      <td>2.2</td>
      <td>Graphs of conditional dependence and associated optimal coloring</td>
      <td>48</td>
    </tr>
    <tr>
      <td>2.3</td>
      <td>Illustration of proposed p-value in a Gaussian and a non-concave log-likelihood</td>
      <td>53</td>
    </tr>
    <tr>
      <td>2.4</td>
      <td>Illustration of nested sampling principle</td>
      <td>60</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Quality of Gaussian approximation of a lognormal distribution for low and high variance</td>
      <td>67</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Illustration of censorship in likelihood for Gaussian additive noise models</td>
      <td>69</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Examples of molecular emission lines spectra</td>
      <td>78</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Structure of a simple feedforward neural network</td>
      <td>89</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Illustration of methods to generate non-lattice datasets</td>
      <td>92</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>log Rosenbrock function</td>
      <td>93</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>Comparison of popular interpolation and regression methods on the log Rosenbrock function</td>
      <td>94</td>
    </tr>
    <tr>
      <td>4.5</td>
      <td>Illustration of an invalid integrated intensity due to a bistability in temperature profile</td>
      <td>97</td>
    </tr>
    <tr>
      <td>4.6</td>
      <td>Comparison of the error factor and the relative error</td>
      <td>99</td>
    </tr>
    <tr>
      <td>4.7</td>
      <td>Illustration of the impact of outliers on the mean, the 99th percentile and the max estimators of the error factor</td>
      <td>100</td>
    </tr>
    <tr>
      <td>4.8</td>
      <td>Comparison of graph of different error functions</td>
      <td>101</td>
    </tr>
    <tr>
      <td>4.9</td>
      <td>Matrix of absolute Pearson correlation coefficients among the L = 5375 lines predicted by Meudon PDR code</td>
      <td>103</td>
    </tr>
    <tr>
      <td>4.10</td>
      <td>Results of clustering of lines predicted by the Meudon PDR code</td>
      <td>105</td>
    </tr>
    <tr>
      <td>4.11</td>
      <td>Structure of a dense neural network</td>
      <td>107</td>
    </tr>
    <tr>
      <td>4.12</td>
      <td>Regularization function for simultaneous outlier identification and ANN training</td>
      <td>112</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>Illustration of the λ function</td>
      <td>117</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Comparison of samplers on a 2D Gaussian mixture model</td>
      <td>125</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>Illustration of the three-case test for model assessment</td>
      <td>128</td>
    </tr>
    <tr>
      <td>5.4</td>
      <td>Application of Bayesian p-value with three-case test on a simple case</td>
      <td>129</td>
    </tr>
    <tr>
      <td>5.5</td>
      <td>Comparison of samplers on the sensor localization problem</td>
      <td>131</td>
    </tr>
  </tbody>
</table>
<table><tr><td>5.6</td><td>Application of model assessment on sensor localization problem</td><td>132</td></tr><tr><td>5.7</td><td>Structure of the synthetic molecular cloud</td><td>133</td></tr><tr><td>5.8</td><td>Some observation maps of the astrophysical experiment</td><td>133</td></tr><tr><td>5.9</td><td>Inference results on a synthetic yet realistic astrophysics case</td><td>134</td></tr><tr><td>5.10</td><td>Optimization of the parameters of the likelihood approximation</td><td>138</td></tr><tr><td>6.1</td><td>Full inference process workflow</td><td>142</td></tr><tr><td>6.2</td><td>Posterior predictive assessment for NGC 7023</td><td>146</td></tr><tr><td>6.3</td><td>Inference results for NGC 7023</td><td>147</td></tr><tr><td>6.4</td><td>Observation of the Carina nebula</td><td>148</td></tr><tr><td>6.5</td><td>Observations of the Carina nebula</td><td>149</td></tr><tr><td>6.6</td><td>Model assessment on the Carina observations</td><td>150</td></tr><tr><td>6.7</td><td>Posterior predictive assessment two pixels of the Carina nebula</td><td>151</td></tr><tr><td>6.8</td><td>Inference results for the Carina nebula</td><td>153</td></tr><tr><td>6.9</td><td>Inference results for the Car I-S pixel of the Carina nebula</td><td>154</td></tr><tr><td>6.10</td><td>Inference results for the Car I-E pixel of the Carina nebula</td><td>155</td></tr><tr><td>6.11</td><td>Structure of the OMC-1 cloud</td><td>156</td></tr><tr><td>6.12</td><td>Observations of OMC-1</td><td>157</td></tr><tr><td>6.13</td><td>Model assessment for OMC-1</td><td>158</td></tr><tr><td>6.14</td><td>Posterior predictive assessment for four pixels of OMC-1</td><td>159</td></tr><tr><td>6.15</td><td>Inference results for OMC-1</td><td>161</td></tr><tr><td>6.16</td><td>Inference results for OMC-1: Orion bar pixel</td><td>162</td></tr><tr><td>6.17</td><td>Inference results for OMC-1: East PDR pixel</td><td>163</td></tr><tr><td>6.18</td><td>Inference results for OMC-1: North-East edge pixel</td><td>164</td></tr><tr><td>6.19</td><td>Inference results for OMC-1: North-West ridge edge pixel</td><td>165</td></tr><tr><td>6.20</td><td>Posterior predictive assessment for the Orion Bar</td><td>168</td></tr><tr><td>6.21</td><td>Inference results for the Orion Bar</td><td>169</td></tr></table>
# List of Tables

<table><tr><td>1</td><td>Reading options per topic</td><td>7</td></tr><tr><td>1.1</td><td>Average elementary abundances in the ISM in the Milky Way</td><td>14</td></tr><tr><td>1.2</td><td>Comparison of the main phases of the ISM</td><td>15</td></tr><tr><td>1.3</td><td>Secondary input parameters in the Meudon PDR code and their default values</td><td>25</td></tr><tr><td>4.1</td><td>Comparison of popular interpolation and regression methods on the log Rosenbrock function</td><td>95</td></tr><tr><td>4.2</td><td>Input parameters in the Meudon PDR code and structure of training dataset.</td><td>96</td></tr><tr><td>4.3</td><td>Performance of interpolation methods and of the proposed ANNs with and without removal of outlier from the training set</td><td>108</td></tr><tr><td>5.1</td><td>Samplers comparison on 2D Gaussian mixture model</td><td>124</td></tr><tr><td>5.2</td><td>Effective Sample Size (ESS) on the sensor localization problem</td><td>131</td></tr><tr><td>5.3</td><td>Reconstruction metrics and relative size of credible intervals for the astrophysics experiment</td><td>136</td></tr><tr><td>6.1</td><td>Inference results on NGC 7023</td><td>148</td></tr><tr><td>6.2</td><td>Inference results on the Orion Bar</td><td>168</td></tr><tr><td>A.1</td><td>Differential entropy for a few common distributions</td><td>179</td></tr><tr><td>A.2</td><td>Comparison of true MSE and Fano bound in Gaussian linear inverse problem</td><td>182</td></tr></table>
# List of Algorithms

<table><tr><td>2.1 Gradient descent (GD) algorithm</td><td>35</td></tr><tr><td>2.2 Metropolis-Hastings (MH) algorithm</td><td>43</td></tr><tr><td>2.3 Metropolis-within-Gibbs sampling</td><td>47</td></tr><tr><td>2.4 Multiple-Try Metropolis (MTM) algorithm</td><td>49</td></tr><tr><td>2.5 Affine-invariant sampler</td><td>57</td></tr><tr><td>2.6 Sequential Monte Carlo (SMC)</td><td>58</td></tr><tr><td>2.7 Nested sampling</td><td>60</td></tr><tr><td>5.1 PMALA kernel K<sub>1</sub> at step t</td><td>120</td></tr><tr><td>5.2 MTM kernel K<sub>2</sub> at step t</td><td>122</td></tr><tr><td>5.3 Proposed sampler: PMALA and MTM</td><td>123</td></tr><tr><td>5.4 MCMC sampling with p-value computation</td><td>127</td></tr><tr><td>5.5 Sampling from the smooth distribution in Eq. 5.39</td><td>139</td></tr></table>
Notation

$\Theta = (\theta_n)_{n=1}^N \in \mathbb{R}^{N \times D}$ set of $N$ physical parameter vectors $\theta_n$. This is the parameter to be inferred in an inverse problem. When $N=1$, this set is a vector and the notation is simplified to $\theta \in \mathbb{R}^D$. When $N=D=1$, this set is a scalar and the notation is simplified to $\theta \in \mathbb{R}$. When computing gradients with respect to $\Theta$, it is handled as a vector with conversion in lexicographic order, such that $\Theta \in \mathbb{R}^{ND}$.

$\theta \in \mathbb{R}^D$ physical parameter vector. A physical parameter $\theta_d \in \mathbb{R}$ can be for instance a thermal pressure, a visual extinction, or the intensity of a radiative field.

$\hat{\Theta} = (\theta_n)_{n=1}^N \in \mathbb{R}^{N \times D}$ generic notation for an estimator of $\Theta$. When $N=1$, this set is a vector and the notation is simplified to $\hat{\theta} \in \mathbb{R}^D$. When $N=D=1$, this set is a scalar and the notation is simplified to $\hat{\theta} \in \mathbb{R}$.

$\mathbf{Y} = (\mathbf{y}_n)_{n=1}^N \in \mathbb{R}^{N \times L}$ set of $N$ individual observation vectors $\mathbf{y}_n$, from which $\mathbf{\Theta}$ is inferred. When $N=1$, this set is a vector and the notation is simplified to $\mathbf{y} \in \mathbb{R}^L$. When $N=L=1$, this set is a scalar and the notation is simplified to $\mathbf{y} \in \mathbb{R}$.

$\mathbf{y} \in \mathbb{R}^L$ individual observation vector. An observation element $\mathbf{y}_l \in \mathbb{R}$ can be e.g., the integrated intensity of an ionic, atomic or molecular emission line.

$\tilde{\mathbf{Y}} = (\tilde{\mathbf{y}}_n)_{n=1}^N \in \mathbb{R}^{N \times L}$ set of $N$ individual reproduced observation vectors $\tilde{\mathbf{y}}_n$. This reproduced observation follows the posterior predictive distribution Eq. 2.46, and is used for posterior predictive assessment. Introduced in Chapter 2 (Section 2.3.2). Uses in interstellar medium studies are reviewed in Chapter 3 (Section 3.3.2). Our approach is presented in Chapter 5 (Section 5.3). When $N=1$, this set is a vector and the notation is simplified to $\tilde{\mathbf{y}} \in \mathbb{R}^L$. When $N=L=1$, this set is a scalar and the notation is simplified to $\tilde{\mathbf{y}} \in \mathbb{R}$.

$D$ dimension of an individual physical parameter vector $\theta$.

$L$ dimension of an individual observation vector $\mathbf{y}$.

$N$ number of individual observations $\mathbf{y}_n$ in the full observation set $\mathbf{Y}$. Also number of individual physical parameter vector $\theta_n$ in the full parameter set $\mathbf{\Theta}$.

$\ln x$ natural logarithm of a scalar $x \in \mathbb{R}$, i.e., $\ln x = \log_e x$. This notation is preferred to the usual log to avoid any confusion between logarithms in base $e$ and 10.

$\pi(\cdot)$ probability density function (pdf). The distribution associated to a pdf is indicated by the variables in the parentheses. In this thesis, we use an abuse of notation to draw samples from a distribution: for instance, $\mathbf{\Theta}^{(1)} \sim \pi(\mathbf{\Theta})$ corresponds to a parameter $\mathbf{\Theta}^{(1)}$ drawn from the prior distribution.

$\pi(\mathbf{Y}|\mathbf{\Theta})$ likelihood function. The likelihood function is defined from an observation model $\mathcal{M}$. It is sometimes written $\pi(\mathbf{Y}|\mathbf{\Theta}, \mathcal{M})$ to highlight this dependence. In this thesis, the negative log-likelihood is assumed to be twice differentiable, i.e., $\mathbf{\Theta} \mapsto -\ln \pi(\mathbf{Y}|\mathbf{\Theta}) \in \mathcal{C}^2$.

$\pi(\mathbf{\Theta})$ pdf of the prior distribution. In this thesis, the negative log-prior is assumed to be twice differentiable, i.e., $\mathbf{\Theta} \mapsto -\ln \pi(\mathbf{\Theta}) \in \mathcal{C}^2$.

$\pi(\mathbf{\Theta}|\mathbf{Y})$ pdf of the posterior distribution, defined in Eq. 2.1. Like the likelihood function, the posterior is defined from an observation model $\mathcal{M}$. It is sometimes written $\pi(\mathbf{\Theta}|\mathbf{Y}, \mathcal{M})$ to highlight this dependence. In this thesis, the negative log-posterior is assumed to be twice differentiable, i.e., $\mathbf{\Theta} \mapsto -\ln \pi(\mathbf{\Theta}|\mathbf{Y}) \in \mathcal{C}^2$.

$-\nabla \ln \pi(\mathbf{\Theta}|\mathbf{Y}) \in \mathbb{R}^{ND}$ gradient of the negative log-pdf of the posterior distribution with respect to $\mathbf{\Theta}$.
$-\nabla^2 \ln \pi(\Theta|\mathbf{Y}) \in \mathbb{R}^{ND \times ND}$ Hessian matrix of the negative log-pdf of the posterior distribution with respect to $\Theta$.

$\mathcal{L}$ generic loss function, to be minimized. In inverse problems, it is often defined as the negative log-likelihood or the negative log-pdf of the posterior distribution. For minimization algorithms, see Chapter 2 (Section 2.2.1). For possible definitions in machine learning, see Chapter 4. In this thesis, this function is assumed to be twice differentiable, i.e., $\mathcal{L} \in \mathcal{C}^2$.

$\mathcal{M}$ observation model. Includes both the forward model $\mathbf{f}$ and the noise model $\mathcal{A}$.

$\mathbf{f}$ true forward model. It is a vector function $\mathbf{f} : \theta \in \mathbb{R}^D \mapsto \mathbf{f}(\theta) = (f_\ell(\theta))_{\ell=1}^L \in \mathbb{R}^L$. In this thesis, this function is assumed to be twice differentiable, i.e., $\mathbf{f} \in \mathcal{C}^2$.

$\tilde{\mathbf{f}}$ approximation of the true forward model, also called emulator. Like the true forward model $\mathbf{f}$, this function is assumed to be twice differentiable, i.e., $\tilde{\mathbf{f}} \in \mathcal{C}^2$. Using an emulator is common when the evaluation time of the true forward model is long. For a review of emulators in interstellar medium and cosmology, see Chapter 3 (Section 3.1.1). Defining this emulator is the core of Chapter 4.

$\psi$ parameters of the approximation of the forward model.

$\mathcal{A}$ general noise model. See Chapter 3 (Section 3.1.2) for a review of noise models in interstellar medium studies. Chapter 3 (Section 3.4) introduces the noise model considered in the class of inverse problems addressed in this thesis.

$\varepsilon^{(m)} = (\varepsilon_{nl}^{(m)}) \in \mathbb{R}^{N \times L}$ multiplicative noise.

$\varepsilon^{(a)} = (\varepsilon_{nl}^{(a)}) \in \mathbb{R}^{N \times L}$ additive noise.

$\Sigma^{(m)} \in \mathbb{R}^{NL \times NL}$ covariance matrix of multiplicative noise.

$\Sigma^{(a)} \in \mathbb{R}^{NL \times NL}$ covariance matrix additive noise.

$\omega_{nl} \in \mathbb{R}$ censorship lower bound on observations $y_{nl}$.

$\mathcal{C} \subset \mathbb{R}^D$ validity set for the physical parameter vector $\theta$. It is usually defined as a product of validity intervals for each component $\theta_d$. This set is used to define the uniform prior from Chapter 5 (Section 5.1.2).

$\tau \in \mathbb{R}^D$ regularization weight for the spatial prior defined in Chapter 5 (Section 5.1.2).

$\eta > 0$ step size, used in gradient descent algorithms, e.g., in Eq. 2.14, or Langevin sampling algorithms, e.g., in Eq. 2.36.

$\mathbf{G} \in \mathbb{R}^{ND \times ND}$ preconditioning matrix, also called preconditioner. This matrix is usually a positive definite approximation of the Hessian matrix. Defined in Eq. 2.17. Our choice of preconditioner is presented in Chapter 5 (Section 5.2.1).

$q$ proposal distribution in sampling algorithms.

$\rho \in [0, 1]$ acceptance probability of Metropolis-Hastings algorithm. Defined in Eq. 2.34.

$\tilde{\rho} \in [0, 1]$ generalized acceptance probability of multiple-try Metropolis algorithm. Defined in Eq. 2.44.

$P_{\text{th}} \ge 0$ thermal pressure, expressed in K cm$^{-3}$.

$G_0 \ge 0$ intensity of a UV radiative field. In this thesis, its is expressed in reference to the Habing field (Habing, 1968).

$A_V^{\text{tot}} \ge 0$ visual extinction, expressed in magnitudes (mag).
# Introduction

"Every adventure requires a first step."

Cheshire Cat in Lewis Carroll's "Alice's
Adventures in Wonderland"

## Contents

<table><tr><td>Scientific context</td><td>1</td></tr><tr><td>PhD project summary</td><td>2</td></tr><tr><td>Structure of the manuscript</td><td>6</td></tr><tr><td>List of publications</td><td>7</td></tr></table>

## Scientific context

This PhD project is a collaboration between statisticians and astrophysicists. It was co-supervised by Pierre Chainais, researcher in signal processing in the Centre de Recherche en Informatique, Signal et Automatique de Lille (CRISTAL), and Franck Le Petit, astronomer in the Observatory of Paris. It was also co-advised by Pierre-Antoine Thouvenin, researcher in signal processing in CRISTAL, and Emeric Bron, assistant astronomer at the Observatory of Paris. This project was funded by the Mission pour les Initiatives Transverses et Interdisciplinaires (MITI) of the Centre National de la Recherche Scientifique (CNRS), within the 80|Prime fund OrionStat.

At a larger scale, this project is part of the ORION-B consortium, led by Jérôme Pety and Maryvonne Gérin and founded in 2015. This consortium gathers experts in signal processing, machine learning, and astrophysics, both on the theoretical side and the observational side. Its members meet weekly for updates in ongoing research projects and twice a year for strategic and in-depth discussions on current and future projects. In particular, this consortium analyzes the observations of the IRAM-30m Large Program "Orion B", described in Pety et al. (2017). The observation program covered about 5 square degrees of the celestial sphere to map the Orion B giant molecular cloud (GMC). It produced a hyper-spectral image with 1 million pixels and 200 000 spectral channels in the radio frequency domain, allowing to map the emission of dozens of molecules over the whole cloud.

Figure 1 shows the resulting Orion B map in radio wavelengths along with an observation of the same region in the visible domain. Multiple statistical analyses were then led from this rich dataset, e.g., principal component analysis (PCA) on integrated line intensities (Gratier et al., 2017) and clustering of the intensities to segment the map (Bron et al., 2018b).
<figure><img src="image_2.png" /><figcaption>Figure 1 – Orion B giant molecular cloud. On each image, the moon is shown for scale.</figcaption></figure>

PhD project summary

This PhD project aims at developing statistical tools to study the physics of the interstellar medium (ISM). The ISM is a very diffuse medium that fills the extraordinarily large volume in a galaxy between celestial objects such as stars and black holes. It comprises a wide variety of environments. Most of the ISM volume is hot, ionized, and diffuse. Conversely, most of its mass is cold, neutral and dense enough for hydrogen to be in its molecular form, and constitutes relatively small regions called *molecular clouds*.

**The study of the ISM** carries fundamental questions such as star formation or the development of molecular complexity possibly leading to the formation of prebiotic molecules. This thesis focuses on star formation and their feedback on molecular clouds. Stars are born from the gravitational collapse of a part of these clouds. Newborn stars impact their parent cloud with their ultraviolet (UV) irradiation and stellar winds, as well as through supernovae explosions at the end of their lives (for the most massive ones). The overall impact of these feedback processes on the remains of the parent cloud is to this day only partially understood. The feedback of newborn stars might dissipate their parent cloud, which would prevent the formation of other stars. It might as well locally compress parts of the surrounding cloud, and thus favor the formation of other stars.

**Observations and physical parameters** – The observations considered in this work are hyperspectral maps of molecular clouds in the far infrared and millimeter wavelength domains. We focus on clouds that are illuminated and heated by nearby massive stars emitting UV photons. The surface layer of such clouds, where the UV irradiation heats and dissociates the molecular gas, is called a photodissociation region (PDR). The ions, atoms and molecules present in the cloud cool mostly through observable radiative emission associated with quantum transitions. These hyperspectral maps are reduced to multispectral maps that we denote $Y = (y_n) \in \mathbb{R}^{N \times L}$, where each pixel $y_n$ contains the integrated intensities of $L \sim 5-30$ emission lines. We work on observation maps that contain from $N=1$ to $N = O(10^4)$ pixels. These maps are therefore much smaller than that of the Orion B cloud which contains other types of environments than PDRs such as dense cores.

These multispectral maps **Y** can be compared with predictions of a numerical model of the ISM such as the Meudon PDR code (Le Petit et al., 2006), that we denote **f**. Such models can produce predictions of the line intensity maps for any provided map of the local physical parameters, denoted $\Theta = (\theta_n) \in \mathbb{R}^{N \times D}$, where each vector $\theta_n$ contains $D \lesssim 10$ physical parameters such as the gas density, the thermal pressure, and the total thickness of the cloud along the
line of sight. The parameters $\Theta$ thus live in high dimension, but show a simple pixel structure. Estimating a map of physical parameters $\Theta$ from an observation map $Y$ and the Meudon PDR code $f$ is an instance of a general class of problems called inverse problems. Figure 2 illustrates the principle of an inference procedure.

<figure><img src="image_3.png" /><figcaption>Figure 2 – General principle of an inference procedure, repeated until a stopping criterion is satisfied.</figcaption></figure>

**Optimization-based inference and existence of multiple solutions** – In astrophysics, the inverse problem is often formulated as an optimization problem with a loss function $\mathcal{L}$, which measures the difference between the model prediction $f(\Theta)$ and the observation $Y$. The estimated physical parameter $\hat{\Theta}$ is then the minimum of a loss function,

$$ \hat{\Theta} \in \arg\min_{\Theta} \mathcal{L}(\Theta; Y), \qquad (1) $$

where $\arg\min$ is the set of values of $\Theta$ that minimize the loss function. For instance, the squared loss

$$ \mathcal{L}(\Theta; Y) = \|f(\Theta) - Y\|_2^2 \qquad (2) $$

is a widespread loss function. It turns out that this inverse problem combines many challenges. One example of difficulty is the potential existence of multiple solutions.

Figure 3 shows two cases where multiple solutions reconstruct equally well the observations, i.e., where the loss function has multiple local minima. The signal-to-noise ratio (SNR) can cover multiple decades in ISM observations. In low SNR regions, multiple physical parameter values $\Theta$ fit the observations equally well. Figure 3a shows that such solution degeneracy is problematic for an optimization problem, as all the values in the wide valley are equally valid. Returning only one value hides this degeneracy. In high intensity regions, the SNR may be high enough for accurate estimations. However, ISM numerical models yield non-linear relations between the physical parameters and the observables so that the cost function is generally non-convex and may thus contain multiple local minima. Figure 3b shows a cost function with two equivalent local minima. Returning only one value hides the existence of the other local minimum.

**Uncertainty quantification with Bayesian sampling-based inference** – A Bayesian sampling approach quantifies the uncertainty associated with an inference, and can identify these two cases – degeneracy and multiple local minima. In such an approach, the unique estimator $\hat{\Theta}$ in the optimization problem from Eq. 1 is replaced by a random variable $\Theta$. The cost function $\mathcal{L}(\Theta; Y)$ is then replaced by a probability distribution $\pi(\Theta|Y)$, called the *posterior distribution*.

Figure 4 shows the two cost functions of Figure 3 replaced by posterior probability distributions. The posterior distribution contains the information on the physical parameters $\Theta$ provided an observation $Y$. In particular, it gives access to estimators, to high probability regions and credibility intervals, and permits to identify local minima, all at once.
<figure><img src="image_4.png" /><figcaption>Figure 3 – Illustration of two types of solution multiplicity in an optimization problem. In both cases, an optimization procedure would return either $\hat{\Theta}^{(1)}$ or $\hat{\Theta}^{(2)}$, but would not detect the existence of the other.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4 – Replacing a loss function by a posterior distribution. The gray areas represent high probability regions that define credibility intervals. They can be obtained with a Bayesian sampling approach.</figcaption></figure>
Extracting information from the posterior distribution often involves integrals over the physical parameters $\Theta$, e.g., for the posterior expectation

$$ \mathbb{E}[\Theta|\mathbf{Y}] = \int \Theta \pi(\Theta|\mathbf{Y}) d\Theta. \qquad (3) $$

In general, such integrals are intractable and require numerical evaluation. Using Riemannian integration is possible when the dimension of $\Theta$ is low, i.e., when $ND \lesssim 10$. However, this method does not scale well to high dimensions as the number of necessary points grows exponentially with $ND$. In the considered astrophysical inverse problems, the number of pixels $N$ ranges from $O(1)$ to $O(10^4)$, and the number of parameters per pixel ranges from 1 to 10. Integrals are thus approximated using Monte Carlo estimators (Robert and Casella, 2004, chapter 3), with $T_{MC}$ samples of the distribution $\Theta^{(t)} \sim \pi(\Theta|\mathbf{Y})$. For instance, the expectation of the posterior distribution in Eq. 3, can be estimated with the empirical mean of samples, also called minimum mean square error (MMSE),

$$ \hat{\Theta}_{\text{MMSE}} = \frac{1}{T_{\text{MC}}} \sum_{t=1}^{T_{\text{MC}}} \Theta^{(t)} \approx \mathbb{E}[\Theta|\mathbf{Y}]. \qquad (4) $$

A general algorithm that efficiently draws independent and identically distributed (i.i.d.) samples from any probability distribution does not exist. For a large set of distributions, Markov chain Monte Carlo (MCMC) algorithms (Robert and Casella, 2004, chapters 6 and 7) generate correlated chains of samples $\Theta^{(t)}$. Such algorithms are often called samplers.

**Tasks addressed in this thesis** – The aforementioned inverse problem is addressed within a Bayesian framework to address the absence of ground truth inherent to astrophysics. Uncertainty quantification, e.g., with credibility intervals, is derived along with point estimates. Solving this inverse problem involved a few main tasks. In each task, the uncertainty associated with proposed approximations is controlled.

- The Meudon PDR code requires a few hours per evaluation, which is prohibitively slow for inference that relies on many evaluations of the likelihood function. It is assumed to be a non-linear and twice differentiable function. Besides, it is a strictly positive function that covers multiple decades, which makes it non gradient Lipschitz continuous or with an extremely large Lipschitz constant. It will be approximated by a light, fast and accurate neural network-based emulator in Chapter 4.

- The full observation model described in Chapter 3 is complex as it involves two sources of noise – one additive and Gaussian, the other multiplicative and lognormal –, and censorship. Without neglecting any source, it leads to a likelihood function whose expression is challenging to address as is. The likelihood function will be approximated with a small and controlled error in Chapter 5.

- The signal-to-noise ratio (SNR) in observation maps greatly varies across lines and pixels. When the SNR is low, the physical parameters usually are poorly constrained. A spatial regularization prior will be included to improve the quality of estimations in these regions. It will enable low SNR pixels to access the information contained in the neighboring pixels.

- The resulting posterior distribution is difficult to sample from. In particular, it is non-log-concave and potentially multimodal because of the non-linearity of the forward PDR model. Besides, the log-posterior is twice differentiable but non gradient Lipschitz continuous. A dedicated MCMC algorithm will be proposed in Chapter 5 as a combination of two sampling kernels: one will efficiently explore the posterior distribution locally, the other will allow escaping from local minima.

- Due to the complexity of the ISM, the Meudon PDR code includes multiple microphysical processes and relies on simplifying assumptions. Its compatibility with the observation will be assessed thanks to a Bayesian hypothesis testing approach presented in Chapter 5.
— The resulting overall method will be applied to synthetic observations in Chapter 5 and to multiple real observations in Chapter 6.

An additional task was also addressed, that is the determination of the best emission lines to observe to obtain low uncertainties on physical parameter estimates. However, this work is still in progress. We currently work on it in collaboration with other members of the ORION-B consortium. Preliminary results are presented in Appendix A.

Addressing these tasks led to contributions in both the statistical and ISM communities.

## Structure of the manuscript and reader's guide

This manuscript targets both the signal processing and astrophysics communities. Though the backgrounds of these two communities share some common references, they strongly differ in mindset and vocabulary. Key notions of both statistics and astrophysics are thus re-introduced and discussed. The “Notation” section gathers the main mathematical objects used throughout this thesis, along with the main hypotheses and properties. Along the manuscript, some boxes separated from the main text target one specific community to provide examples or a detailed explanation. In addition, sections whose title begins with “Illustration: ...” provide illustrative examples. These examples are detailed walkthroughs, applications on simplified cases or effect illustrations, but do not contain essential content. The reader may very well choose not to read these sections. Finally, this thesis is structured such that it can be read in two ways.

**Option 1: full reading with the default structure** – This manuscript is divided into six chapters and one appendix.

**Chapter 1** provides some background on the ISM and the associated scientific questions. Then, it focuses on photodissociation regions (PDRs), a particular type of ISM environment. The model used throughout this work to simulate these regions, the Meudon PDR code, is also presented.

**Chapter 2** provides some background in statistics. Necessary notions for modeling are introduced. Inference methods from optimization and Bayesian approaches, as well as existing alternatives, are then reviewed. Some Markov chain Monte Carlo (MCMC) algorithms that are necessary to introduce the proposed sampler of Chapter 5 are described.

**Chapter 3** is a state-of-the-art that reviews applications of statistical models and inference algorithms to ISM studies.

**Chapter 4** describes a model reduction step that is necessary for the proposed sampler of Chapter 5. An artificial neural network (ANN) is proposed to emulate the Meudon PDR code. Some specificities of the Meudon PDR code make traditional off-the-shelf ANN tools inefficient. We design and train ANNs that address these specificities. The corresponding paper, Palud et al. (2023c), was published in the international journal Astronomy & Astrophysics (A&A).

**Chapter 5** details the statistical methods we developed to address the inverse problem and evaluate the results. The proposed model and sampler are described. This sampler was presented in a journal article published in the international journal Transactions on Signal Processing (TSP), Palud et al. (2023b). It was also described in two conference papers, one at the French Gretsi national conference (Palud et al., 2022a) and one at the international EUSIPCO conference (Palud et al., 2022b). A proposed model checking strategy enables validating the results of an inversion. This model checking approach is discussed in a conference Gretsi paper, Palud et al. (2023a).
**In Chapter 6**, the MCMC algorithm presented in Chapter 5 and the reduced model of Chapter 4 are applied to real observations. Inversion results are presented and analyzed. A journal article is currently in preparation, Palud et al. (in prep[a]).

**Appendix A** proposes a general method to determine which lines are most informative for inference using the differential entropy of a probability distribution. This method paves the way for a new variable selection method. An associated article is currently in preparation, Palud et al. (in prep[b]).

**Option 2: reading according to topics** – A second possibility is to read this document by topic. Table 1 lists the sections associated with each topic.

Table 1 – Reading options per topic. “full” means that a given topic is covered in the whole chapter.

<table><thead><tr><th>Topic</th><th>Ch. 1<br>ISM</th><th>Ch. 2<br>statistics</th><th>Ch. 3<br>interface<br>ISM and<br>statistics</th><th>Ch. 4<br>ANNs</th><th>Ch. 5<br>MCMC<br>algorithm</th><th>Ch. 6<br>real<br>data</th><th>App. A<br>line<br>selection</th></tr></thead><tbody><tr><td><strong>ISM</strong></td><td>full</td><td>–</td><td>full</td><td>–</td><td>5.4.2</td><td>full</td><td>full</td></tr><tr><td><strong>Statistical model</strong></td><td>–</td><td>2.1</td><td>3.1, 3.4</td><td>full</td><td>5.1</td><td>full</td><td>full</td></tr><tr><td>forward model</td><td>1.2, 1.3</td><td>–</td><td>3.1.1</td><td>full</td><td>–</td><td>full</td><td>full</td></tr><tr><td>noise model</td><td>–</td><td>–</td><td>3.1.2, 3.4</td><td>–</td><td>–</td><td>full</td><td>full</td></tr><tr><td>prior</td><td>–</td><td>–</td><td>3.1.3</td><td>–</td><td>5.1.2</td><td>–</td><td>–</td></tr><tr><td><strong>Inference</strong></td><td>–</td><td>2.2, 2.A</td><td>3.2</td><td>–</td><td>5.2</td><td>full</td><td>full</td></tr><tr><td>local exploration</td><td>–</td><td>2.2,</td><td>3.2</td><td>–</td><td>5.2.1</td><td>–</td><td>–</td></tr><tr><td>escaping local modes</td><td>–</td><td>2.2, 2.A</td><td>3.2</td><td>–</td><td>5.2.2</td><td>–</td><td>–</td></tr><tr><td><strong>Model relevance</strong></td><td>–</td><td>2.3</td><td>3.3</td><td>–</td><td>5.3</td><td>full</td><td>–</td></tr><tr><td>model selection</td><td>–</td><td>2.3.1</td><td>3.3.1</td><td>–</td><td>–</td><td>–</td><td>–</td></tr><tr><td>model checking</td><td>–</td><td>2.3.2</td><td>3.3.2</td><td>–</td><td>5.3</td><td>full</td><td>–</td></tr></tbody></table>

# List of publications

## International journals

Palud, P., P.-A. Thouvenin, P. Chainais, E. Bron, and F. Le Petit (2023b). « Efficient sampling of non log-concave posterior distributions with mixture of noises ». *IEEE Transactions on Signal Processing* 71, pp. 2491–2501.

Palud, P. et al. (2023c). « Neural network-based emulation of interstellar medium models ». *A&A* 678.A198.

## International conferences

Palud, P. et al. (2022b). « Mixture of noises and sampling of non-log-concave posterior distributions ». *2022 30th European Signal Processing Conference (EUSIPCO)*, pp. 2031–2035.

## National conferences

Palud, P., P. Chainais, F. Le Petit, P.-A. Thouvenin, and E. Bron (2023a). « Problèmes inverses et test bayésien d'adéquation du modèle ». *29° Colloque sur le traitement du signal et des images*. Grenoble: GRETSI - Groupe de Recherche en Traitement du Signal et des Images, p. 697–700.
Palud, P. et al. (2022a). « Mélange de bruits et échantillonnage de posterior non log-concave ». 28° *Colloque sur le traitement du signal et des images*. 001-0176. Nancy: GRETSI - Groupe de Recherche en Traitement du Signal et des Images, p. 705-708.

**International journal articles in preparation**

Palud, P., P.-A. Thouvenin, P. Chainais, E. Bron, and F. Le Petit (in prep[a]). « Bayesian inversion of large interstellar medium observation maps ».

Palud, Pierre, Pierre-Antoine Thouvenin, Pierre Chainais, Emeric Bron, and Franck Le Petit (in prep[b]). « Entropy-based selection of most informative observables for inference from interstellar medium observations ».

The above published conference articles were published in statistical conferences. I also participated to astrophysical international conferences:

- the IRAM<fnref n="1" /> 2022 conference (https://iram2022nice.sciencesconf.org/),

- the PCMI<fnref n="2" /> 2022 conference (https://pcmi2022paris.sciencesconf.org/),

- the EAS<fnref n="3" /> 2023 conference (https://eas.unige.ch/EAS_meeting/index.jsp).

In the ISM community, conferences are based on abstracts and not on articles.

**Software**

The software developed during this thesis was published in two PYTHON packages.

The package associated with the neural network approximation presented in Chapter 4 is called NNBMA (for Neural Network-Based Model Approximation). The code can be found at https://github.com/einigl/ism-model-nn-approximation. It can be installed via PIP – see https://pypi.org/project/nnbma/. The associated documentation can be found at https://ism-model-nn-approximation.readthedocs.io/en/latest/? badge=latest.

The package associated with the Bayesian inference algorithm presented in Chapter 5 is called BEETROOTS (for BayEsian infErence with spaTial Regularization of nOisy multi-line ObservaTion 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 documentation can be found at https://beetroots.readthedocs.io/en/latest/index.html.

<footnote n="1">1. Institut de radioastronomie millimétrique</footnote>

<footnote n="2">2. Physique et chimie du milieu interstellaire</footnote>

<footnote n="3">3. European astronomical society</footnote>
# Part I

**Backgrounds: interstellar medium, statistical inference and interactions**
<figure><img src="image_6.png" /></figure>
# Chapter 1

# Background on the interstellar medium
(ISM)

> " The Milky Way is largely empty. Stars are separated by some 2 pc in the solar neighborhood. [...] If we take our Solar System as a measure, with a heliosphere radius of ≈235 AU, stars and their associated planetary systems fill about 3 × 10⁻⁸ % of the available space. This [work] deals with what is in between these stars: the interstellar medium (ISM). "

Tielens (2005, introduction of chapter 1)

## Contents

<table><tr><td>1.1</td><td>The interstellar medium (ISM)</td><td>13</td></tr><tr><td>1.1.1</td><td>A short overview of the ISM</td><td>13</td></tr><tr><td>1.1.2</td><td>Star formation, stellar feedback and photodissociation regions (PDRs)</td><td>15</td></tr><tr><td>1.1.3</td><td>Chemical complexity and planet formation</td><td>17</td></tr><tr><td>1.1.4</td><td>PDRs: structure and physical parameters</td><td>18</td></tr><tr><td>1.2</td><td>Linking physical conditions and observations: astrophysical numerical models</td><td>20</td></tr><tr><td>1.2.1</td><td>Radiative transfer models</td><td>20</td></tr><tr><td>1.2.2</td><td>Astrochemical models</td><td>21</td></tr><tr><td>1.2.3</td><td>Dust models</td><td>22</td></tr><tr><td>1.2.4</td><td>Holistic models</td><td>22</td></tr><tr><td>1.3</td><td>The Meudon PDR code</td><td>23</td></tr><tr><td>1.3.1</td><td>Physics and chemistry taken into account</td><td>23</td></tr><tr><td>1.3.2</td><td>Input parameters and considered environments</td><td>25</td></tr><tr><td>1.3.3</td><td>Limitations</td><td>26</td></tr><tr><td>1.4</td><td>Conclusion</td><td>26</td></tr></table>

This chapter provides a general overview of the ISM and of the associated scientific questions. It briefly introduces the physics of the ISM and the current state-of-the-art numerical models. The elements and notions covered in this chapter will be used throughout this thesis. Two fundamental notions necessary for the non astrophysicist are first introduced: distance units and scales.
**Distance units** – The solar system that we live in hosts one of the ~10<sup>11</sup> stars contained in our Galaxy, the Milky Way (MW). Then again, the Milky Way is one of the ~10<sup>11</sup> galaxies of the observable Universe. The study of the Universe at its largest scales falls into the domain of cosmology. In this work, we consider “smaller” astrophysical scales, between the scales of stars and of galaxies. Dedicated distance units are necessary to describe such scales. Three distance units are commonly used in astrophysics:

— The astronomical unit (au): 1 au corresponds to the distance between the Sun and the Earth, i.e., 1.5 × 10<sup>6</sup> km.

— The light-year (ly): 1 ly corresponds to the distance traveled in one year at the speed of light, i.e., 9.5 × 10<sup>12</sup> km.

— The parsec (pc): 1 pc is defined as <sup>648 000</sup>&frasl;<sub>&pi;</sub> au, i.e., about 3.26 ly. Figure 1.1 illustrates its definition. One parsec corresponds to the order of magnitude of the average distance between a star and the closest other star.

<figure><img src="image_7.png" /><figcaption>Figure 1.1 – Illustration of the parsec (pc) definition. Distances are not at scale.</figcaption></figure>

**Scales** – Figure 1.2 lists some typical astrophysical distances with the corresponding values in the three aforementioned units. It covers more than 16 decades. This work focuses on molecular clouds such as Orion-B. Molecular clouds correspond to one type of environment of the so-called ISM, which by definition lies between stars. They are extremely large: their sizes range from a few parsecs to a hundred parsecs. For instance, Orion-B is a “neighbor” molecular cloud, about 400 pc away from us. It is 10 pc large, which is about 10<sup>10</sup> times larger than the diameter of the Earth. This size ratio corresponds to that of a person and a hydrogen atom.

In this chapter, we provide a general description of the ISM. Section 1.1 provides an overview of the ISM, of the associated physics and scientific questions. Section 1.2 describes how models of the ISM and observations can improve our understanding of star formation. We also list some existing ISM models. Finally, Section 1.3 focuses on the Meudon PDR code, which is used in the remainder of this work.
<figure><img src="image_8.png" /><figcaption>Figure 1.2 – Scales in astrophysics. α-Centauri is the closest star to the Sun. The Andromeda galaxy is the closest galaxy to the Milky Way. Orion B is a giant molecular cloud (GMC) of the Milky Way.</figcaption></figure>

# 1.1 The interstellar medium (ISM)

In the Milky Way, stars and other celestial bodies such as black holes, planets, and comets represent about 90% of the total baryonic mass, i.e., disregarding dark matter and dark energy. The remaining 10% correspond to the extremely diffuse medium that lies in the volume between stars: the ISM. The study of the ISM addresses fundamental questions on the Universe including the formation of stars and planets, and the development of molecular complexity possibly leading to the formation of prebiotic molecules. In this section, we provide a general overview of the ISM with its composition, structure, and role in the cycle of matter in a galaxy. We also detail the aforementioned fundamental questions.

## 1.1.1 A short overview of the ISM

**Composition** – Hydrogen remaining from the Big Bang represents 90% in mass of the ISM. Helium corresponds to most of the remaining 10%. It comes in part from the Big Bang and in part from nuclear reactions inside stars. Other elements, from carbon to uranium, are also produced by nuclear reactions inside stars, and by supernovae for the most massive elements. They account for less than 1% of the total mass. These other elements, in particular carbon C and oxygen O, play a key role in the ISM physics and observations. The concentration of elements heavier than helium is called the *metallicity*. Table 1.1 lists estimated abundances of the most common elements in the ISM relatively to hydrogen in the solar neighborhood. These elements can exist in ionic, atomic, or molecular form, depending on the environment conditions.

In addition to the gas, the ISM contains small dust grains whose size typically ranges from 1 to 100 nm, and larger in the ISM densest regions and in protoplanetary disks. It is estimated that these dust grains represent about 1% of the ISM mass. They play a key role in the physics and chemistry of the ISM. They extinguish UV radiation and thus shield the gas from its ionizing and dissociating effects. Dust grains also convert a small fraction of the absorbed UV energy into gas heating through the photoelectric effect, which is one of the dominant heating processes in the ISM. Some species exist in both gas and dust phase, and are thus depleted in the gas phase – a phenomenon called the *interstellar depletion* (Draine, 2011, chapter 9). Heavy elements such as
<figure><img src="image_9.png" /><figcaption>(a) 2014 Hubble observation, in visible light.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) 2022 JWST observation, in NIR. Credits: NASA, ESA, CSA, STScI; Joseph DePasquale (STScI), Anton M. Koekemoer (STScI), Alyssa Pagan (STScI).</figcaption></figure>

Figure 1.3 – observation of the pillars of creation, in the Eagle Nebula, by the Hubble telescope and the James Webb spatial telescope (JWST). Unlike the observation in the visible light, the near infrared (NIR) observation shows the stars inside and behind the pillars.

silicon Si and iron Fe are extremely depleted, as their gas phase abundances correspond to only a few percent of the total expected abundance. These elements, along with lighter elements such as carbon C and oxygen O, are part of the components of the solid cores of grains. In addition, in the dense cold regions, water and CO molecules from the gas phase freeze out on the dust surface to form ice mantles. The depletion of C and O and other thus increases strongly in the dense cold regions where ice mantles can form. Many chemical reactions then occur at the surface of grains and within ice mantles.

<table><caption>Table 1.1 – Solar abundance of elements in the gas phase of the interstellar medium relative to hydrogen. Adapted from Draine (2011, table 1.4).</caption><thead><tr><th>Element</th><th>Abundance (relative to H)</th></tr></thead><tbody><tr><td>He</td><td>9.55 × 10<sup>-2</sup> ±2.3%</td></tr><tr><td>C</td><td>2.95 × 10<sup>-4</sup> ±12.2%</td></tr><tr><td>N</td><td>7.41 × 10<sup>-5</sup> ±12.2%</td></tr><tr><td>O</td><td>5.37 × 10<sup>-4</sup> ±12.2%</td></tr><tr><td>Si</td><td>3.55 × 10<sup>-5</sup> ±9.6%</td></tr><tr><td>S</td><td>1.45 × 10<sup>-5</sup> ±7.2%</td></tr><tr><td>Fe</td><td>3.47 × 10<sup>-5</sup> ±9.6%</td></tr></tbody></table>

**Phases in the ISM** – The ISM is a very inhomogeneous medium. It goes from very diffuse hot ionized environments with a density of about $3 \times 10^{-3}$ particles cm$^{-3}$ that represent most of the volume of the ISM, to dense molecular clouds with $10^6$ particles cm$^{-3}$ that represent most of its mass. For comparison, at a standard pressure, the Earth atmosphere contains about $10^{19}$ particles cm$^{-3}$ and the vacuum at the surface of the Moon about $10^5$ particles cm$^{-3}$ (Öpik, 1962).

Table 1.2, adapted from Galliano (2022), lists the main phases of the ISM along with their respective volume densities, average temperature and volume filling factor. The main distinctions between these phases are the state of hydrogen – ionized H$^{+}$, atomic H or molecular H$_2$ –, their temperature and their volume density, i.e., the average number of particles per cubic centimeter. The hot ionized medium (HIM), warm ionized medium (WIM), and HII regions correspond to hot plasmas. In particular, HII regions are located around a massive star where the gas is
ionized by the extreme ultraviolet (EUV) photons of the star. The EUV field contains photons with wavelengths smaller than 91.2 nm, i.e., 13.6 eV, the hydrogen ionization potential. This thesis focuses on the remaining four phases, where hydrogen is mostly neutral and where stars are formed. Hydrogen is mostly in its atomic form in warm neutral medium (WNM) and cold neutral medium (CNM), and in its molecular form in diffuse H₂ and dense H₂ regions. These regions form molecular clouds such as Orion-B, shown in Figure 1, or the pillars of creation, shown in Figure 1.3.

<table><caption>Table 1.2 – Comparison of the main phases of the ISM. The sum of volume filling factors exceeds 100% because the associated uncertainties are large. Adapted from Galliano (2022).</caption><thead><tr><th>Hydrogen state</th><th>Phase</th><th>Volume density n<br>(in cm<sup>-3</sup>)</th><th>Temperature<br>(in K)</th><th>Volume filling factor<br>(in %)</th></tr></thead><tbody><tr><td rowspan="3">ionized</td><td>HIM</td><td>≈ 3 × 10<sup>-3</sup></td><td>≈ 10<sup>6</sup></td><td>≈ 50</td></tr><tr><td>WIM</td><td>≈ 0.1</td><td>≈ 10<sup>4</sup></td><td>≈ 25</td></tr><tr><td>HII regions</td><td>≈ 1 – 10<sup>5</sup></td><td>≈ 10<sup>4</sup></td><td>≈ 1</td></tr><tr><td rowspan="2">atomic</td><td>WNM</td><td>≈ 0.3</td><td>≈ 10<sup>4</sup></td><td>≈ 30</td></tr><tr><td>CNM</td><td>≈ 30</td><td>≈ 100</td><td>≈ 1</td></tr><tr><td rowspan="2">molecular</td><td>Diffuse H<sub>2</sub></td><td>≈ 100</td><td>≈ 50</td><td>≈ 0.1</td></tr><tr><td>Dense H<sub>2</sub></td><td>≈ 10<sup>3</sup> – 10<sup>6</sup></td><td>≈ 10</td><td>≈ 0.01</td></tr></tbody></table>

**The cycle of matter in a galaxy** – The ISM plays a key role in the evolution of matter in a galaxy. Figure 1.4 shows the cycle of baryonic matter. The ISM provides the building material for star formation. In average in the Milky Way, about 1.3 $M_\odot$ yr$^{-1}$ of the ISM is converted into stars, with $1M_\odot = 2 \times 10^{30}$ kg corresponding to one solar mass. This process is detailed in the next section. This figure also shows that matter can go back from stars to the ISM in diverse ways, closing the cycle of matter in a galaxy. For instance, during its life, a star emits stellar winds that continuously feed the ISM in matter and in mechanical energy. Some stars also affect the ISM with violent deaths such as supernovae, feeding it both in energy and in elements heavier than helium.

<figure><img src="image_11.png" /><figcaption>Figure 1.4 – Observed flow of baryons in the Milky Way. Adapted from Draine (2011, chapter 1).</figcaption></figure>

### 1.1.2 Star formation, stellar feedback and photodissociation regions (PDRs)

Star formation starts with the gravitational collapse of a cold and dense region of a molecular cloud. When this region collapses, the associated gas and dust form an infalling envelope and a rotating disk. Planetary systems form from this gas and dust. Figure 1.5 illustrates the star formation process. The first three steps typically last a few million years. During this dramatic compression process, the gas collapses from scales of ~1 pc to ~10$^{-7}$ pc, which corresponds to a density increase of a factor ~10$^{21}$ (Draine, 2011, chapter 41).
<figure><img src="image_12.png" /><figcaption>Figure 1.5 – Key steps of the star formation process. Artist impression from Bill Saxton (NRAO/AUI/NSF), taken from Dishoeck and Bergin (2021).</figcaption></figure>

The gravitational collapse starts when the self-gravity of a region of the cloud exceeds a cer-
tain threshold. Different definitions of this threshold such as the Jeans instability rely on different
assumptions and geometries (Draine, 2011, chapter 41). Multiple physical phenomena resist the
cloud self gravity such as the internal thermal pressure, the magnetic pressure, turbulence at the
lowest scales, and the region angular moment. In this adversarial context, two main mechanisms
contribute to attain the density threshold that triggers the collapse. First, large scale turbulence
or other compressive events, such as Supernovae shock waves, or HII region expansion, can locally
compress the cloud. Second, atomic and molecular radiative emission evacuates energy, and thus
cools the cloud and helps reduce the thermal pressure of the cloud, allowing self-gravity to take
over.

Despite good observational evidence and advanced models, the star formation process is to this day only partially understood. Observations lead to a star formation rate (SFR) averaged over the past 3 Myr in the Milky Way of about $1.3 M_{\odot} yr^{-1}$ – as shown in Figure 1.4. A simple theoretical estimate based on free-falling gas and using typical values of the Milky Way – total molecular gas $\sim 10^9 M_{\odot}$ with density $\sim 50 \text{ cm}^{-3}$ – leads to a maximum SFR of about $200 M_{\odot} yr^{-1}$, i.e., two orders of magnitude larger than the measured value (Draine, 2011, chapter 42). Therefore, the star formation process is highly inefficient.

A first explanation of this inefficiency is that only the densest regions of the ISM are expected
to actually collapse. This hypothesis is confirmed by the Kennicutt-Schmidt empirical law, which
relates the SFR to the surface density of gas with a power law with exponent around 1.4 (Draine,
2011, chapter 42). The exponent of this power law is yet to be physically explained.

A second reason for this inefficiency is the stellar feedback that disrupts the parent cloud and thus decreases its star forming capability. We call stellar feedback all the processes by which a star injects energy into its parent cloud. It includes mechanical energy, such as stellar winds and supernovae, and radiative energy emitted in the UV domain. However, to this day, it is unclear whether stellar feedback favors or impedes star formation. Overall, stellar feedback tends to dissipate the parent cloud. For instance, the expansion of the HII bubble ionized by massive stars tends to disperse the surrounding molecular cloud. However, this expansion also causes the apparition of a compressed layer at the frontier between the molecular parent cloud and the HII region. This compressed layer at the edge of the molecular cloud and of the HII region is called a photodissociation region (PDR). This local compression of the gas might favor the formation of additional stars. Overall, the global parent cloud disruption and the local compression at the PDR are competing effects. Understanding the influence of stellar feedback requires studying the
structure of PDRs and the local physical conditions such as the thermal pressure, the volume
density, and the cloud depth along the line of sight.

**1.1.3 Chemical complexity and planet formation**

In 2018, over 200 molecules had been detected in the ISM (McGuire, 2018), organic or not, including simple diatomic molecules such as H₂ and CO, more complex molecules such as water H₂O and methanol CH₃OH, and complex ions such as H₃O⁺ and HCO⁺. Molecules from the ISM can be observed in the optical/UV domain for electronic transitions, in the IR domain for vibrational transitions, and in the millimeter domain for rotational transitions. Dust grains play a key role in the formation of some complex molecules, acting as a catalyst. For instance, methanol mainly forms on the surface of grains in CO-rich ice. Water is more complex. It is the product of a variety of reactions, some occurring in the gas phase and some taking place on the surface of dust grains, like H₂ + O. It can freeze onto the grain, forming water ice mantles. This process occurs in the deep, cold, and dense regions of a molecular cloud. These ice mantles can serve as reservoirs for water and other simple molecules and allow for a rich chemistry that would be extremely slow in the gas phase in dense and cold conditions.

As already mentioned, stars and the associated planets form from the molecular gas and dust
that gravitationally collapses. A significant part of the chemical composition of the planet-building
material appears to be set in the cold pre- and protostellar stages and preserved to planet and
comet construction sites. In the early phases of a star formation, the gas and dust spirals in a
protoplanetary disk. Figure 1.6 shows two spatially well resolved protoplanetary disks observed by
the ALMA telescope in the FIR and millimeter domains.

<figure><img src="image_13.png" /><figcaption>(a) Dust continuum observation at 1.3 mm of the disk surrounding the star HL Tau. Adapted from Andrews et al. (2016). (b) Observations of the 870 µm continuum emission from the TW Hya disk. Taken from ALMA Partnership et al. (2015).</figcaption></figure>

<caption>Figure 1.6 – ALMA observations of two protoplanetary disks.</caption>

Protoplanetary disks models study stellar system formation from the initial molecular cloud
collapse. Planets form with the aggregation of gas and dust grains. One of the key questions such
models try to address is the origin of the chemical composition of planets, in particular for the
necessary bricks for the apparition of life, such as water, simple sugars and peptide bonds. This
composition might be inherited from the parent molecular cloud and the original dust grains and
ice mantles, or might be reset with additional chemistry in the disk, either fully or partially. The
study of comets and meteorites in the solar system indicates both a full reset of the chemistry
and inheritance from the ISM (Dishoeck and Bergin, 2021). Therefore, a significant fraction of
the water on Earth might come from the ISM.
### 1.1.4 PDRs: structure and physical parameters

Modeling and observing molecular clouds permits to understand the mechanisms of star formation and its feedback, and, within protoplanetary disks, of the formation of complex molecules. In this thesis, we focus on the process of star formation feedback in molecular clouds. More precisely, we focus on the surface layer of molecular clouds irradiated by stellar UV photons, i.e. photodissociation regions (PDRs). This section introduces the key parameters that describe the physical conditions of PDRs and thus quantify the impact of stellar feedback.

Figure 1.7 illustrates the structure of a plane-parallel PDR along with the associated main physical parameters. It also shows a PDR in the observation of the Carina nebula from the JWST. The gas around the star is impacted by its UV radiation. The ISM close to the star is ionized by the EUV part of the spectrum of a star. Such regions are called HII regions (in blue in the JWST observation). Once the ionizing photons are absorbed, the hydrogen can exist in neutral atomic form. Deeper in the cloud, hydrogen forms molecules. A PDR corresponds to the region with neutral hydrogen (red regions in the JWST observation, and shaded regions in the diagram). The physics and chemistry of a PDR are controlled by far ultraviolet (FUV) dissociating photons – with energy lower than the ionization potential of H, 13.6 eV, i.e., with wavelengths greater than 91.2 nm.

<figure><img src="image_14.png" /><figcaption>Figure 1.7 – Structure of a PDR region. The image on top is extracted from the JWST Carina nebula observation. Credits to NASA, ESA, CSA, and STScI.</figcaption></figure>

**Temperature, volume density and thermal pressure** – The physical state of the gas at any point in a PDR is usually described by its temperature $T_{gas}$, expressed in K, and its proton volume density $n_{H} = n(H) + 2n(H_{2}) + n(H^{+})$, expressed in protons cm$^{-3}$. Note that the proton volume density is not equal to the volume density of particles in the gas, which is $n_{tot} = n(H) + n(H_{2}) + n(H^{+}) + n(C) + n(O) + n(CO) + ...$, expressed in particle cm$^{-3}$. The gas temperature greatly varies in a PDR. In particular, it drops at the photodissociation front
from ~ $10^2 - 10^3$ K to ~ 10 K. The proton volume density typically ranges from $10^2$ to $10^8$ cm$^{-3}$. The particle volume density and the temperature determine the thermal pressure of the gas $P_{th} = n_{tot} \times T_{gas}$, expressed in K cm$^{-3}$. In a PDR, the thermal pressure ranges from $10^5$ to $10^9$ K cm$^{-3}$.

**Amount of matter: optical depth, total visual extinction and column density – Astrophysical observations are based on emission or absorption of photons in specific directions. However, one cannot locate photon sources along the line of sight. A natural quantity to describe the ISM – including PDRs – is thus the amount of matter along the line of sight.**

## Details for statisticians 1.1: The radiative transfer equation and optical depth

The radiative transfer equation relies on the specific intensity $I_ν(x, k, t)$, expressed in erg cm$^{-2}$s$^{-1}$Hz$^{-1}$sr$^{-1}$. $I_ν(x, k, t)$ quantifies the intensity of a radiation field at frequency $ν$, for a position $x$, a direction of propagation $k$ and an instant $t$. Light and matter can interact in three ways: a photon can either be absorbed, scattered, or emitted. The radiative energy absorbed and scattered at position $x$ are quantified by an absorption coefficient $κ_ν$ and a scattering coefficient $s_ν$, respectively, both expressed in cm$^{-1}$. The energy emitted along the direction $k$ is quantified by an emissivity term $ϵ_ν(k)$. The energy scattered from any direction $k'$ back to $k$ within a solid angle $dΩ$ is quantified with a term $\tilde{s}_ν(x, k', k)$. Neglecting time dependency terms, the radiative transfer equation reads

$$ \frac{\partial I_ν}{\partial x}(x, k) = -I_ν(x, k)[\kappa_ν(x) + s_ν(x)] + \epsilon_ν(x, k) + \frac{1}{4\pi} \int I_ν(x, k')\tilde{s}_ν(x, k', k)dΩ. \quad (1.1) $$

To simplify Eq. 1.1, the variable change $dτ_ν = (κ_ν + s_ν)dx$ is usually performed, where $τ_ν$ is called the optical depth. Neglecting the integral term, Eq. 1.1 then becomes:

$$ \frac{dI_ν}{dτ_ν} = -I_ν + S_ν, \quad (1.2) $$

with $S_ν = ϵ_ν/(κ_ν + s_ν)$ the source function. For some initial position $x_0$ – usually the back of the cloud – it reads (Draine, 2011, Eq. 7.14)

$$ \tau_{\nu}(x) = \int_{x_0}^{x} d\tau_{\nu}(x') dx' = \int_{x_0}^{x} [\kappa_{\nu}(x') + s_{\nu}(x')] dx'. \quad (1.3) $$

The amount of matter along the line of sight can be quantified through radiative transfer and extinction. The optical depth quantifies the attenuation of the radiative field at frequency $ν$ after passing through a cloud. It is positive except in case of optical pumping and MASER effect. When $τ_ν \gtrsim 1$, the cloud is said to be *optically thick*. A photon emitted inside the cloud is very likely to be re-absorbed – a phenomenon called *radiative trapping* (Draine, 2011, chapter 19). Therefore, an observer only detects photons emitted by the border of the cloud. Conversely, when $τ_ν \ll 1$, the cloud is said to be *optically thin*, and a photon is very likely to escape. In this case, an observer detects photons emitted in the entire cloud. The V (“visual”) band is a reference frequency band that corresponds to a filter with a 551 nm effective wavelength and an 88 nm full width at half maximum. One then commonly uses the optical depth in the V band $τ_V$ to measure the interstellar reddening, i.e., the extinction of visible light by interstellar dust. As astronomers commonly express intensities as magnitudes – which correspond to the decimal logarithm of physical fluxes – the extinction is defined in terms of a visual extinction

$$ A_V = [2.5 \times \log_{10} e] \quad \tau_V \approx 1.086 \tau_V. \quad (1.4) $$

Clouds with $A_V < 1$ are not deep enough to form molecules, as they get dissociated by the UV
photons. Dense cores can reach values of $A_V$ of the order $10^2$.

The proton column density, denoted $N_H$ and expressed in $cm^{-2}$, counts the protons along the line of sight. For a cloud of depth $x^{\text{tot}}$ (in cm),

$$N_H = \int_{0}^{x^{\text{tot}}} n_H dx. \qquad (1.5)$$

In a cloud with constant volume density $n_H$, $N_H = x^{\text{tot}} \times n_H$. As dust is the dominant source of extinction in the V band and as one commonly assumes that dust is uniformly mixed with the gas, the column density is proportional to the visual extinction

$$N_H = \frac{R_V}{C_D} A_V \approx 1.086 \frac{R_V}{C_D} \tau_V, \qquad (1.6)$$

where $C_D = 5.8 \times 10^{21} cm^{-2}$ (Bohlin et al., 1978), and $R_V$ is the so called total-to-selective extinction ratio and represents the slope of the extinction curve in the visible range. Typically, $R_V = 3.1$ in local diffuse clouds (Fitzpatrick and Massa, 1990).

**Radiative field and $G_0$** – As already stated, PDRs are generally illuminated by a massive star or a cluster of stars. The shape and intensity of the UV radiative field reaching the surface of the PDR depends on the properties of the star and on the distance between the PDR and the star. The $G_0$ parameter quantifies the intensity of the FUV radiative field regardless of its shape. It is defined as the ratio between the integral of the spectrum of the radiative field over a selected wavelength range and of the same integral computed for a reference interstellar radiation field (ISRF). For instance, prescriptions of ISRF can be found in Habing (1968), Draine (1978), and Mathis et al. (1983). Conversions of $G_0$ exist between the different ISRF prescriptions. In this thesis, we measure $G_0$ in reference to the Habing field. Far from any star in our Galaxy, $G_0 \approx 1$. In a moderately illuminated cloud such as the Horsehead nebula, $G_0 \approx 10^2$. In a highly illuminated cloud such as the Orion bar, $G_0 \approx 10^4$.

## 1.2 Linking physical conditions and observations: astrophysical numerical models

Astrophysical codes for ISM environments can model observed regions and link numerous observables such as line intensities to physical conditions such as the thermal pressure or the total visual extinction. There exist two main categories of codes: 3D (magneto)hydrodynamics simulations that model the dynamics at large scales, and detailed models of the local physico-chemistry of the gas and dust and of their interaction with the radiation field. In this thesis, we consider the latter category because they solve the microphysics of the cloud, compute chemical abundances and line intensities that can be compared to observations. In this section, we list different types of such numerical models of the ISM.

### 1.2.1 Radiative transfer models

FUV photons heat the gas by photoelectric effect on grains and excite atoms and molecules in electronic states. Molecules cascade in part in their rovibrational levels by emitting photons in FIR and radio domains. Radiative transfer thus carries two fundamental questions regarding radiative feedback in PDRs:

- What fraction of the stellar UV photons reaches a given position inside the cloud?

- What fraction of the photons emitted at a given position of the cloud escape and contribute to the observable emission and cooling, both for the dust continuum and in ionic, atomic, and molecular lines?
Answering these two questions requires accounting for absorption, scattering, and emission processes for the dust and the gas at each point of the cloud. Taking into account all these processes is challenging. First, the radiative transfer is coupled with quantum level populations of the emitting species for lines, and with the dust temperature for the continuum. Besides, the continuum and the lines are not independent. For instance, the dust continuum can pump the water molecules. Approximations are thus often used in radiative transfer models. A few remarkable models include RADEX, MOLPOP-CEP and Monte Carlo simulators.

RADEX<fnref n="1" /> (van der Tak et al., 2007) is a simple code that considers the cloud as spatially uniform in terms of density, temperature, chemical abundances and level populations, and is thus called a 0D code. This code solves the populations in quantum states in non local thermal equilibrium (non-LTE). It takes into account collisional and radiative transitions, but neglects interspecies interactions or pumping by any other line or by the dust continuum. As it models the cloud as completely uniform, it has to rely on a simple approximation for radiative trapping: a photon emitted within the cloud either escapes from the cloud or is absorbed on-the-spot. The code proposes different escape probability approximations. As this code is easy to use, it is widespread in the ISM community despite the coarse approximations it relies on.

MOLPOP-CEP (Asensio Ramos and Elitzur, 2018) is a 1D code that smartly formulates the problem only in terms of level populations. This reformulation allows it to solve the problem of radiative transfer in the lines exactly, although neglecting overlaps between lines and pumping by the dust continuum.

Finally, Monte Carlo simulators approach the radiative transfer differently. Instead of solving approximately the radiative transfer equation, they simulate the stochastic propagation of many individual photons in the cloud, and compute statistics. Such codes include RATRAN<fnref n="2" /> (Hogerheijde and van der Tak, 2000), LIME (Brinch and Hogerheijde, 2010), MCFOST (Pinte et al., 2022), or RADMC-3D (Dullemond et al., 2012). As they address more complex geometries, they usually simplify the physical processes solved at each point.

### 1.2.2 Astrochemical models

Astrochemistry codes are time-dependent 0D simulators that compute the evolution of the chemical abundance of up to hundreds of species accounting for thousands of reactions. A chemical network groups the list of photo-ionization, photodissociation, two-body or more rarely three-body reactions<fnref n="3" />, either in gas phase or surface phase, i.e., on the surface of dust grains. The evolution of the volume density $n_X$ of a species X is then controlled by a differential equation:

$$ \frac{dn_X}{dt} = F_X - D_X \qquad (1.7) $$

where $F_X$ and $D_X$ are formation rate and destruction rates, respectively, i.e., they sum the rates of all reactions in the network that contribute to the formation or destruction of species X. Combined with initial conditions, this set of equations permits to compute the volume density of all the species at any time t. The stationary state can be obtained by solving $F_X = D_X$.

Solving such a system of equations is numerically simple. The complexity lies in building the chemical network and modeling the reactions in ice and surface phases with the associated kinetic constants. Databases such as the KIDA database<fnref n="4" /> (Wakelam et al., 2012) and the UMIST database<fnref n="5" /> (McElroy et al., 2013) gather the kinetic constants of large chemical networks.

There exist many astrochemical codes based on these databases. NAHOON is associated with the KIDA database and solves networks for the gas phase. This model is not suited to dense

<footnote n="1">1. https://home.strw.leidenuniv.nl/~moldata/radex.html</footnote>

<footnote n="2">2. https://home.strw.leidenuniv.nl/~michiel/ratran/</footnote>

<footnote n="3">3. Three-body reactions are rare in the gas phase of the ISM, as it is not dense enough.</footnote>

<footnote n="4">4. https://kida.astrochem-tools.org/</footnote>

<footnote n="5">5. http://udfa.ajmarkwick.net/</footnote>
molecular regions, which require more processes involving dust grains. More advanced models such as ASTROCHEM (Maret and Bergin, 2015) and UCLCHEM (Holdship et al., 2017) include e.g., freeze-out of species onto dust grains or non-thermal desorption of species from dust grains due to UV photons and cosmic rays. These two models can be used to model a wide variety of environments, including molecular gas and dense cores. More recent and sophisticated models include NAUTILUS<fnref n="6" /> (Ruaud et al., 2016). The results of such time-dependent astrochemical models highly depend on the initial conditions that are poorly known. These codes can yield unexpected results such as oscillations, as demonstrated with the CHIMES<fnref n="7" /> code (Roueff and Le Bourlot, 2020).

The main limitation of such codes is their 0D geometry. This is problematic for photo-reactions as they do not solve radiative transfer. These models are suited to deep, cold and dark clouds that are protected from UV radiation. In such environments, the chemistry is driven by cosmic rays. The associated characteristic times are hundreds of million of years (Indriolo et al., 2007). Therefore, such environments cannot be considered in a stationary state.

### 1.2.3 Dust models

Dust grains play three central roles in the physics and chemistry of the PDRs: they absorb UV radiations and thus protect from photoionization and photodissociation, they heat the gas through the photoelectric effect, and they permit the formation of complex molecules on their surface. They also emit a continuum in the FIR domain. An observation of these emissions is called a spectral energy density (SED).

Understanding the many properties of dust grains is an active area of research. There exist different grain population models. Widespread models such as THEMIS<fnref n="8" /> (Jones et al., 2013) or the Astrodust+PAH model (Hensley and Draine, 2023) involve grains mostly made of silicon, carbon, or a mixture. Their structure is either amorphous, crystalline, porous, or aggregated. Finally, the grain size is usually described with a distribution, either a power law, a lognormal distribution or a mixture thereof.

From a grain population model, dust emission models can compute the extinction and emissivity of grains and predict their SED. The most widespread dust emission model is DUSTEM<fnref n="9" /> (Compiègne et al., 2011). Figure 1.8 shows a comparison of an observed SED (in gray) in the diffuse high galactic latitude medium with predictions from DUSTEM. The contribution to the predicted SED from different types of grain is represented, including polycyclic aromatic hydrocarbons, small amorphous carbon grains, large amorphous carbon grains, and amorphous silicates.

### 1.2.4 Holistic models

The aforementioned numerical models focus on a specific subpart of the physics and chemistry of the ISM. Some other models adopt a more holistic approach and account for many physical phenomena at once as well as their couplings. For instance, CLOUDY<fnref n="10" /> (Ferland et al., 2017) simulates HII regions. The Paris-Durham code<fnref n="11" /> (Godard et al., 2019) and the MAPPINGS code<fnref n="12" /> (Sutherland et al., 2018) simulate shock-dominated regions, i.e., regions crossed by a shock wave caused e.g., by protostellar outflows, cloud collisions, supernovae, or galactic outflows. PRODIMO<fnref n="13" /> (Woitke et al., 2009) simulates protoplanetary disks. The Meudon PDR code (Le Petit et al., 2006) and KOSMA-τ<fnref n="14" /> (Röllig and Ossenkopf-Okada, 2022) describe UV-irradiated

<footnote n="6">6. https://kida.astrochem-tools.org/codes.html</footnote>

<footnote n="7">7. https://ism.obspm.fr/chimes.html</footnote>

<footnote n="8">8. https://www.ias.u-psud.fr/themis/THEMIS_model.html</footnote>

<footnote n="9">9. https://www.ias.u-psud.fr/DUSTEM/index.html</footnote>

<footnote n="10">10. https://trac.nublado.org/</footnote>

<footnote n="11">11. https://ism.obspm.fr/shock.html</footnote>

<footnote n="12">12. https://mappings.anu.edu.au/code/index.html</footnote>

<footnote n="13">13. https://prodimo.iwf.oew.ac.at/</footnote>

<footnote n="14">14. https://astro.uni-koeln.de/riechers/research/kosma-tau</footnote>
<figure><img src="image_15.png" /></figure>

Figure 1.8 – Illustration of a spectral energy density for dust. The contributions of dust grain subpopulations are represented, such as polycyclic aromatic hydrocarbonss (PAHs), small amorphous carbon grains (SamC), large amorphous carbon grains (LamC), and amorphous silicates (aSi). Taken from Compiègne et al. (2011).

medium at the edge of molecular clouds in star forming regions or diffuse gas, i.e., PDRs.

All these models compute the chemical and physical structure resulting from a coupled treatment of chemistry, radiative transfer, and thermal processes. To account for all these phenomena and their coupling while maintaining reasonable evaluation times, these models rely on approximations such as the stationary state or a 1D geometry – except for PRODIMO that handles a 2D geometry, but resorts to other approximations. Despite these approximations, these models are usually slower and more computationally expensive than the aforementioned codes that address only a specific subpart of the physics and chemistry.

## 1.3 The Meudon PDR code

In this section, we give a short overview of the Meudon PDR code<sup>15</sup> (Le Petit et al., 2006). It is considered state-of-the-art PDR model as it includes many physical phenomena, which also causes each evaluation to last a few hours. This code is used in the remainder of this thesis to model each pixel individually in map inversion procedures. It can simulate a large variety of physical environments, such as diffuse clouds, dense PDRs, diffuse gas in the galactic center, nearby galaxies, damped Lyman alpha systems, circumstellar disks, etc. Indeed, part of the observable emissions of each of these environments comes from molecular gas irradiated by UV photons. In particular, this code permits to study effects such as the radiative feedback of a newborn star on its parent molecular cloud.

### 1.3.1 Physics and chemistry taken into account

The Meudon PDR code simulates the physics and chemistry of neutral interstellar gas illuminated with a FUV radiation field and computes its stationary state. The cloud is assumed to be a plane-parallel slab of neutral gas and dust with finite thickness, and is thus modeled with a one-dimensional geometry – the Meudon PDR code is called a 1D code. From a set of

<footnote n="15">15. https://ism.obspm.fr</footnote>
physical conditions of the cloud, the code iteratively solves large systems of coupled equations to compute the gas temperature, the chemical abundances and the excitation of species at each point of an adaptive grid. These equations include radiative transfer, thermal balance and chemistry. These intermediate and local results permit for example to compute integrated intensities comparable with observations. A full evaluation is computationally intensive and typically lasts a few hours. Figure 1.9 lists some physical phenomena taken into account for each physical aspect.

<figure><img src="image_16.png" /><figcaption>Figure 1.9 – Summary of Meudon PDR code modeling components.</figcaption></figure>

Radiative transfer equations are solved at each position on an adaptive spatial grid in a fully wavelength-dependent approach. In addition, the angular dependence is handled using a Legendre polynomial decomposition. Multiple phenomena are accounted for such as the redistribution across angles by dust scattering and the continuum absorption by dust and gas, including absorption in UV lines of H and H₂. For emission lines, the escape probability approximation of de Jong et al. (1980) is used. An emitted photon is therefore either absorbed on the spot or leaves the cloud.

The temperature of the gas at each position is solved at thermal balance,

$$ \Gamma(T) = \Lambda(T), \qquad (1.8) $$

where Γ(T) and Λ(T) are the heating rates of all heating and cooling processes, respectively. This equation is solved at each position of the grid to compute the gas temperatures from the specific intensity of the obtained radiation field. The heating rate Γ(T) includes a large range of processes such as photoelectric heating, collisional de-excitation of UV-pumped H₂, photoionization and photodissociation heating, cosmic ray heating, and exothermic chemical reactions. The cooling rate Λ(T) is estimated from the non local thermal equilibrium (non-LTE) excitation in the quantum levels of the main species by considering radiative and collisional processes as well as chemical formation and destruction. In total, the rovibrational states are considered for a few tens of species, and electronic states for some species such as H₂. Additional processes can either heat or cool the gas, such as H₂ heating or gas-grain collisions.

Finally, the chemical composition of the gas is computed at each position by solving for the chemical stationary state of a network of 200 species and 3000 reactions. In the gas phase, the
code accounts for two body reactions, photo-reactions – with primary and secondary photons –, reactions induced by cosmic rays, and some three body reactions. On the surface of dust grains, adsorption and thermal and non-thermal desorption are included. Finally, in ices, the code includes two body reactions due to thermal agitation and tunnel effect diffusions. The chemical reaction network was built combining different sources including data from the KIDA database and the UMIST database as well as data from articles. For key photo-reactions, cross-sections are taken from Ewine van Dishoeck's photodissociation and photoionization database<fnref n="16" /> (Heays et al., 2017).

The successive resolution of these three highly coupled aspects – radiative transfer, thermal balance, and chemistry equations – is iterated until a global stationary state is reached. After a run, the code provides density profiles of the chemical species and the temperature profiles of both the grains and the gas. It also outputs the line intensities emerging from the cloud that can be compared to observations. As of version 7, yet to be released, a total of 5 409 line intensities are predicted from 40 species such as H₂, HD, C⁺, C, CO, ¹³CO, C¹⁸O, ¹³C¹⁸O, SO, HCO⁺, HCN, HNC, CH⁺, CN or CS.

### 1.3.2 Input parameters and considered environments

For an input set of physical parameters, the Meudon PDR code computes the structure of the cloud, emitted line integrated intensities, and the dust continuum intensities. The predicted and observable intensities can be compared with actual observations to constrain the input physical parameters. The physical parameters to be estimated in this thesis characterize the radiative feedback. These parameters are the thermal pressure that describes local conditions, the visual extinction that quantifies the total amount of matter along the line of sight, and the intensity of the incident UV radiative field illuminating the surface of the PDR. Parameters associated with cosmic rays and dust grains are set to default values. Table 1.3 lists these default values.

<table><caption>Table 1.3 – Secondary input parameters in the Meudon PDR code and their default values.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Unit</th><th>Note</th></tr></thead><tbody><tr><td>cosmic rays ionization rate</td><td>10<sup>-16</sup></td><td>s<sup>-1</sup></td><td>Le Petit et al. (2004), Indriolo et al. (2007)</td></tr><tr><td>Dust extinction curve</td><td>Galaxy</td><td>–</td><td>Fitzpatrick and Massa (1990)</td></tr><tr><td>R<sub>V</sub></td><td>3.1</td><td>–</td><td>Fitzpatrick and Massa (1990)</td></tr><tr><td>C<sub>D</sub> = N<sub>H</sub>/E(B-V)</td><td>5.8 &times; 10<sup>21</sup></td><td>cm<sup>-2</sup></td><td>Bohlin et al. (1978)</td></tr><tr><td>Mass grain/Mass gas</td><td>0.01</td><td>–</td><td>–</td></tr><tr><td>Distribution on grain size a</td><td>&prop; a<sup>-3.5</sup></td><td>–</td><td>Mathis et al. (1977)</td></tr><tr><td>Min grain radius a<sub>min</sub></td><td>10<sup>-7</sup></td><td>cm</td><td>–</td></tr><tr><td>Max grain radius a<sub>max</sub></td><td>3 &times; 10<sup>-5</sup></td><td>cm</td><td>–</td></tr><tr><td>Turbulent velocity</td><td>2.0</td><td>km s<sup>-1</sup></td><td>–</td></tr></tbody></table>

**Constant density and constant pressure models** – The Meudon PDR code handles constant density and constant pressure models. Constant density models result in large pressure gradients in the cloud, as the temperature strongly drops in the H / H₂ transition region. They are therefore not very consistent with the stationarity hypothesis. However, they are quite common in astrophysics as they are easier to interpret. Constant pressure models are more consistent with the stationarity hypothesis and generally considered more realistic (Marconi et al., 1998; Lemaire et al., 1999; Allers et al., 2005; Goicoechea et al., 2016; Joblin et al., 2018; Wu et al., 2018). We restrict ourselves to constant pressure models. However, the thermal pressure is harder to physically interpret than volume density, which plays a role both in chemistry and radiative

<footnote n="16">16. https://home.strw.leidenuniv.nl/~ewine/photo/index.html</footnote>
transfer – for instance with critical densities.

**Total visual extinction** – As already mentioned, the Meudon PDR code relies on an adaptive one-dimensional grid of a plane-parallel slab of gas. The grid of the Meudon PDR is a grid on the visual extinction $A_V$, which quantifies the amount of matter along the line of sight. The total depth of the cloud is thus defined with a total visual extinction $A_V^{tot}$.

**Radiation field** – The Meudon PDR code permits the use of any incident radiation field, including a stellar radiation field. For simplicity, in this thesis, we consider the radiation field to be the Mathis ISRF Mathis et al. (1983) scaled at the front of the cloud by one scalar parameter denoted radm. The conversion of radm to $G_0$ is done with

$$G_0 = \frac{1.2786}{2} \text{radm}, \qquad (1.9)$$

where the 1.2786 factor comes from the conversion from the Mathis to Habing ISRF, and the $\frac{1}{2}$ factor from the fact that the front face of the cloud is illuminated from the outside of the cloud, and not from the inside.

### 1.3.3 Limitations

The geometry simplification can affect integrated intensities and the effective depth of penetration of UV photons. For instance, under an isotropic external radiation field, UV photons should penetrate deeper in a spherical cloud than in a plane-parallel cloud. Other PDR codes assume different cloud geometries, e.g., the KOSMA-$\tau$ numerical code (Röllig and Ossenkopf-Okada, 2022) assumes a 1D spherical geometry.

The stationarity assumption neglects gas dynamics and non-equilibrium effects. In other words, the timescales of all physical and chemical processes taken into account must be smaller than that of the changes in the external conditions, and of the gas dynamics. In a PDR, the physical processes are dominated by the photodissociation of $H_2$ by UV photons. The $H_2$ photodissociation rate at the surface of a cloud is $2.9 \times 10^{-11} G_0 s^{-1}$. For $G_0 = 10^4$, the characteristic time of $H_2$ photodissociations is about 40 days. As this duration is very short compared to gas dynamics timescales, the stationarity hypothesis is realistic. However, the UV photons do not penetrate much beyond the photodissociation front. In deeper clouds, the $H_2$ dissociation is dominated by cosmic rays. As indicated in Table 1.3, the cosmic ray ionization rate is $\sim 10^{-16} s^{-1}$ (Indriolo et al., 2007). The corresponding characteristic time is about 300 Myr. Therefore, such clouds require time-dependent chemistry models.

Finally, the neutral gas assumption means that in the modeled region, hydrogen is either in its atomic form ($H$) or its molecular form ($H_2$) but never in its ionized form ($H^+$), even in the atomic layer, close to the ionization front. Other PDR codes that relax these hypotheses exist. For instance, the Hydra PDR code (Bron et al., 2018a) relaxes the stationarity and neutrality assumptions by including gas hydrodynamics and photo-evaporation. However, this code simplifies the modeling of other phenomena compared to the Meudon PDR code in order to maintain reasonable computational durations, i.e., a few hours per evaluation.

## 1.4 Conclusion

The ISM is a complex and inhomogeneous medium that fills the unimaginably large volume that lies between stars. It is made mostly of hydrogen and helium. It also contains traces of other elements and small dust grains that play a key role in its chemistry and physics. In this thesis, we focus on one particular type of environment called photodissociation region (PDR) that lies at the border between a massive star or cluster of stars and a molecular cloud. Studying PDRs permits to understand better the impact of radiative feedback from newborn stars on their parent cloud and its star forming capability.
Observations with high spatial resolution such as the Orion B map described in Pety et al. (2017), the OMC-1 map described in Goicoechea et al. (2019), and future observations from ALMA or the JWST contain much physical information regarding the radiative feedback. In this thesis, our goal is to extract this information. Our aim is to infer maps of the key physical parameters from such large multiline PDR observations. We model PDRs with the Meudon PDR code. As it iteratively solves large, coupled, and complex non-linear sets of equations, each evaluation of the Meudon PDR code requires from a few hours to a few days. In Chapter 4, we derive a fast and accurate neural-network based approximation of the Meudon PDR code to perform fast inference on large observation maps. In Chapter 5, we propose a new statistical inference algorithm that exploits this Meudon PDR code approximation, and validate it on synthetic data. In Chapter 6, this statistical inference algorithm is applied to multiple real observations.
# Chapter 2

# Background on Bayesian statistical modeling and inference

> "The theory of probabilities is at bottom nothing but common sense reduced to calculus; it enables us to appreciate with exactness that which accurate minds feel with a sort of instinct for which of times they are unable to account."

Pierre-Simon Laplace, in Théorie Analytique Des Probabilités, 1814

## Contents

<table><tr><td>2.1 Bayesian statistical modeling</td><td>30</td></tr><tr><td>2.1.1 Main random variables and probability distributions</td><td>30</td></tr><tr><td>2.1.2 Estimators</td><td>32</td></tr><tr><td>2.2 Statistical inference</td><td>33</td></tr><tr><td>2.2.1 Evaluating estimators defined as solutions of optimization problems</td><td>34</td></tr><tr><td>2.2.1.1 Gradient descent (GD) algorithm</td><td>34</td></tr><tr><td>2.2.1.2 Preconditioned variants</td><td>36</td></tr><tr><td>2.2.1.3 Block coordinate variant</td><td>38</td></tr><tr><td>2.2.2 Approximating integrals with Monte Carlo estimators</td><td>39</td></tr><tr><td>2.2.2.1 Monte Carlo (MC) estimators</td><td>40</td></tr><tr><td>2.2.2.2 From Monte Carlo to Markov chain Monte Carlo (MCMC) methods</td><td>40</td></tr><tr><td>2.2.2.3 Metropolis-Hastings (MH) algorithm</td><td>43</td></tr><tr><td>2.2.2.4 Using the gradient with Metropolis adjusted Langevin algorithm (MALA)</td><td>44</td></tr><tr><td>2.2.2.5 Preconditioning with MALA to improve mixing properties</td><td>45</td></tr><tr><td>2.2.2.6 Gibbs sampling: a block coordinate MH variant</td><td>46</td></tr><tr><td>2.2.2.7 Generating multiple candidates</td><td>48</td></tr><tr><td>2.3 Comparing and checking observation models</td><td>49</td></tr><tr><td>2.3.1 Model selection and evaluation</td><td>50</td></tr><tr><td>2.3.2 Model checking using Bayesian hypothesis testing</td><td>52</td></tr><tr><td>2.4 Conclusion</td><td>54</td></tr><tr><td>Appendix 2.A Samplers widespread in astrophysics dedicated to multimodal distributions</td><td>56</td></tr><tr><td>2.A.1 Adapting meta-heuristics to MCMC</td><td>56</td></tr><tr><td>2.A.2 Prior or parallel identification of the modes</td><td>58</td></tr><tr><td>2.A.3 Nested sampling</td><td>59</td></tr></table>
In Chapter 1, we described the interstellar medium (ISM) and some numerical models that can be used to infer physical parameters from observations. This chapter describes notions and tools necessary to formulate and solve the inverse problem considered in this thesis – see Chapter 5. Some of these methods are already popular in ISM studies, as we will show in Chapter 3.

Section 2.1 describes the general Bayesian statistical modeling approach, the associated random variables, probability distributions and estimators. Section 2.2 presents two common approaches for estimator derivation, namely optimization methods and sampling methods. In particular, methods based on first and second order derivatives are covered to exploit the twice differentiable log-posterior and address its absence of gradient Lipschitz continuity. Block coordinate algorithms are also described to exploit the natural structure of the parameter $\Theta$. The sampler proposed in Chapter 5 (Section 5.2) builds on the methods presented in this section. Finally, Section 2.3 depicts two methods that assess the compatibility of the model with the observations, namely model selection and model checking.

## 2.1 Bayesian statistical modeling

This section covers fundamental notions of Bayesian statistical modeling. First, the main variables and probability distributions involved in statistical modeling are reviewed, including the likelihood function and posterior distribution. Then, the main estimators are described.

### 2.1.1 Main random variables and probability distributions

The Bayesian framework is a statistical approach that encodes and describes uncertainty using random variables and probability distributions. In this context, a probability distribution indicates a of belief and not a limit of frequencies as in frequentist statistics. Bayesian statistics exploits random variables to model uncertainty in any quantity of interest such as observations, physical parameters, or model hyperparameters. For a review on random variables and probability distributions, see Le Gall (2022).

There are two main random variables of interest in Bayesian statistics: the physical parameters to infer $\Theta \in \mathbb{R}^{N \times D}$ and the observation $\mathbf{Y} \in \mathbb{R}^{N \times L}$. Assuming an observation model $\mathcal{M}$, the Bayes theorem states that the conditional probability density functions (pdfs) of these variables verify

$$ \pi(\Theta|\mathbf{Y}, \mathcal{M}) = \frac{\pi(\mathbf{Y}|\Theta, \mathcal{M}) \pi(\Theta)}{\pi(\mathbf{Y}|\mathcal{M})}, \qquad (2.1) $$

where $\pi(\mathbf{Y}|\Theta, \mathcal{M})$ is the likelihood function, $\pi(\Theta)$ is the pdf of the prior distribution, $\pi(\Theta|\mathbf{Y}, \mathcal{M})$ is the pdf of the posterior distribution, and $\pi(\mathbf{Y}|\mathcal{M})$ is the Bayesian evidence.

The likelihood function $\pi(\mathbf{Y}|\Theta, \mathcal{M})$ evaluates how well the parameter $\Theta$ and the model $\mathcal{M}$ reconstruct the observations $\mathbf{Y}$. It is the pdf of a distribution on the observation $\mathbf{Y}$ given a fixed value for the physical parameters $\Theta$. In inference, the observation $\mathbf{Y}$ is fixed and the free variable is $\Theta$. The denomination “likelihood function” emphasizes that $\Theta \mapsto \pi(\mathbf{Y}|\Theta, \mathcal{M})$ is not a pdf. The likelihood function is derived from an observation model $\mathcal{M} = (\mathbf{f}, \mathcal{A})$ that groups a forward model $\mathbf{f}$, generally a numerical simulator in astrophysics, and a model noise $\mathcal{A}$ such that

$$ \mathbf{Y} = \mathcal{A}(\mathbf{f}(\Theta)). \qquad (2.2) $$
## Details for astrophysicists 2.1: Example: Gaussian additive and uncorrelated noise

For a Gaussian additive and uncorrelated noise, the observation model in Eq. 2.2 reads

$$ \forall n \in [1, N], \quad \ell \in [1, L], \quad y_{n\ell} = f_{\ell}(\theta_n) + \varepsilon_{n\ell}^{(a)}, \quad \varepsilon_{n\ell}^{(a)} \sim \mathcal{N}(0, \sigma_{n\ell}^2). \quad (2.3) $$

Using the pdf of a Gaussian distribution, this observation model translates to the likelihood function

$$ \pi(\mathbf{Y}|\mathbf{\Theta}, \mathcal{M}) \propto \prod_{n=1}^{N} \prod_{\ell=1}^{L} \exp \left[ -\frac{1}{2\sigma_{n\ell}^2} (f_{\ell}(\theta_n) - y_{n\ell})^2 \right], \quad (2.4) $$

and to the negative log-likelihood function

$$ -\ln \pi(\mathbf{Y}|\mathbf{\Theta}, \mathcal{M}) = \sum_{n=1}^{N} \sum_{\ell=1}^{L} \frac{1}{2\sigma_{n\ell}^2} (f_{\ell}(\theta_n) - y_{n\ell})^2, \quad (2.5) $$

up to an additive constant. The $\chi^2$ loss function, widespread in ISM studies, is equivalent to the negative log-likelihood function in Eq. 2.5 for this specific case of Gaussian additive and uncorrelated noise. The negative log transform turns the product to a sum and removes the exponential, which makes derivative computations easier.

Note that in this thesis, the natural logarithm of $x \in \mathbb{R}$ is denoted $\ln x$. This notation is preferred to the usual log to avoid any confusion between logarithms in base e and 10.

In this thesis, the forward model $f$ is set to the Meudon PDR code, introduced in Chapter 1 (Section 1.3). This forward model maps physical parameter vectors $\theta \in \mathbb{R}^D$ to observables $f_{\ell}(\theta) \in \mathbb{R}^L$. Therefore, evaluating the Meudon PDR code for a map $\mathbf{\Theta} = (\theta_n)_{n=1}^N$ requires $N$ evaluations: $f(\mathbf{\Theta}) = (f(\theta_n))_{n=1}^N$. To accelerate the inference procedure, we proposed an approximation of the Meudon PDR code based on an artificial neural network with controlled error. The derivation of this approximation is the core of Chapter 4. For a review of how numerical models are handled in inference procedures in ISM studies, see Chapter 3 (Section 3.1.1). Similarly, Section 3.4 introduces the noise model considered in this thesis. For a review on noise models used in astrophysics, see Section 3.1.2.

The **prior distribution** $\pi(\mathbf{\Theta})$ encodes prior information on the physical parameters of interest $\mathbf{\Theta}$. It may be non-informative, such as a uniform distribution on a set of parameters, or more informative, e.g., to include the results of other trusted studies or to constrain acceptable solutions to be physically meaningful. Similarly to the likelihood function, using the negative log-prior $-\ln \pi(\mathbf{\Theta})$ simplifies computations. In the optimization framework, $-\ln \pi(\mathbf{\Theta})$ is often called a regularization function.

Chapter 3 (Section 3.1.3) reviews prior distributions used in ISM studies. Chapter 5 (Section 5.1.2) describes the prior distribution considered in this thesis.

The **posterior distribution** $\pi(\mathbf{\Theta}|\mathbf{Y}, \mathcal{M})$ combines the likelihood function and the prior distribution. Like the prior distribution, it is a distribution on $\mathbf{\Theta}$. The posterior describes the parameters of interest knowing the observations. It is the distribution from which estimators of the physical parameters $\mathbf{\Theta}$ are extracted, as well as the associated uncertainty. The posterior distribution associated with the inverse problem considered in this thesis is introduced in Chapter 5 (Section 5.1.3).

In some simple cases, the likelihood function and prior distribution enjoy a *conjugacy relation* (Robert and Casella, 2004, chapter 1). The posterior distribution is then an element of a
parametric family of distributions with closed-form formulae for its parameters. However, these simple cases are very limited and do not occur in astrophysical inference problems. For instance, there exists no conjugacy relation when the forward model is non-linear.

The Bayesian evidence $\pi(\mathbf{Y}|\mathcal{M})$ is a normalization constant that does not depend on the physical parameter $\Theta$ and that is generally challenging to evaluate. It is also called marginal likelihood, as it marginalizes over the physical parameter $\Theta$

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

Many approaches that solve inverse problems require the assumption of a model $\mathcal{M} = (\mathbf{f}, \mathcal{A})^1$. The Bayesian evidence permits to assess the consistency of the model with the observations. It is a crucial quantity for Bayesian model selection, as we will discuss in Section 2.3.1. However, it is not always necessary for statistical inference. In inference, Eq. 2.1 is therefore often simplified to

$$ \pi(\Theta|\mathbf{Y}, \mathcal{M}) \propto \pi(\mathbf{Y}|\Theta, \mathcal{M}) \pi(\Theta), \qquad (2.7) $$

where $\propto$ means “proportional to”, i.e., equal up to a multiplicative constant independent of $\Theta$. The model $\mathcal{M}$ remains unchanged during statistical inference. Writing explicitly the dependence to the model $\mathcal{M}$ is useful when evaluating its relevance or comparing multiple models. Otherwise, this dependence is not explicitly written to avoid heavy notations. Therefore, Eq. 2.7 will often be simplified to

$$ \pi(\Theta|\mathbf{Y}) \propto \pi(\mathbf{Y}|\Theta) \pi(\Theta). \qquad (2.8) $$

**Hyperparameters** can have a dramatic impact on the shape of the posterior distribution. They can appear in both the prior distribution and the likelihood function. For instance, a Gaussian prior on $\Theta$ would be defined with a mean vector and a covariance matrix. Similarly, the likelihood function obtained with a Gaussian additive noise model is defined with a noise covariance matrix. These hyperparameters greatly influence the shape of the posterior distribution. Hierarchical models (Gelman et al., 2015, chapter 5) enable inferring simultaneously both $\Theta$ and hyperparameters from the observations $\mathbf{Y}$. For now, we assume hyperparameters set to constant values, and restrict ourselves to non-hierarchical models.

## 2.1.2 Estimators

Two main types of estimators can be defined from the likelihood function and the posterior distribution: some are solutions of an optimization problem, while others are defined as an integral.

**Estimators solutions of an optimization problem** – The point that maximizes the likelihood function is by definition the point that best reproduces observations. It is called the maximum likelihood estimator (MLE), and denoted $\hat{\Theta}_{MLE}$. It is a natural estimator with good theoretical properties (Robert and Casella, 2004, chapter 1). Similarly, the posterior mode is the point that best balances observation reproduction and consistency with prior knowledge. It is called the maximum a posteriori (MAP), and denoted $\hat{\Theta}_{MAP}$. The MLE and the MAP are defined as solutions of a maximization problem, often converted to a minimization problem in log-scale:

$$ \hat{\Theta}_{MLE} \in \arg\max_{\Theta \in \mathbb{R}^{N \times D}} \pi(\mathbf{Y}|\Theta) = \arg\min_{\Theta \in \mathbb{R}^{N \times D}} -\ln \pi(\mathbf{Y}|\Theta), \qquad (2.9) $$

and

$$ \hat{\Theta}_{MAP} \in \arg\max_{\Theta \in \mathbb{R}^{N \times D}} \pi(\Theta|\mathbf{Y}) = \arg\min_{\Theta \in \mathbb{R}^{N \times D}} -\ln \pi(\Theta|\mathbf{Y}). \qquad (2.10) $$

<footnote n="1">1. An example of approach that does not rely on the assumption of a model resorts to end-to-end neural networks – see e.g., conditional invertible neural networks (Ardizzone et al., 2019).</footnote>
## Details for astrophysicists 2.2: Well-posed and ill-posed inverse problems

Inverse problems are usually divided in two classes: *well-posed problems* and *ill-posed problems*.

In well-posed problems, the observations alone are sufficiently informative to yield physical parameter accurate estimations. In astrophysics, these problems are said to be *well constrained*. The MLE is very useful in such problems.

Conversely, in ill-posed problems, the observations may come in low SNR regime, or the observables may not be good tracers of the physical parameter of interest. In astrophysics, such problems are called *not constrained or badly constrained*. In such cases, the MLE leads to unstable results. By exploiting regularizing prior knowledge, the MAP can be more robust to noise and to non-linearities.

**Estimators defined as an integral** – Many estimators are defined as an integral, including the posterior mean and covariance. The posterior mean $E[\Theta|\mathbf{Y}, \mathcal{M}]$ is the point estimate that minimizes the mean squared error $E[\|\Theta - \hat{\Theta}\|_2^2 | \mathbf{Y}, \mathcal{M}]$, i.e., the average squared distance with other points $\Theta$ drawn from the posterior. This estimator is therefore usually called the minimum mean square error (MMSE), and denoted $\hat{\Theta}_{\text{MMSE}}$. The posterior covariance $\text{Cov}(\Theta|\mathbf{Y}, \mathcal{M})$ permits the detection of degeneracies between physical parameters<fnref n="2" />. In general, multiple estimators can be written as the integral of a function $g$ on $\Theta$ with respect to the posterior,

$$E[g(\Theta)|\mathbf{Y}] = \int g(\Theta) \pi(\Theta|\mathbf{Y}) d\Theta. \qquad (2.11)$$

For instance, the posterior mean $E[\Theta|\mathbf{Y}]$ corresponds to $g: \Theta \mapsto \Theta$, i.e., the identity function. Similarly, the posterior covariance matrix corresponds to $g: \Theta \mapsto (\Theta - E[\Theta|\mathbf{Y}])(\Theta - E[\Theta|\mathbf{Y}])^T$.

A more accurate uncertainty quantification on physical parameters than the posterior covariance can be obtained with a credible region (Pereyra, 2017). A credible region $C_\alpha$ with confidence level $1 - \alpha$ verifies

$$\mathbb{P}[\Theta \in C_\alpha] = \int 1_{C_\alpha}(\Theta) \pi(\Theta|\mathbf{Y}) d\Theta = 1 - \alpha. \qquad (2.12)$$

A posterior distribution admits infinitely many credible regions $C_\alpha$. The *highest posterior density region* $C_\alpha^*$ is the credible region with minimum volume. It is defined as

$$C_\alpha^* = \{\Theta | -\ln \pi(\Theta|\mathbf{Y}) \le c_\alpha\} \qquad (2.13)$$

where $c_\alpha \in \mathbb{R}$ is set so that Eq. 2.12 is verified with $C_\alpha^*$.

For one parameter $\theta_{nd}$, a credibility interval either corresponds to the highest posterior density region $\pi(\theta_{nd}|\mathbf{Y})$ – with all other parameters being marginalized – or to an interval defined from percentiles. In this work, the $1 - \alpha$ credibility interval is defined with the $\alpha/2$ and $1 - \alpha/2$ percentiles to simplify the presentation of results. These percentiles can be obtained with the inverse cumulative density function (cdf) of $\pi(\theta_{nd}|\mathbf{Y})$. Note that when the marginalized posterior $\pi(\theta_{nd}|\mathbf{Y})$ is unimodal and symmetric, both definitions coincide.

## 2.2 Statistical inference

This section covers fundamental methods of statistical inference that numerically compute the aforementioned estimators. Estimators defined as a minimum of a loss function are evaluated

<footnote n="2">². A degeneracy between physical parameters corresponds to an *ill-conditioned* Hessian matrix. The notion of condition number of a matrix will be introduced in Section 2.2.1.2.</footnote>
through optimization methods. Estimators defined with potentially high dimensional integrals are approximated with Monte Carlo (MC) estimators, which requires sampling from the posterior distribution. Unlike the optimization approach, this sampling approach thus permits the use of one algorithm to evaluate multiple estimators. It can therefore extract more information from the posterior distribution, in particular regarding uncertainty. In this thesis, due to the absence of ground truth in ISM studies, we address inverse problems with a sampling approach.

The optimization approach, that we do not use for inference, is presented for three reasons. First, it is more intuitive than the sampling approach, and many links can be made between their respective methods. We believe that knowing optimization methods simplifies the understanding of sampling methods for the unfamiliar reader. Second, some proposed sampling methods are inspired from optimization methods that thus need to be introduced. Finally, as we will show in Chapter 3 (Section 3.2.1), optimization methods are widely exploited in ISM studies.

This section only covers methods that are used to derive the inference algorithm proposed in Chapter 5. In particular, variational Bayes inference and proximal methods are not reviewed. Variational Bayes inference is an intermediate approach between optimization and sampling. It consists in projecting the true posterior onto a class of simple distributions. This target class is generally chosen so that it is easy to extract information. For instance, for a Gaussian approximation, the mean vector and covariance matrix are adjusted to approximate the posterior distribution. Other insights such as credibility intervals are then very simple to derive from the approximating distribution. See e.g., Pereyra et al. (2016) for a review. These methods require prior knowledge of the shape of the posterior distribution to select a relevant class of approximating distribution. Because of the non-linear and non gradient Lipschitz continuous forward model, and of the complex noise model, we do not have such prior knowledge. Besides, the quality of the approximation greatly impacts the relevance of the provided uncertainty quantification. Therefore, we decided not to consider these methods, and do not introduce them formally in this thesis.

Proximal gradient descent and proximal sampling algorithms (Pereyra et al., 2016) are now widespread methods that permit to address cases where the posterior distribution is not differentiable. In the inverse problem considered in this thesis, the log-posterior is twice differentiable. Algorithms based on the proximal operator are therefore not covered in this thesis.

## 2.2.1 Evaluating estimators defined as solutions of optimization problems

Some estimators $\hat{\Theta}$ are defined as the minimum of a loss function, denoted $\mathcal{L}$, as described in the previous section. For instance, this loss function $\mathcal{L}$ corresponds to the negative log-likelihood for the MLE and to the negative log-posterior for the MAP. In most cases, such estimators cannot be written in closed-form. They are therefore numerically approximated. The most naive method is to look for the precomputed model that leads to the lowest loss value among a set of precomputed models. This method is usually called a grid search. It scales very poorly as the number of required precomputed models increases exponentially with the dimension. A wide variety of methods permit better and cheaper evaluation of these two point estimations. In the following, we describe the main alternative method, the gradient descent (GD) algorithm, and its preconditioned and block coordinate variants.

### 2.2.1.1 Gradient descent (GD) algorithm

In its simplest version, the gradient descent (GD) algorithm starts at some state $\Theta^{(0)}$, and updates it with

$$ \Theta^{(t)} = \Theta^{(t-1)} - \eta \nabla \mathcal{L}(\Theta^{(t-1)}), \quad (2.14) $$

with $\eta > 0$ a step size hyperparameter.

The matrices $\Theta \in \mathbb{R}^{N \times D}$ are converted to vectors of $\mathbb{R}^{ND}$ in lexicographic order to evaluate the gradient. In the remainder of this thesis, this vector form of the physical parameter $\Theta$ is
considered whenever a first or second order derivative is used. Iterations of the GD algorithm are repeated until some stopping criterion is satisfied. This stopping criterion typically corresponds to a maximum number T of iterations or to a lower threshold on the gradient norm ||∇L (Θ(t)) || or, equivalently, on the distance of two successive iterates ||Θ(t) − Θ(t−1)||. Algorithm 2.1 summarizes the GD algorithm.

<caption>**Algorithm 2.1:** Gradient descent (GD) algorithm</caption>

**Input:** Starting point Θ⁽¹⁾, step size η > 0.
1 t = 0
2 while Stopping criterion not satisfied do
3 Update Θ⁽t⁾ with Θ⁽t⁻¹⁾ and ∇L // using Eq. 2.14, Eq. 2.17 or Eq. 2.22
4 t = t + 1

**Output:** Last iterate Θ⁽t⁾

The gradient descent (GD) algorithm (Shalev-Shwartz and Ben-David, 2014, chapter 14) is generally preferred to grid search, as it generally yields better estimates $\hat{\Theta}$, i.e., estimates with lower loss value $L(\hat{\Theta})$, with less function evaluations. This advantageous comparison is due to two strengths of the GD algorithm. First, it exploits more information from the loss function $L$, with its first and potentially second order derivative information. Second, its iterative nature exploits past evaluations of the loss functions to progressively converge to a critical point, potentially a local or global minimum.

**Convergence properties** – GD algorithms converge to critical points, i.e., points Θ such that $\nabla L(\Theta) = 0$, under different sets of assumptions. Theoretical guarantee of convergence can be obtained for a loss function whose gradient $\nabla L$ is a Lipschitz continuous function, i.e., if there exists a constant $\beta > 0$ such that

$$ \forall \Theta_1, \Theta_2, \quad \|\nabla L(\Theta_1) - \nabla L(\Theta_2)\| \le \beta \|\Theta_1 - \Theta_2\|. \quad (2.15) $$

Such loss functions $L$ are called *gradient Lipschitz continuous*, with *Lipschitz constant* $\beta$.

**Details for astrophysicists 2.3: Two examples for gradient Lipschitz continuity**

**Example 1: gradient Lipschitz continuous** – The gradient of the $x \in \mathbb{R} \mapsto x^2$ is $x \mapsto 2x$, which admits 2 as a Lipschitz constant. The function $x \in \mathbb{R} \mapsto x^2$ is therefore gradient Lipschitz continuous with constant $\beta = 2$.

**Example 2: non gradient Lipschitz continuous** – The gradient of $x \in \mathbb{R} \mapsto \exp(x^2)$ is $x \mapsto 2x \exp(x^2)$, which does not admit a Lipschitz constant. This second function is therefore not gradient Lipschitz continuous.

For non-convex and gradient Lipschitz continuous loss functions, convergence to a critical point is guaranteed when the step size $\eta$ satisfies $\eta < 2/\beta$ (Beck, 2017, Theorem 10.15)<fnref n="3" />. Smaller step sizes lead to slower convergence, and the GD algorithm does not converge for values larger than the $2/\beta$ limit. In many imaging inverse problems, the log-posterior is gradient Lipschitz continuous with a known Lipschitz constant. In general, however, the log-likelihood and

<footnote n="3">3. This theorem holds for proximal gradient descent. As the gradient descent algorithm is a special case of proximal gradient descent, the result also holds for gradient descent.</footnote>
log-posterior are not gradient Lipschitz continuous, or the constant is either hard to compute or too large to be of any practical use.

### 2.2.1.2 Preconditioned variants

The update step from Eq. 2.14 can be prohibitively slow to converge, especially when $\Theta \in \mathbb{R}^{N \times D}$ lives in high dimension or when the dynamics of its components $\theta_{nd} - n \in [1, N]$, $d \in [1, D]$ – are very different. In this section, we introduce the condition number, that indicates how slowly the standard GD will converge. Then, we describe preconditioning techniques that accelerate convergence by locally reshaping the loss function.

**Condition number and preconditioning** – The condition number $\kappa$ of the Hessian matrix $\nabla^2 \mathcal{L}(\Theta)$ quantifies the difference in dynamic ranges between components $\theta_{nd}$ or their combinations. It is defined as

$$ \kappa(\nabla^2 \mathcal{L}(\Theta)) = \frac{\sigma_{\max}(\nabla^2 \mathcal{L}(\Theta))}{\sigma_{\min}(\nabla^2 \mathcal{L}(\Theta))} \geq 1, \qquad (2.16) $$

where $\sigma_{\max}$ and $\sigma_{\min}$ denote the maximal and minimal singular values of a matrix, respectively. These singular values are the absolute value of the eigenvalues. When $\kappa \approx 1$, the Hessian matrix is said to be *well-conditioned*. It equals one, e.g., for the identity matrix. Conversely, when $\kappa \gg 1$, the Hessian matrix is said to be *ill-conditioned*.

Figure 2.1 shows two ellipses, each corresponding to a contour level of a quadratic function, i.e., with constant Hessian matrices. In Figure 2.1a, the condition number is close to 1. There is a slight anisotropy, but the dynamic ranges are roughly the same for both variables. Conversely, in Figure 2.1b, the condition number is 8. There is a strong anisotropy, and the dynamic ranges are not the same for the two variables. The GD algorithm performs best in the first case, and converges slowly in the second case – see e.g., Nocedal and Wright (2006, theorem 3.4).

<figure><img src="image_17.png" /><figcaption>Figure 2.1 – Illustration of the condition number of the Hessian matrix for two two-dimensional quadratic functions. The ellipse corresponds to a contour line of the quadratic function.</figcaption></figure>

The goal of preconditioning is to transform Figure 2.1b into Figure 2.1a. To do so, a preconditioning matrix $G^{(t-1)} = G(\Theta^{(t-1)}) \in \mathbb{R}^{ND \times ND}$ added to Eq. 2.14

$$ \Theta^{(t)} = \Theta^{(t-1)} - \eta G^{(t-1)} \nabla \mathcal{L}(\Theta^{(t-1)}). \qquad (2.17) $$

In presence of a preconditioning matrix $G^{(t-1)}$, the condition number is measured on the new Hessian $G^{(t-1)} \nabla^2 \mathcal{L}(\Theta^{(t-1)})$. The goal when resorting to a preconditioner is thus to decrease as much as possible the condition number. In Guggenheimer et al. (1995), the authors provide a general upper bound on the condition number of an invertible matrix. They also include an analysis of situations where their bound is tight.
The inverse Hessian matrix $G^{(t-1)} = [\nabla^2 \mathcal{L}(\Theta^{(t-1)})]^{-1}$ is the natural choice of preconditioner. It yields the Newton algorithm (Nocedal and Wright, 2006, chapter 10) applied to $\nabla \mathcal{L}$. By definition the Newton algorithm looks for zeros of a function. Applying it on $\nabla \mathcal{L}$ boils down to finding critical points of the loss function $\mathcal{L}$. With this preconditioner, the condition number equals exactly 1. The modified Hessian matrix is thus perfectly conditioned. However, this preconditioner suffers from three main drawbacks. First, it can only be used when the Hessian matrix is invertible. Second, the Newton algorithm requires evaluating and inverting the full Hessian matrix, and then applying a matrix-vector product with the gradient. These three steps are very expensive, especially in high dimensions. A faster and more stable implementation directly evaluates the vector $\theta^{(t)} - \theta^{(t-1)}$ by considering it as the solution of a linear system:

$$[\nabla^2 \mathcal{L}(\Theta^{(t-1)})](\theta^{(t)} - \theta^{(t-1)}) = -\nabla \mathcal{L}(\Theta^{(t-1)}), \quad (2.18)$$

but remains expensive in high dimensions. Lastly, it can suffer from instabilities for loss functions that are not gradient Lipschitz continuous.

The Gauss-Newton and Levenberg-Marquardt algorithms (Nocedal and Wright, 2006, chapter 10) are preconditioned GD algorithms. As we will show in Chapter 3, Levenberg-Marquardt is quite common in ISM studies. These two algorithms resort to a cheaper and invertible approximation of the inverse Hessian that does not require the evaluation of second order derivatives. They enjoy similar convergence results and convergence rate as the Newton algorithm. However, they are specialized to non-linear least squares problems, i.e., to loss functions $\mathcal{L}$ of the form

$$\mathcal{L}: \Theta \in \mathbb{R}^{ND} \mapsto \sum_{n=1}^{N} \sum_{\ell=1}^{L} \frac{1}{2\sigma_{n\ell}^2} (f_{\ell}(\theta_n) - y_{n\ell})^2 \quad (2.19)$$

The Gauss-Newton algorithm uses $G^{(t-1)} = [\nabla \mathcal{L}(\Theta) \nabla \mathcal{L}(\Theta)^T]^{-1}$ to approximate the Hessian matrix. This approximation is symmetric by construction, and positive definite as long as the matrix $\nabla \mathcal{L}(\Theta) \nabla \mathcal{L}(\Theta)^T \in \mathbb{R}^{ND \times ND}$ is full rank. This condition requires $L \ge D$, i.e., the number $L$ of observables per pixel needs to be greater or equal to the dimension $D$ of the parameter $\theta$ per pixel. Though this condition is usually verified in inverse problems, it is not sufficient to guarantee this matrix is invertible. In some cases, the Gauss-Newton algorithm thus suffers from numerical instabilities, especially far from a local minimum. The Levenberg-Marquardt algorithm “damps” this approximation of the Hessian matrix: $G^{(t-1)} = [\nabla \mathcal{L}(\Theta) \nabla \mathcal{L}(\Theta)^T + \epsilon I_{ND}]^{-1}$, with $I_{ND}$ the $ND \times ND$ identity matrix and $\epsilon \ge 0$ a damping parameter. Using large enough values of $\epsilon$ ensures better numerical stability than Gauss-Newton. For large $\epsilon$ values, this preconditioner boils down to the identity matrix, which is equivalent to the standard GD in Eq. 2.14. For low values, one recovers the Gauss-Newton preconditioner. In practice, the damping parameter $\epsilon$ is often adjusted during the optimization procedure to benefit from the gradient descent stability in the first iterations and from the Gauss-Newton algorithm speed in the last iterations. Similarly to the Newton algorithm, the Gaussian-Newton and Levenberg-Marquardt algorithms involve either the inversion of a dense matrix and matrix-vector products or solving a linear system, which are both expensive in high dimensions.

**Quasi-Newton optimization algorithms** such as BFGS or limited memory BFGS (L-BFGS), propose efficient approximations of the inverse Hessian matrix for more general loss functions (Nocedal and Wright, 2006, chapter 6). In particular, the limited memory BFGS preconditioner does not require any matrix inversion and rely on vector-vector products only. Therefore, it can be used even in high dimensional settings.

**Diagonal preconditioners for neural network training** – Finally, neural networks training involves multiple dedicated diagonal preconditioners – see Chapter 4 for more details on neural
networks. In the deep learning community, these preconditioners are known as “adaptive learning rates”. They are dedicated to non-convex optimization in extremely high dimensions, i.e., $O(10^6)$ and higher. Such preconditioners include the adaptive gradient algorithm (AdaGrad) (Duchi et al., 2011)

$$G^{(t-1)} = \operatorname{diag} \left( \frac{1}{\epsilon + \sqrt{v^{(t-1)}}} \right), \quad \text{with} \quad v^{(t-1)} = \sum_{\tau=1}^{t-1} \left[ \nabla \mathcal{L}(\Theta^{(\tau)}) \right]^2, \qquad (2.20)$$

and the root mean squared propagation (RMSProp) (Tieleman and Hinton, 2012)

$$G^{(t-1)} = \operatorname{diag} \left( \frac{1}{\epsilon + \sqrt{v^{(t-1)}}} \right), \quad \text{with} \quad v^{(t-1)} = a v^{(t-2)} + (1-a) \left[ \nabla \mathcal{L}(\Theta^{(t-1)}) \right]^2, \qquad (2.21)$$

where all exponentiations are taken element-wise, $a \in ]0, 1[$ is an exponential decay rate, and $\epsilon > 0$ is a damping factor that prevents divisions by zero – typically, $\epsilon = 10^{-5}$. The AdaGrad preconditioner accumulates the squared gradients, causing the components of $v$ to diverge. Conversely, RMSProp progressively forgets past gradients with exponential decay rate $\epsilon$. As these two preconditioners are diagonal, they are simple to invert and only involve vector-vector products. In Dauphin et al. (2015), the authors propose another diagonal preconditioner that considerably decreases the upper bound on the condition number from Guggenheimer et al. (1995). They also show that RMSProp is very similar to their preconditioner, which partially explains the performance of RMSProp observed in practice.

As the inverse problem we consider in this thesis is high dimensional, the full rank preconditioners from the Newton, Gauss-Newton and Levenberg-Marquardt algorithms are too computationally expensive to be considered. L-BFGS and RMSProp scale better while enforcing low condition numbers. Section 2.2.2.5 describes how these two preconditioners can be exploited in sampling algorithms.

### 2.2.1.3 Block coordinate variant

In high dimensional settings, evaluating the full gradient at once may be quite costly, both in time and memory. In such cases, one can adopt a divide-and-conquer strategy by dividing the large physical parameter $\Theta$ into $J$ smaller blocks, e.g., $\Theta = (\Theta_j)_{j=1}^J$ (Pereyra et al., 2016). This division replaces the computation of the full gradient vector $\nabla \mathcal{L}$ by the smaller gradient vector $\nabla_j \mathcal{L}$ with respect to part $\Theta_j$. At step $t$, instead of being performed at once, the updated of iterate $\Theta^{(t-1)}$ is divided into a succession of $J$ updates of its parts, with

$$\Theta_j^{(t)} = \Theta_j^{(t-1)} - \eta \nabla_j \mathcal{L}(\Theta_1^{(t)}, \Theta_2^{(t)}, \dots, \Theta_{j-1}^{(t)}, \Theta_j^{(t-1)}, \dots, \Theta_j^{(t-1)}). \qquad (2.22)$$

Considering these faster and lighter individual and alternating updates typically speeds up individual updates and can lighten memory requirements. The updates for each part $j$ can be combined with preconditioning for faster convergence. The update order of the blocks $j$ can be defined either deterministically or stochastically. Besides, a random activation of blocks at each step can lead to a stochastic estimate of the full gradient, yielding a stochastic gradient descent algorithm (Pereyra et al., 2016). Such division of the parameter space typically comes at the price of slower convergence to a critical point in terms of number of iterations, especially in case of degeneracy between components $\Theta_j$. See for instance Chouzenoux et al. (2016) for an application.

#### Limitations of optimization-based inference

**Convergence to the global minimum for non-convex loss function** – Until now, we have described the GD algorithm, its preconditioned variant and its block coordinate variant. At best, these algorithms are guaranteed to converge to a critical point – see e.g., Theorem 10.15 from Beck
(2017). However, the MLE and MAP are defined as global minima of their twice differentiable loss functions $\mathcal{L}$. When $\mathcal{L}$ is in addition convex, any critical point is a global minimum. However, the loss functions considered in this thesis are non-convex. Such functions have potentially many saddle points and local minima with higher loss values than the global minimum. In practice, some non-convex preconditioned variants such as RMSProp escape quickly from the neighborhood of saddle points and reach a local minimum (Dauphin et al., 2015). However, the reached local minimum may have a high loss value. GD algorithms usually fail to escape from such local minima.

A first solution to compensate for this effect is to launch the chosen minimization algorithm many times with different initial states. For high dimensional $\Theta$, this solution requires numerous and potentially costly runs. This is also true in low dimensions when the loss function contains many local minima. Meta-heuristics form an alternative class of algorithms that addresses this multimodality issue. By proposing random candidates, they can escape from local minima with high loss value. In general, they converge to a local minimum with a loss value close to that of the global minimum. See Section 2.A.1 for more details on these methods.

**Lack of uncertainty quantification** – Estimators defined as solutions of an optimization problem do not naturally come with an uncertainty description. Some dedicated methods such as the Cramér-Rao bound propose approximate uncertainty quantification from a point estimate, but rely on strong hypotheses on the posterior, namely unimodality or Gaussianity.

With a point estimate $\hat{\Theta}$, one can derive uncertainties on predictions of reproduced observations $\tilde{Y}|\hat{\Theta}$. These uncertainties can be defined using the noise model from the likelihood in inverse problems. In machine learning, these uncertainties can also be defined without assuming any noise model e.g., with quantile regression (Koenker, 2005) or conformal prediction (Angelopoulos and Bates, 2022). We do not cover these distribution-free uncertainty quantification approaches in this thesis, as our goal is to provide uncertainties on the physical parameters $\Theta$.

### 2.2.2 Approximating integrals with Monte Carlo estimators

As we showed in Section 2.1.2, some estimators can accurately quantify the uncertainty on physical parameters $\Theta$. These estimators are defined as an integral of a function $g$ on $\Theta$ over the posterior distribution. In this section, we first introduce Monte Carlo (MC) estimators that approximate the integral using samples of the posterior distribution. Then, we present Markov chain Monte Carlo (MCMC) algorithms, including the widespread Metropolis-Hastings (MH) algorithm and Metropolis adjusted Langevin algorithm (MALA). Preconditioned and block coordinate variants are also presented. Finally, the multiple-try Metropolis (MTM), which extends the MH algorithm by generating multiple candidates at each iteration instead of one, is presented.

The sampler proposed in Chapter 5 (Section 5.2) builds on the algorithms presented in this section. It combines two sampling kernels. The first addresses the log-posterior gradient Lipschitz regularity issue. It builds on preconditioning to explore the parameter space more efficiently. The second addresses the posterior multimodality difficulty. It performs block coordinate updates and generates multiple candidates for each block to improve mixing properties.

Additional sampling algorithms dedicated to multimodal distributions are presented in Appendix 2.A. Three classes of algorithms are considered: adaptations of meta-heuristics to MCMC algorithms, exploitation of prior knowledge of the modes localization, and nested sampling. We do not use these algorithms to derive the proposed sampler. However, they are popular in interstellar medium studies, as we will show in Chapter 3 (Section 3.2.2). Besides, we will compare these algorithms to the proposed sampler in Chapter 5 (Section 5.4).
### 2.2.2.1 Monte Carlo (MC) estimators

In case of no conjugacy relation, two main approaches permit to numerically evaluate the integral in Eq. 2.11, reminded here,

$$ \mathbb{E}[g(\Theta)|\mathbf{Y}] = \int g(\Theta) \pi(\Theta|\mathbf{Y}) d\Theta. \qquad (2.23) $$

The first approach relies on a deterministic quadrature of the parameter space, i.e., on a finite set of points with deterministic positions. The posterior pdf is then evaluated on each point. The integration is performed using Riemannian integration, or more accurate methods such as the trapezoidal rule or Simpson's rule. This approach is efficient in low dimensions (Robert and Casella, 2004, chapter 1). However, the number of points in the quadrature increases exponentially with the dimension. This method has already been applied in ISM studies, as detailed in Chapter 3 (Section 3.2.2).

The second approach relies on a stochastic quadrature of the parameter space. It first generates samples $\Theta^{(t)}$ from the posterior distribution. Estimators are then evaluated as empirical means of the $g(\Theta^{(t)})$:

$$ \mathbb{E}[g(\Theta)|\mathbf{Y}] \approx \frac{1}{T_{\text{MC}}} \sum_{t=1}^{T_{\text{MC}}} g(\Theta^{(t)}), \quad \Theta^{(t)} \sim \pi(\Theta|\mathbf{Y}), \quad t=1, \dots, T_{\text{MC}}. \qquad (2.24) $$

Such estimators are called Monte Carlo (MC) estimators (Robert and Casella, 2004, chapter 3) – hence the notation “MC” for $T_{\text{MC}}$.

For independent and identically distributed (i.i.d.) samples drawn from the posterior, the law of large numbers guarantees that

$$ \lim_{T_{\text{MC}} \to \infty} \frac{1}{T_{\text{MC}}} \sum_{t=1}^{T_{\text{MC}}} g(\Theta^{(t)}) = \mathbb{E}[g(\Theta)|\mathbf{Y}], \qquad (2.25) $$

and the central limit theorem (CLT) describes the distribution on the error of this empirical mean compared to the true one

$$ \left( \frac{1}{T_{\text{MC}}} \sum_{t=1}^{T_{\text{MC}}} g(\Theta^{(t)}) - \mathbb{E}[g(\Theta)|\mathbf{Y}] \right) \xrightarrow{D} \mathcal{N} \left( 0, \frac{1}{T_{\text{MC}}} \text{Cov}[g(\Theta)|\mathbf{Y}] \right), \qquad (2.26) $$

where the convergence is in distribution and $\text{Cov}[g(\Theta)|\mathbf{Y}]$ is a fixed covariance matrix, generally unknown. Essentially, the CLT shows that an MC estimator built with i.i.d. samples converges with speed $1/\sqrt{T_{\text{MC}}}$. The variances in $\text{Cov}[g(\Theta)|\mathbf{Y}]$ scale polynomially with the dimension of the parameter space. The sampling approach therefore scales much better than deterministic integration, at the condition of efficient sampling of the posterior distribution.

### 2.2.2.2 From Monte Carlo to Markov chain Monte Carlo (MCMC) methods

There exists no efficient algorithm that draws independent samples from an arbitrary posterior distribution. Rejection sampling (Robert and Casella, 2004, chapter 2) can in principle draw independent samples from any posterior distribution. It generates candidates according to a proposal distribution and accepts or rejects them with a certain acceptance probability. In practice, the acceptance probabilities are prohibitively small in high dimensional settings with a complex model.

Some alternative methods generate correlated samples from which MC estimators can be evaluated. Despite the correlation between samples, MC estimators built from Markov chains enjoy similar convergence results as MC estimators built from independent samples. Such algorithms
include nested sampling, sequential MC (SMC) and Markov chain Monte Carlo (MCMC) methods. Nested sampling and SMC are widespread in astrophysics, but are not used directly in this thesis. They are presented in Appendix 2.A. They both build on or rely on MCMC methods.

An MCMC algorithm (Robert and Casella, 2004, chapters 6 and 7) generates a Markov chain whose stationary distribution is the posterior distribution. The goal is to generate a correlated set of samples from the posterior distribution to evaluate MC estimators (Eq. 2.24). To do so, each MCMC algorithm relies on a different transition kernel $\mathcal{K} : \mathbb{R}^{N \times D} \times \mathbb{R}^{N \times D} \to \mathbb{R}_+$. The MCMC class gathers a wide variety of algorithms that can tackle problems with diverse sets of properties.

**Burn-in phase** – A Markov chain generated with an MCMC algorithm contains two regimes. The chain is initialized with a starting point $\Theta^{(0)}$. At first, it is in a transient regime, usually called the *burn-in phase*, where iterates are close to $\Theta^{(0)}$ and evolve towards high probability regions. When it reaches a high probability region, the Markov chain enters a stationary regime, in which the elements $\Theta^{(t)}$ of the chain can be considered as samples of the posterior. In particular, these samples are considered independent of $\Theta^{(0)}$. The MC estimators are evaluated with these samples. The first $T_{BI} \ge 1$ iterates associated with the burn-in phase are removed to avoid biasing estimators. A good transition kernel $\mathcal{K}$ thus reaches quickly the posterior high probability region and explores it efficiently.

**Convergence properties**$^4$ – The convergence of the Markov chain to the posterior distribution and of the associated MC estimators to the true values are difficult to establish. Three properties are central. First, a Markov chain admits an *invariant distribution*. An invariant distribution $\pi$ verifies

$$ \left\{ \begin{array}{l} \Theta^{(1)} \sim \pi \\ \Theta^{(2)} \sim \mathcal{K}(\Theta^{(1)}, \cdot) \end{array} \right. \implies \Theta^{(2)} \sim \pi. \qquad (2.27) $$

The *detailed balance* is a sufficient condition to establish that the posterior distribution is the invariant distribution of a Markov chain. A transition kernel $\mathcal{K}$ verifies the detailed balance property with the posterior if and only if

$$ \forall \Theta^{(1)}, \Theta^{(2)} \in \mathbb{R}^{N \times D}, \quad \mathcal{K}(\Theta^{(1)}, \Theta^{(2)}) \pi(\Theta^{(1)}|\mathbf{Y}) = \mathcal{K}(\Theta^{(2)}, \Theta^{(1)}) \pi(\Theta^{(2)}|\mathbf{Y}) \quad (2.28) $$

A Markov chain with such a transition kernel is reversible and admits the posterior as invariant density (Robert and Casella, 2004, Theorem 6.46).

Second, the *ergodicity* of a Markov chain guarantees the asymptotic convergence of the Markov chain to its invariant with respect to a distance on probability distributions, the total variation. In particular, this result holds for any starting point $\Theta^{(0)}$. It also guarantees that MC estimators asymptotically converge to their true values, i.e., that the law of large numbers (Eq. 2.25) applies. The *ergodicity theorem* (Robert and Casella, 2004, Theorem 6.63) states that a sufficient condition for ergodicity is for the chain to be Harris recurrent. We now provide an intuition of the definition of Harris recurrence – see Robert and Casella (2004, Definition 6.32) for the exact definition. The *irreducibility* of a Markov chain means that its transition kernel $\mathcal{K}$ allows for a transition from any $\Theta^{(1)} \in \mathbb{R}^{N \times D}$ to any $\Theta^{(2)} \in \mathbb{R}^{N \times D}$, potentially with multiple steps. Besides, for any subset of $\mathbb{R}^{N \times D}$ of non-zero measure, if the transition kernel $\mathcal{K}$ permits to remain in the subset with non-zero probability, the chain is said to be *aperiodic*. The *recurrence* of a Markov chain means that the expected number of visits of any subset of $\mathbb{R}^{N \times D}$ of non-zero measure is infinite. From Robert and Casella (2004, Theorem 6.30), an irreducible chain is either recurrent or transient. The *Harris recurrence* is a stronger recurrence property for an irreducible chain in which the probability of visiting any subset of $\mathbb{R}^{N \times D}$ of non-zero measure an infinite number of times is exactly 1.

<footnote n="4">4. This paragraph contains mathematical details that are not mandatory for the remainder of the thesis.</footnote>
Third, the speed of convergence of a Markov chain to its invariant distribution permits the derivation of stopping rules for required precision levels. There are two main characterizations of the speed of convergence for Markov chains. The first is geometric ergodicity (Robert and Casella, 2004, Definition 6.54). In essence, it ensures that an extended total variation distance between the Markov chain and its invariant distribution decreases at least geometrically. The second, uniform ergodicity (Robert and Casella, 2004, Definition 6.58), enforces a stronger condition on the rate of geometric convergence. Both conditions can be associated with a CLT (Eq. 2.26) - see Robert and Casella (2004, Theorems 6.67 and 6.77). A third condition permits to derive a CLT requires the Markov chain to be aperiodic, irreducible and reversible with the posterior as invariant distribution.

**CLT and effective sample size (ESS)** – In the context of Markov chains, the three previous CLTs yield, for a scalar component *i* of function *g*,

$$ \left( \frac{1}{T_{\text{MC}}} \sum_{t=1}^{T_{\text{MC}}} g_t(\Theta^{(t)}) - \mathbb{E}[g_t(\Theta)|\mathbf{Y}] \right) \rightarrow \mathcal{N} \left( 0, \frac{\gamma_{g_t}^2}{T_{\text{MC}}} \right), \quad (2.29) $$

with

$$ \frac{\gamma_{g_t}^2}{T_{\text{MC}}} = \frac{1}{T_{\text{MC}}} \left( \text{Var}[g_t(\Theta)|\mathbf{Y}] + 2 \sum_{\tau=1}^{\infty} \text{Cov}[g_t(\Theta^{(t)}), g_t(\Theta^{(t+\tau)})|\mathbf{Y}] \right). \quad (2.30) $$

The first element in the definition of $\gamma_{g_t}^2$ corresponds to the variance in the CLT for i.i.d. samples (Eq. 2.26). The second term sums the autocovariances in the Markov chain with lag $\tau \ge 1$. These autocovariances can be converted into more interpretable autocorrelations,

$$ \frac{\gamma_{g_t}^2}{T_{\text{MC}}} = \frac{1}{T_{\text{MC}}} \text{Var}[g_t(\Theta)|\mathbf{Y}] \left( 1 + 2 \sum_{\tau=1}^{\infty} r_t^{(\tau)} \right), \quad (2.31) $$

where $r_t^{(\tau)}$ is the autocorrelation of the chain with time lag $\tau \ge 0$ for $g_t(\Theta)$. Finally, one can rewrite Eq. 2.31

$$ \frac{\gamma_{g_t}^2}{T_{\text{MC}}} = \frac{\text{Var}[g_t(\Theta)|\mathbf{Y}]}{\text{ESS}^{(T_{\text{MC}})}}, \quad (2.32) $$

where the ESS stands for the effective sample size (ESS) (Robert and Casella, 2004, Section 12.3.5). It permits to rewrite the CLT in Eq. 2.26 in the context of Markov with a size ESS($T_{\text{MC}}$), instead of $T_{\text{MC}}$. Therefore, the ESS measures the number of independent samples with the same estimation power as a correlated chain. Ideally, the correlation between iterates $\Theta^{(t)}$ should be as small as possible. In the extreme case where the correlation in the chain is negligible, ESS($T_{\text{MC}}$) $\simeq T_{\text{MC}}$, i.e., the Markov chain can be considered as a set of independent samples. Conversely, highly correlated chains typically “mix” poorly, as they need a large amount of iterations to explore the posterior.

In the general definition above, the ESS depends on the function $g$. In practice, for simplicity, it is often reduced to the case $g: \Theta \mapsto \Theta$, i.e.,

$$ \text{ESS}_{nd}^{(T_{\text{MC}})} = \frac{T_{\text{MC}}}{1 + 2 \sum_{\tau=1}^{\infty} r_{nd}^{(\tau)}}, \quad (2.33) $$

with $r_{nd}^{(\tau)}$ the autocorrelation of the chain with time lag $\tau \ge 0$ for the parameter $\theta_{nd}$. As the above ESS is defined for a scalar parameter, for a *ND*-dimensional parameter space, one obtains *ND* estimations of the ESS. In this thesis, we implement the ESS estimation presented in Gelman et al. (2015). This definition accounts for results of potentially *M* $\ge 1$ Markov chains. It can therefore be used when an MCMC algorithm is run *M* times.
### 2.2.2.3 Metropolis-Hastings (MH) algorithm

The Metropolis-Hastings (MH) algorithm (Metropolis et al., 1953; Hastings, 1970) is historically the first proposed MCMC algorithm, and remains arguably the most famous and widespread to this day. It is very general and versatile, and enjoys strong theoretical properties. The two kernels proposed in Chapter 5 (Section 5.2) build on this algorithm.

**The algorithm** – The MH algorithm can start from any point $\Theta^{(0)}$. Then, at each step $t$, a candidate $\Theta_c^{(t)}$ is sampled from a proposal distribution $q(\Theta_c|\Theta^{(t-1)})$. Unlike rejection sampling, the proposal distribution can depend on the current iterate $\Theta^{(t-1)}$. This candidate is accepted with probability

$$ \rho^{(t)} = \min \left( 1, \frac{\pi(\Theta_c^{(t)}|\mathbf{Y}) q(\Theta^{(t-1)}|\Theta_c^{(t)})}{\pi(\Theta^{(t-1)}|\mathbf{Y}) q(\Theta_c^{(t)}|\Theta^{(t-1)})} \right). \quad (2.34) $$

When the candidate is not accepted, the new iterate $\Theta^{(t)}$ is set to $\Theta^{(t-1)}$. Therefore, the transition kernel $K$ of the MH algorithm is defined as a combination of the proposal $q$ and the accept-reject step. Algorithm 2.2 summarizes the MH algorithm.

<caption>**Algorithm 2.2: Metropolis-Hastings (MH) algorithm**</caption>

**Input:** Starting point $\Theta^{(0)}$, proposal $q$, numbers of samples $T_{BI}$ and $T_{MC}$
**for** $t = 1, \dots, T_{MC}$ **do**
2 &nbsp;&nbsp;&nbsp;&nbsp; $\Theta_c^{(t)} \sim q(\Theta_c|\Theta^{(t-1)})$ &nbsp;&nbsp;&nbsp;&nbsp; // generate candidate
3 &nbsp;&nbsp;&nbsp;&nbsp; $\rho^{(t)}$ &nbsp;&nbsp;&nbsp;&nbsp; // compute acceptance probability using Eq. 2.34
4 &nbsp;&nbsp;&nbsp;&nbsp; $\zeta \sim \text{Unif}(0, 1)$ &nbsp;&nbsp;&nbsp;&nbsp; // update iterate
5 &nbsp;&nbsp;&nbsp;&nbsp; $\Theta^{(t)} = \Theta_c^{(t)}$ if $\zeta \le \rho^{(t)}$ else $\Theta^{(t-1)}$
**Output:** Chain of samples $\{\Theta_c^{(t)}\}_{t=T_{BI}+1}^{T_{MC}}$ &nbsp;&nbsp;&nbsp;&nbsp; // The first $T_{BI}$ are rejected

**Theoretical convergence properties** – The theoretical convergence properties of the MH algorithm are simple. The following is a short summary of Robert and Casella (2004, Sections 7.3.1 and 7.3.2). First, the accept-reject step allows $\Theta^{(t)} = \Theta^{(t-1)}$. This condition is sufficient for the Markov chain to be aperiodic. Then, the main (sufficient) condition for convergence relies on the positivity of the proposal $q$ on $\mathbb{R}^{N \times D}$, i.e.,

$$ \forall \Theta_c, \Theta^{(t-1)} \in \mathbb{R}^{N \times D}, \quad q(\Theta_c|\Theta^{(t-1)}) > 0. \quad (2.35) $$

This condition is verified e.g., for a Gaussian proposal with any mean vector and covariance matrix. When it is verified, the transition kernel satisfies the detailed balance property (Eq. 2.28). Therefore, the Markov chain generated by the MH algorithm admits the posterior distribution as an invariant probability distribution (Robert and Casella, 2004, Theorem 7.2). Besides, it is also a sufficient condition for the chain ($\Theta^{(t)}$) to be irreducible. Therefore, with such a proposal $q$, the law of large numbers (Eq. 2.25) applies to the MH Markov chain (Robert and Casella, 2004, Theorem 7.4). Note that this condition is only sufficient. Similar results hold for less restrictive conditions (Robert and Casella, 2004, Corollary 7.7).

**Common proposal distributions** $q$ – As the theoretical results show, the choice of the proposal distribution $q$ is crucial to efficiently explore the posterior distribution.

Choosing a proposal $q$ that is independent of current iterate $\Theta^{(t-1)}$, i.e., $q(\cdot|\Theta^{(t-1)}) = q(\cdot)$, seems a simple and natural choice. Such a proposal yields the so-called independent Metropolis-Hastings (I-MH). As stated in Robert and Casella (2004, Theorem 7.8), this algorithm produces a uniformly ergodic chain if there exists a constant $M \ge 1$ such that $Mq(\Theta) \ge \pi(\Theta|\mathbf{Y})$ for all
$\Theta \in \mathbb{R}^{N \times D}$. In this case, the expected probability acceptance is at least $1/M$ when the chain is stationary. As for rejection sampling, expected probability acceptance may thus be extremely low in practice for posteriors whose high probability region are concentrated in a region with much smaller volume than that of the proposal. This is typically the case in high dimensions, due to the curse of dimensionality.

Random walk Metropolis-Hastings (RWMH) is another very simple proposal that is more widely spread. It generates candidates $\Theta_c^{(t)} = \Theta^{(t-1)} + Z^{(t)}$, with generally $Z^{(t)} \sim N(0, \Sigma)$ for a covariance matrix $\Sigma$. In this case, $q(\cdot|\Theta^{(t-1)}) = N(\Theta^{(t-1)}, \Sigma)$. As this proposal $q$ is symmetric, i.e., $q(\Theta_c^{(t)}|\Theta^{(t-1)}) = q(\Theta^{(t-1)}|\Theta_c^{(t)})$, the acceptance probability in Eq. 2.34 simplifies to the ratio of posterior pdfs of the candidate and of the current iterate. RWMH is never uniformly ergodic. It can be shown to be geometrically ergodic for log-concave distributions (Robert and Casella, 2004, Theorem 7.15). However, as the Markov chain is reversible, the CLT still applies. In practice, the covariance matrix $\Sigma$ plays a crucial role in RWMH: a matrix $\Sigma$ closer to the covariance of the posterior will lead to larger steps with higher acceptance probabilities. However, even with an adapted covariance matrix, RWMH does not scale up well due to its blind nature (Pereyra et al., 2016).

As RWMH performs local steps, its optimal expected acceptance probability is not 1. Indeed, a covariance matrix $\Sigma$ with small variances will generate candidates very close to the current iterate. In case of smooth log-posterior, these candidates will have an expected acceptance probability close to 1, but will lead to slow exploration of the parameter space. Conversely, if the variances are large, the candidates may fall far from the current iterate, potentially out of the posterior high probability region. In high dimensions, such candidates have an expected acceptance probability close to 0. However, each accepted candidate corresponds to a potentially large step, which improves the exploration. In general, the optimal mixing for RWMH is achieved with a trade-off between large steps and large acceptance probability. Optimal mixing for RWMH is achieved with an expected acceptance probability of 20% - 25% (Pereyra et al., 2016).

### 2.2.2.4 Using the gradient with Metropolis adjusted Langevin algorithm (MALA)

Metropolis adjusted Langevin algorithm (MALA) (Roberts and Stramer, 2002) and Hamiltonian Monte Carlo (HMC) (Neal, 2011) are both MH algorithms that exploit gradient information in the mean of the proposal $q$. Thanks to the gradient information, they both scale better than RWMH (Pereyra et al., 2016).

HMC relies on the introduction of a momentum auxiliary variable and exploits Hamiltonian dynamics to explore the posterior. MALA is similar to a gradient descent algorithm on the negative log-posterior. HMC generally mixes better than MALA for a fixed number of samples. However, HMC requires tuning more parameters and has a higher computational cost per iteration. Besides, MALA can be easily associated with advanced optimization methods (Pereyra et al., 2016). It is therefore simpler to transfer advanced optimization techniques to MALA than to HMC. For these reasons, we favor MALA over HMC in this thesis.

MALA can be defined from a particular Langevin diffusion process. This continuous diffusion process admits the posterior as an invariant distribution and is geometrically ergodic. First, the diffusion process is discretized, usually following the Euler-Maruyama scheme. The resulting random process is called unadjusted Langevin algorithm (ULA), and reads

$$ \Theta_c^{(t)} = \Theta^{(t-1)} + \eta \nabla \ln \pi(\Theta^{(t-1)}|\mathbf{Y}) + \sqrt{2\eta} Z^{(t)}, \quad Z^{(t)} \sim N(0_{ND}, \mathbf{I}_{ND}). \quad (2.36) $$

ULA is a GD algorithm perturbed with an additive Gaussian noise. It therefore enjoys some similar theoretical properties. The discretization of the diffusion process introduces an error in ULA. Due to this error causes, the posterior is generally not the invariant distribution of ULA. In Durmus and Moulines (2017), the authors propose three upper bounds on the distance between an ULA Markov chain and the posterior distribution. Each upper bound is associated with a
class of posterior distribution: super-exponential outside a ball around the mode – i.e., posteriors with heavy tails –, log-concave, and strongly log-concave. In all three cases, the log-posterior is assumed to be gradient Lipschitz continuous with Lipschitz constant $\beta$. The three upper bounds depend on the Lipschitz constant $\beta$, either directly or indirectly through the step size $\eta \in ]0, 1/\beta[$. Similarly to GD in optimization, step sizes close to $1/\beta$ lead to faster convergence to the invariant distribution.

MALA corrects the discretization error in ULA thanks to an accept-reject step. This accept-reject steps makes MALA an MH algorithm. Its proposal distribution can thus be written

$$q(\cdot|\Theta^{(t-1)}) = \mathcal{N}(\Theta^{(t-1)} + \eta \nabla \ln \pi(\Theta^{(t-1)}|\mathbf{Y}), 2\eta \mathbf{I}_{ND}), \quad (2.37)$$

which is an equivalent form of Eq. 2.36. Similarly to the covariance matrix in RWMH, the step size $\eta$ defines a trade-off between large steps and high expected acceptance probability. As it exploits gradient information, MALA produces better candidates than RWMH. In case of a gradient Lipschitz continuous log-posterior of Lipschitz constant $\beta$, a step size $\eta$ close to $1/\beta$ leads to good mixing properties. Overall, the MALA optimal expected acceptance probability is 50% - 60% (Pereyra et al., 2016), i.e., two to three times larger than RWMH.

**Discretization error** – As already mentioned, the Euler-Maruyama scheme introduces a discretization error. In MALA, this error is corrected with the accept-reject step. However, in case of no discretization error, all candidates would be accepted, which would accelerate the posterior exploration. Therefore, a lower discretization error leads to more relevant candidates. For instance, in Durmus et al. (2017), the authors exploit an Ozaki discretization scheme, which is more accurate than the usual Euler-Maruyama. They obtain considerable improvement of the mixing performance of MALA, and derive an optimal expected acceptance probability around 70%. In general, using a better discretization scheme reduces this discretization error and therefore generates better candidates.

### 2.2.2.5 Preconditioning with MALA to improve mixing properties

As stated in Section 2.2.1 for optimization methods, preconditioning permits to reduce the condition number of the Hessian matrix. This reduction permits a more efficient exploration of the parameter space. As ULA is defined as a GD algorithm perturbed with Gaussian noise, it seems natural to extend it with preconditioning.

**Preconditioning and Riemannian manifold** – In Girolami and Calderhead (2011), the authors associate preconditioning with sampling on a Riemannian manifold. The preconditioner is then closely related to the metric of the manifold. The authors extended MALA and HMC to Riemannian manifolds, yielding the Manifold Metropolis adjusted Langevin algorithm (MMALA) and Riemannian manifold HMC (RMHMC) algorithms. In Xifara et al. (2014), a minor correction was applied to MMALA, leading to a renamed preconditioned Metropolis adjusted Langevin algorithm (PMALA). We use PMALA to derive one of the two sampling kernels of the MCMC algorithm proposed in Chapter 5 (Section 5.2). The associated proposal distribution $q$ reads

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

with

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

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)}}. \qquad (2.40) $$

Depending on the complexity of the preconditioner G, this additional drift term can be quite complicated to evaluate. For instance, for a G defined with second order derivative information, evaluating this additional drift term will require accessing third order derivatives.

**Stochastic gradient MCMC (SG-MCMC) and preconditioner choice** – SG-MCMC (Welling and Teh, 2011) is a special class of MCMC algorithms that makes an extensive use of preconditioning and of the results from Girolami and Calderhead (2011). They extend MALA and HMC for stochastic estimates of the gradient. As in optimization, the stochasticity is generally included to avoid a costly evaluation of the full gradient. A first example of SGD algorithm was presented in Section 2.2.1.3. A second example of application for these methods is to sample a posterior distribution on the parameters of a neural network – see Chapter 4 for more details on neural networks. Such posterior distributions are generally not log-concave and very high-dimensional – $O(10^6)$. Like in optimization methods for neural network training, preconditioning is necessary in this context.

We consider three main preconditioners proposed in SG-MCMC for a stochastic gradient MALA: the expected Fisher information matrix (Patterson and Teh, 2013), L-BFGS (Simsekli et al., 2016) and RMSProp (Li et al., 2016).

Using the full expected Fisher information matrix leads in general to a dense covariance matrix, i.e., with mostly non-zero terms. In high dimensions such as $O(10^3)$ or $O(10^4)$, inverting such a matrix is prohibitively expensive.

The L-BFGS preconditioner and its inverse are defined with sequences of vectors and inner products. They actually never require to compute a matrix or a matrix-vector product. This preconditioner thus scales very well. In addition, in Simsekli et al. (2016), the authors slightly modify the definition of L-BFGS to ensure $\gamma^{(t-1)} = 0$ for all $t \ge 1$. However, this preconditioner was designed for log-concave posterior distribution. A positive regularization parameter is proposed to ensure positive definiteness of the preconditioner, but this parameter can become extremely large for highly non log-concave distributions. In the limit case of a very high regularization parameter, the preconditioner boils down to the identity matrix. Therefore, although this preconditioner is able to capture correlations between parameters, it is not applicable in the inverse problem considered in this thesis.

Finally, RMSProp is a diagonal preconditioner and thus is trivial to invert. As we mentioned in Section 2.2.1.2, this preconditioner shows good performance in practice. This performance in practice was partially explained in Dauphin et al. (2015). For this reason, in Chapter 5 (Section 5.2.1), we propose a PMALA sampling kernel equipped with the RMSProp preconditioner. Note that RMSProp being diagonal, it is not efficient in case of degeneracies, i.e., of high correlation between two parameters. Other preconditioners with particular structure such as a band matrix or block diagonal may still be used to address a limited number of degeneracies. However, as we will show in Chapter 5, RMSProp yields a drift term $\gamma$ that can be evaluated with limited additional cost.

### 2.2.2.6 Gibbs sampling: a block coordinate MH variant

One of the current main challenges in MCMC is sampling from high dimensional posterior distributions, e.g., distributions on images. In Section 2.2.1.3, we showed that block coordinate optimization divides a high dimensional problem into smaller ones. In images, each block $j$ can correspond to one pixel or to a group of pixels. Similarly, Gibbs sampling (Geman and Geman, 1984) is an MCMC algorithm that performs component-wise updates.
**Gibbs and Metropolis-within-Gibbs algorithms** – In Gibbs sampling, at each step $t$, the components $j = 1, \dots, J$ are updated sequentially. Each component $j$ is updated by sampling from the conditional posterior

$$ \Theta_j^{(t)} \sim \pi \left( \Theta_j | \mathbf{Y}, \Theta_{\setminus j}^{(t-1+\frac{j-1}{J})} \right), \qquad (2.41) $$

with $\Theta_{\setminus j}^{(t-1+\frac{j-1}{J})} = (\Theta_1^{(t)}, \dots, \Theta_{j-1}^{(t)}, \Theta_{j+1}^{(t-1)}, \dots, \Theta_J^{(t-1)})$. The $\Theta_{\setminus j}$ notation corresponds to $\Theta$ with the $j$th block removed. As shown in Geman and Geman (1984), the Gibbs sampler leads to ergodic Markov chains for mild assumptions on the conditional posterior.

As we already discussed, sampling directly from the conditional posterior may be out of reach. The so-called Metropolis-within-Gibbs algorithm (Gelman et al., 2015, chapter 11) resorts to MH steps for each of the $J$ components, i.e., using proposals $q_j$ and performing an accept-reject step. Algorithm 2.3 summarizes the Metropolis-within-Gibbs algorithm. Both Gibbs sampling and Metropolis-within-Gibbs can produce geometric and uniform ergodic chains (Johnson et al., 2013). At each step $t$, Gibbs sampling and Metropolis-within-Gibbs generate $J$ samples from conditional posterior distributions with lower dimension, instead of generating one sample from the full posterior distribution. Such a divide-and-conquer approach can be very useful for distributions that show a particular structure, such as images or time series.

<caption>**Algorithm 2.3: Metropolis-within-Gibbs sampling**</caption>

**Input:** Starting point $\Theta^{(0)}$, number of samples $T_{MC}$
**for** $t = 1, \dots, T_{MC}$ **do**
1 &nbsp;&nbsp;&nbsp;&nbsp;**for** $j = 1, \dots, J$ **do**
2 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**if** $\Theta_j^{(t)} \sim q_j (\Theta_j | \mathbf{Y}, \Theta_{\setminus j}^{(t-1+\frac{j-1}{J})})$ **then**
3 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**Generate candidate for component j**
4 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**Compute acceptance probability using Eq. 2.34**
5 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**Update iterate**
6 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end if**
7 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;**end for**
**Output:** Chain of samples $\{\Theta^{(t)}\}_{t=1}^{T_{MC}}$

**Chromatic Gibbs sampling** – One constraint of the Gibbs and Metropolis-within-Gibbs algorithms is the sequential sampling of each component $n$. This sequential sampling is necessary to take into account the latest value of components $n_0$ and $n_1$ to update some other $n_2 \neq n_1$. However, in case of conditional independence between $\theta_{n_1}$ and $\theta_{n_2}$ when $\theta_{n_0}$ is given, then $\theta_{n_1}$ and $\theta_{n_2}$ can be sampled in parallel. In general, the set of conditional dependencies between individual components $\theta_n$ is represented as a graph. This graph can be colored such that two neighbor nodes never share the same color. Each color then corresponds to a set of components $n$ that can be sampled in parallel using the chromatic Gibbs algorithm (Gonzalez et al., 2011). This parallelization, combined with vectorized computations, can considerably accelerate the sampling. Besides, this algorithm also produces ergodic Markov chains.

Figure 2.2 shows two examples of graphs of conditional dependence. The graph on the left shows the extreme case of a fully connected graph. Such a graph requires one color per node. In general, finding the minimum number of colors for a given graph is a NP-hard problem. However, in many practical use cases, one can easily find the smallest number of colors. The graph on the right shows a graph of pixels in an image. In this graph, a pixel only interacts with its direct neighbors. In this case, it is easy to see that the graph can be partitioned into two colors. Therefore, in this case, at each step $t$, only two steps are required, whatever the size of the image. In Chapter 5 (Section 5.1.2), we set a spatial regularization prior on maps of physical parameters that yields this exact graph. The sampling kernel that we propose in Chapter 5 (Section 5.2.2) thus exploits chromatic Gibbs sampling.
<figure><img src="image_18.png" /><figcaption>Figure 2.2 – Graphs of conditional dependence and associated optimal coloring.</figcaption></figure>

### 2.2.2.7 Generating multiple candidates

As mentioned in Sections 2.2.2.3 and 2.2.2.4, the acceptance rate that leads to the best mixing is roughly 20%-25% for RWMH and 50%-60% for MALA. A proposal distribution that both achieves such acceptance rates and takes relatively large steps can be out of reach, especially in high dimensions. In such cases, very long chains are necessary to fully explore the posterior. An alternative to taking long chains is to generate multiple candidates at each step, which leads to an increase of the acceptance probability. This idea is at the core of multiple-try Metropolis (MTM) algorithms (Liu et al., 2000).

At each step $t$, MTM algorithms first generate $K \ge 1$ candidates $\Theta_c^{(k)}$ instead of 1 in MH algorithms. Then, using an importance weight function

$$w(\Theta_c^{(k)}) = \frac{\pi(\Theta_c^{(k)} | \mathbf{Y}, \Theta^{(t-1)})}{q(\Theta_c^{(k)} | \Theta^{(t-1)})}, \qquad (2.42)$$

one candidate $i$ is selected according to the set of selection probabilities $(w_k)_{k=1}^K$ with a multinomial distribution

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

The accept-reject step is then performed with the selected candidate $i$ and the generalized acceptance probability (Liu et al., 2021; Martino, 2018)

$$\tilde{\rho}^{(t)} = \min \left( 1, \frac{w(\Theta_c^{(i)}) + \sum_{j=1, j \neq i}^{K} w(\Theta_c^{(j)})}{w(\Theta^{(t-1)}) + \sum_{j=1, j \neq i}^{K} w(\Theta_c^{(j)})} \right). \qquad (2.44)$$

Algorithm 2.4 summarizes the MTM sampler. Note that for $K = 1$, MTM becomes equivalent to MH. Like MH, MTM verifies the detailed balance property and thus admits the posterior as invariant distribution. In addition, like MH, MTM produces ergodic Markov chains. Finally, like I-MH, proposal distributions $q$ that are independent of the current iterate $\Theta^{(t-1)}$ – independent multiple-try Metropolis (I-MTM) – leads to uniformly ergodic Markov chains if there exists a finite constant $M$ such that $Mq(\Theta) \ge \pi(\Theta|\mathbf{Y})$ for all $\Theta$. In particular, such proposal distributions do not get trapped in local modes. See Martino (2018) for a review on MTM algorithms.