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

Technical Summary: Estimating Hyperparameters and Instrument Parameters in Regularized Inversion

Problem Statement
The paper addresses the challenge of image reconstruction in astrophysics, specifically for the SPIRE instrument on the Herschel space observatory. The core problem is an "unsupervised and myopic" inverse problem. Standard inversion methods typically require:

1. Supervised tuning: Hyperparameters (which balance data fidelity against prior regularization) must be set manually, often requiring knowledge of the "true" sky or extensive trial-and-error.
2. Known instrument models: Instrument parameters (such as lobe width, response time, and gain) are usually determined via dedicated calibration observations. These values inevitably contain errors due to spectral dependencies and calibration uncertainties.

The authors aim to solve a threefold problem from a single observation: simultaneously estimating the sky map ($x$), the hyperparameters ($\xi$), and the instrument parameters ($\eta$) without prior knowledge of the ground truth or dedicated calibration data.

Methodology
The proposed approach is grounded in a Full-Bayesian framework. The authors construct a posterior density $p(x, \xi, \eta | y)$ for all unknown quantities given the observations $y$.

1. Probabilistic Modeling:

   • Likelihood: The observation model is linear ($y = A_\eta x + n$), where noise $n$ is modeled as Gaussian with a mean $m_n$ and covariance structure defined by a precision parameter $\gamma_n$.
   • Prior for the Map: The sky map $x$ is modeled using a Gaussian prior to enforce spatial regularity (smoothness), controlled by a precision parameter $\gamma_x$. The regularity structure is defined via a quadratic norm of spatial derivatives.
   • Hyperpriors: The hyperparameters ($\xi = [m_n, \gamma_n, \gamma_x]$) are assigned non-informative Jeffreys' priors ($1/\gamma$ for precisions, uniform for means) to minimize prior bias.
   • Instrument Parameters: The instrument parameters $\eta$ are modeled with Gaussian priors centered on nominal values with associated uncertainties.

2. Posterior Density:
   The joint posterior density is derived by combining the likelihood and the priors. Due to the high dimensionality of the map ($M \sim 10^5$) and the complex, non-linear dependence of the observations on the instrument parameters $\eta$, the posterior cannot be analyzed directly or maximized analytically.

3. Stochastic Sampling (MCMC):
   To explore the posterior density, the authors employ a Markov Chain Monte-Carlo (MCMC) algorithm based on a Gibbs loop with a Metropolis-Hastings (MH) step:

   • Map Sampling ($x$): Sampling from the conditional posterior of the map is difficult due to the large size and correlation structure. The authors utilize a perturbed regularized least-squares criterion. By adding specific noise perturbations to the data and the prior mean, the minimizer of this criterion yields a sample distributed exactly according to the desired conditional Gaussian posterior.
   • Hyperparameter Sampling ($\xi$): The conditional posteriors for the precision parameters ($\gamma_n, \gamma_x$) are Gamma distributions, and the mean noise ($m_n$) follows a Gaussian distribution. These are sampled using standard routines.
   • Instrument Parameter Sampling ($\eta$): The conditional posterior for $\eta$ is non-standard and complex. The authors use a Metropolis-Hastings random-walk step. A proposal is generated, and the acceptance ratio is computed based on the ratio of the conditional posterior densities (likelihood $\times$ prior).

Key Contributions

• Joint Estimation Framework: The paper presents a unified Bayesian framework that jointly estimates the map, hyperparameters, and instrument parameters, moving beyond the traditional separation of calibration and inversion.
• Algorithmic Innovation:
  • Development of a specific sampling technique for high-dimensional Gaussian maps using perturbed criteria to avoid direct matrix inversion.
  • Integration of a Metropolis-Hastings step to handle the non-standard, non-linear conditional posterior of instrument parameters within a Gibbs loop.
• Unsupervised and Myopic Capability: The method operates without "sky truth" (unsupervised) and without dedicated calibration observations (myopic), estimating instrument parameters directly from the science data.

Results
The method was evaluated using the SPIRE instrument model with both simulated data (Galactic Cirrus, extended emission, and point sources) and real observations (NGC 7023).

• Map Reconstruction: The unsupervised method produced maps visually and quantitatively comparable to "best" supervised maps (where hyperparameters are tuned by hand using ground truth). The reconstruction error was only slightly higher (3.47% vs. 3.12% for simulated Galactic Cirrus). The method successfully restored both large-scale structures and small-scale details invisible in naive co-added maps.
• Hyperparameter Estimation: The algorithm correctly estimated noise and regularity parameters. For simulated data with known ground truth, the estimated noise precision was within 1% of the true value.
• Instrument Parameter Estimation: In the myopic scenario, the instrument parameter (related to the Point Spread Function width) was estimated with an error of approximately 5–7%, despite a prior uncertainty of 33%. The estimation was more accurate when the source contained high-frequency content (point sources).
• Uncertainty Quantification: The method naturally provides the posterior standard deviation (PSD), offering a map of uncertainty that correctly identifies regions with high confidence (observed areas) versus extrapolated regions.

Significance and Claims
The authors claim that this work demonstrates the feasibility of performing unsupervised and myopic inversion in a complex astrophysical context. They assert that:

• The proposed method can automatically tune hyperparameters and instrument parameters, eliminating the need for manual intervention or dedicated calibration scans.
• The quality of the resulting maps is nearly equivalent to those obtained with perfect knowledge of the instrument and optimal manual tuning.
• The approach provides a consistent framework for handling the trade-off between prior information and observational data while quantifying uncertainty.

The paper concludes that these initial results are promising for the systematic, automatic processing of large data streams from current and future space-based instruments (e.g., Herschel, Planck, JWST), though it acknowledges the computational cost (approx. 10 hours for the myopic case) and the slight underestimation of instrument parameters in certain spectral regimes as areas for future refinement.