Pixellated Posterior Sampling of Point Spread Functions in Astronomical Images

TL;DR

This paper tackles the problem of precise PSF characterization in astronomical images, where the effective PSF varies across time, position, and wavelength. It introduces a pixel-level Bayesian framework that combines a Gaussian likelihood with a score-based diffusion prior trained on HST ePSF templates to sample the posterior distribution $P(M|D)$ over high-resolution PSFs. The key contributions are a flexible prior for PSF morphology via diffusion models, an iterative posterior-sampling scheme yielding 32 samples per cutout, and a pragmatic interpolation method to propagate PSF uncertainty across the detector. The results show residuals consistent with noise and substantial reduction of flux biases, enabling uncertainty-aware measurements in weak lensing, astrometry, and high-contrast imaging.

Abstract

We introduce a novel framework for upsampled Point Spread Function (PSF) modeling using pixel-level Bayesian inference. Accurate PSF characterization is critical for precision measurements in many fields including: weak lensing, astrometry, and photometry. Our method defines the posterior distribution of the pixelized PSF model through the combination of an analytic Gaussian likelihood and a highly expressive generative diffusion model prior, trained on a library of HST ePSF templates. Compared to traditional methods (parametric Moffat, ePSF template-based, and regularized likelihood), we demonstrate that our PSF models achieve orders of magnitude higher likelihood and residuals consistent with noise, all while remaining visually realistic. Further, the method applies even for faint and heavily masked point sources, merely producing a broader posterior. By recovering a realistic, pixel-level posterior distribution, our technique enables the first meaningful propagation of detailed PSF morphological uncertainty in downstream analysis. An implementation of our posterior sampling procedure is available on GitHub.

Figures (10)

• Figure 1: Methodological comparison. For a bright point source data cutout, we present three traditional modeling results alongside a single posterior sample from our procedure. The Moffat model includes only a smooth decline in brightness and leaves many structured residuals. The ePSF and PSFex models include structural details, resulting in unstructured, though still implausibly large residuals. Our Bayesian posterior sample gives detailed structure and pure noise-like residuals. Yellow pixels are masked for various reasons described in the text. A p-value test on the residuals (described in the text) gives 0.08 for our posterior sample, while for the other algorithms it confidently rejects the null hypothesis ($6\times 10^{-17}$, $2\times10^{-22}$, and $4\times10^{-13}$ respectively).
• Figure 2: Example data from a single HST flt frame. Top: The full frame, with selected stars (red) and the science target (blue). Bottom: Cutouts of the data, variance, and mask for each star. These cutouts illustrate the core challenge: each star presents a unique, complex PSF morphology that must be modeled from noisy, incomplete (masked) data.
• Figure 3: Traditional Method 1: Parametric Moffat Fit. From left to right the columns show the data, the best-fit Moffat model, the $4\times$ high-resolution model, the residuals, and a histogram of the residuals. Yellow pixels in the data and residuals are masked pixels (see text). The p-values in the top right of each histogram are from a Kolmogorov-Smirnov test on the central residuals against a Normal distribution Virtanen2020. The key takeaway is that this method fails: the strong, structured patterns in the residuals (column 4) and the non-Gaussian residual histograms (column 5, with p-values $\approx 0$) prove that a simple analytic function is a misspecified model for the true PSF.
• Figure 4: Traditional Method 2: Fixed ePSF Template. The format is identical to \ref{['fig:astrophot']}, yellow pixels in column 3 indicate pixel values below $10^{-5}$. This represents an improvement: the ePSF model (column 2) captures complex features like diffraction spikes, resulting in smaller residuals than the Moffat model. However, significant $>5\sigma$ residuals remain, proving that a single, fixed "point estimate" template is still an imperfect match to the unique PSF of any given star.
• Figure 5: Traditional Method 3: Empirical PSFex Model. Format is identical to \ref{['fig:astrophot']}. This method shows mixed failure: most stars have implausibly large residuals, while others (e.g., Star 4) have residuals that are too small, indicating the flexible model has overfit the data. This demonstrates that PSFex's simple smoothness regularization is a misspecified prior, motivating our new Bayesian approach.