Estimating hyperparameters and instrument parameters in regularized inversion. Illustration for SPIRE/Herschel map making

Résumé Technique : Estimation des hyperparamètres et des paramètres d'instrument dans l'inversion régularisée

Énoncé du problème
L'article traite du défi de la reconstruction d'images en astrophysique, plus précisément pour l'instrument SPIRE de l'observatoire spatial Herschel. Le problème central est un problème d'inversion « non supervisé et myope ». Les méthodes d'inversion standard nécessitent généralement :

1. Un réglage supervisé : Les hyperparamètres (qui équilibrent la fidélité aux données et la régularisation a priori) doivent être définis manuellement, ce qui nécessite souvent la connaissance du « ciel vrai » ou de nombreux essais et erreurs.
2. Des modèles d'instrument connus : Les paramètres de l'instrument (tels que la largeur de la tache, le temps de réponse et le gain) sont habituellement déterminés via des observations de calibration dédiées. Ces valeurs contiennent inévitablement des erreurs dues aux dépendances spectrales et aux incertitudes de calibration.

Les auteurs visent à résoudre un problème tripartite à partir d'une seule observation : estimer simultanément la carte du ciel ($x$), les hyperparamètres ($\xi$) et les paramètres de l'instrument ($\eta$) sans connaissance préalable de la vérité terrain ou de données de calibration dédiées.

Méthodologie
L'approche proposée est fondée sur un cadre bayésien complet. Les auteurs construisent une densité de probabilité postérieure $p(x, \xi, \eta | y)$ pour toutes les quantités inconnues données les observations $y$.

1. Modélisation probabiliste :

   • Vraisemblance : Le modèle d'observation est linéaire ($y = A_\eta x + n$), où le bruit $n$ est modélisé par un processus gaussien avec une moyenne $m_n$ et une structure de covariance définie par un paramètre de précision $\gamma_n$.
   • A priori pour la carte : La carte du ciel $x$ est modélisée par un a priori gaussien pour imposer une régularité spatiale (lissage), contrôlée par un paramètre de précision $\gamma_x$. La structure de régularité est définie via une norme quadratique des dérivées spatiales.
   • Hyperpriorités : Les hyperparamètres ($\xi = [m_n, \gamma_n, \gamma_x]$) se voient assigner des a priori de Jeffreys non informatifs ($1/\gamma$ pour les précisions, uniforme pour les moyennes) afin de minimiser le biais a priori.
   • Paramètres de l'instrument : Les paramètres de l'instrument $\eta$ sont modélisés avec des a priori gaussiens centrés sur des valeurs nominales avec des incertitudes associées.

2. Densité postérieure :
   La densité postérieure jointe est dérivée en combinant la vraisemblance et les a priori. En raison de la haute dimensionnalité de la carte ($M \sim 10^5$) et de la dépendance complexe et non linéaire des observations vis-à-vis des paramètres de l'instrument $\eta$, la postérieure ne peut être analysée directement ou maximisée analytiquement.

3. Échantillonnage stochastique (MCMC) :
   Pour explorer la densité postérieure, les auteurs emploient un algorithme de Monte Carlo par chaînes de Markov (MCMC) basé sur une boucle de Gibbs avec une étape de Metropolis-Hastings (MH) :

   • Échantillonnage de la carte ($x$) : L'échantillonnage de la postérieure conditionnelle de la carte est difficile en raison de sa grande taille et de sa structure de corrélation. Les auteurs utilisent un critère de moindres carrés régularisés perturbé. En ajoutant des perturbations de bruit spécifiques aux données et à la moyenne a priori, le minimiseur de ce critère produit un échantillon distribué exactement selon la densité gaussienne postérieure conditionnelle souhaitée.
   • Échantillonnage des hyperparamètres ($\xi$) : Les postérieures conditionnelles pour les paramètres de précision ($\gamma_n, \gamma_x$) sont des distributions Gamma, et la moyenne du bruit ($m_n$) suit une distribution gaussienne. Ceux-ci sont échantillonnés à l'aide de routines standards.
   • Échantillonnage des paramètres de l'instrument ($\eta$) : La postérieure conditionnelle pour $\eta$ est non standard et complexe. Les auteurs utilisent une étape de marche aléatoire de Metropolis-Hastings . Une proposition est générée, et le rapport d'acceptation est calculé sur la base du rapport des densités postérieures conditionnelles (vraisemblance $\times$ a priori).

Contributions clés

• Cadre d'estimation jointe : L'article présente un cadre bayésien unifié qui estime conjointement la carte, les hyperparamètres et les paramètres de l'instrument, dépassant la séparation traditionnelle entre calibration et inversion.
• Innovation algorithmique :
  • Développement d'une technique d'échantillonnage spécifique pour les cartes gaussiennes de haute dimension utilisant des critères perturbés pour éviter l'inversion directe de matrice.
  • Intégration d'une étape de Metropolis-Hastings pour gérer la postérieure conditionnelle non standard et non linéaire des paramètres de l'instrument au sein d'une boucle de Gibbs.
• Capacité non supervisée et myope : La méthode opère sans « vérité du ciel » (non supervisée) et sans observations de calibration dédiées (myope), en estimant les paramètres de l'instrument directement à partir des données scientifiques.

Résultats
La méthode a été évaluée en utilisant le modèle de l'instrument SPIRE avec des données simulées (Cirrus galactiques, émission étendue et sources ponctuelles) et des observations réelles (NGC 7023).

• Reconstruction de la carte : La méthode non supervisée a produit des cartes visuellement et quantitativement comparables aux « meilleures » cartes supervisées (où les hyperparamètres sont ajustés à la main en utilisant la vérité terrain). L'erreur de reconstruction n'est que légèrement supérieure (3,47 % contre 3,12 % pour les Cirrus galactiques simulés). La reconstruction a réussi à restaurer à la fois les structures à grande échelle et les détails à petite échelle invisibles dans les cartes co-ajoutées naïves.
• Estimation des hyperparamètres : L'algorithme a correctement estimé les paramètres de bruit et de régularité. Pour les données simulées avec une vérité terrain connue, la précision du bruit estimée était à moins de 1 % de la valeur réelle.
• Estimation des paramètres de l'instrument : Dans le scénario myope, le paramètre de l'instrument (lié à la largeur de la fonction d'étalement du point ou PSF) a été estimé avec une erreur d'environ 5 à 7 %, malgré une incertitude a priori de 33 %. L'estimation était plus précise lorsque la source contenait un contenu fréquentiel élevé (sources ponctuelles).
• Quantification de l'incertitude : La méthode fournit naturellement l'écart-type de la postérieure (PSD), offrant une carte d'incertitude qui identifie correctement les régions de haute confiance (zones observées) par rapport aux régions extrapolées.

Signification et affirmations
Les auteurs affirment que ce travail démontre la faisabilité d'une inversion non supervisée et myope dans un contexte astrophysique complexe. Ils soutiennent que :

• La méthode proposée peut ajuster automatiquement les hyperparamètres et les paramètres de l'instrument, éliminant ainsi le besoin d'intervention manuelle ou de balayages de calibration dédiés.
• La qualité des cartes résultantes est presque équivalente à celle obtenue avec une connaissance parfaite de l'instrument et un réglage manuel optimal.
• L'approche fournit un cadre cohérent pour gérer l'équilibre entre l'information a priori et les données d'observation tout en quantifiant l'incertitude.

L'article conclut que ces résultats initiaux sont prometteurs pour le traitement systématique et automatique des flux de données massifs provenant des instruments spatiaux actuels et futurs (par exemple, Herschel, Planck, JWST), tout en reconnaissant que le coût computationnel (environ 10 heures pour le cas myope) et la légère sous-estimation des paramètres de l'instrument dans certains régimes spectraux sont des domaines nécessitant des affinements futurs.