PSF Regularization in Astronomical Imaging

• PSF regularization is a set of methods that stabilize the estimation and inversion of point spread functions to control noise amplification in astronomical data.
• Techniques range from Fourier-based reciprocal operators and starlet-sparse deconvolution to Laplacian-based kernel smoothing and learned diffusion-prior models.
• These approaches enable robust image homogenization, partial deconvolution for enhanced resolution, and a balanced trade-off between bias and noise.

Searching arXiv for recent and foundational papers on PSF regularization, deconvolution, and PSF homogenization. arXiv search query: "PSF regularization astronomical images deconvolution point spread function regularization" PSF regularization denotes a family of methods that stabilize the estimation, homogenization, matching, or inversion of a point spread function (PSF) when direct deconvolution or kernel fitting is ill-posed. In the astronomical literature, the term spans at least four closely related settings: post-facto homogenization of a spatially varying PSF field into a single target PSF by local transfer kernels (Hughes et al., 2022); sparse-regularized deconvolution and PSF reconstruction using starlet priors and partial deconvolution (Michalewicz et al., 2023); regularized PSF-matching kernels for image subtraction and coaddition (Becker et al., 2012); and blind or myopic deconvolution frameworks in which PSF parameters and image regularization are estimated jointly, often by SURE, Bayesian inference, or null-space methods (Sanders, 2022, Orieux et al., 2010, Bunyak et al., 2015). Across these formulations, the common objective is to control noise amplification and overfitting while preserving the physically relevant structure of the optical response.

1. Scope and problem setting

In imaging instruments, point sources are blurred by a PSF that may be spatially invariant or spatially varying. A representative forward model writes the observed image as the convolution of a true scene with a PSF plus additive noise, either in spatially varying form,

$$I_{\rm obs}(x) = \int I_{\rm true}(u)\,{\rm PSF}_x(x-u)\,du + \epsilon(x),$$

or, in the shift-invariant setting,

$y = Hx + n.$

These models appear explicitly in the formulations of PSF homogenization and deconvolution (Hughes et al., 2022, Michalewicz et al., 2023).

The inverse problem is ill-posed because the PSF suppresses high spatial frequencies and noise is amplified by naïve inversion. In the deconvolution setting, STARRED states that “the deconvolution operation is an ill-posed inverse problem due to noise and pixelization of the data,” and that regularization is necessary to guarantee robustness (Michalewicz et al., 2023). In PSF-matching, the same instability appears when a flexible kernel basis is used: delta-function kernel bases are highly expressive but tend to overfit, yielding noisy kernels with large variance (Becker et al., 2012).

A useful distinction is between regularizing the image reconstruction and regularizing the PSF or kernel itself. Some methods penalize the reconstructed object through Tikhonov, starlet, or geometry-aware functionals (Sanders, 2022, Michalewicz et al., 2023, Bunyak et al., 2015). Others regularize the PSF representation, the inverse PSF, or the matching kernel directly through smoothness penalties, basis restrictions, or learned priors (Becker et al., 2012, Bunyak et al., 2015, Stone et al., 24 Nov 2025). PSF regularization in the narrow sense often refers to the latter, but the literature repeatedly couples both.

2. Post-facto PSF homogenization across a field of view

A particularly direct use of the term appears in “Coma Off It: Removing Variable Point Spread Functions from Astronomical Images” (Hughes et al., 2022). There, PSF regularization is a post-processing method that converts a slowly varying PSF field $\{K_x\}$ into a homogeneous effective PSF $\phi$ across the entire field of view. The defining condition is

$K_x \otimes T_x = \phi,$

where $T_x$ is a local transfer PSF. In Fourier space this becomes

$$\mathcal{K}_x(\omega)\,\mathcal{T}_x(\omega)=\Phi(\omega).$$

Because direct division by $\mathcal{K}_x(\omega)$ is unstable when the modulation transfer function approaches zero, the method introduces a regularized reciprocal operator. Using the paper’s notation,

$$\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr) = [\mathcal{K}_x(\omega)]^* \cdot |\mathcal{K}_x(\omega)|^{\alpha-1} \big/ \left(|\mathcal{K}_x(\omega)|^{\alpha+1} + [\epsilon\,|\Phi(\omega)|]^{\alpha+1}\right),$$

and the transfer-PSF MTF is

$$\mathcal{T}_x(\omega)=\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr)\cdot \Phi(\omega).$$

The spatial-domain kernel is then obtained by inverse FFT (Hughes et al., 2022).

The implementation is tile-based. The full image is partitioned into overlapping $y = Hx + n.$0 neighborhoods with overlap $y = Hx + n.$1 in each direction. For each tile, local PSFs are estimated from bright, isolated stars, patches are normalized to unit total flux and combined by a pixel-wise median, and an optional parametric fit may be used to denoise the PSF estimate (Hughes et al., 2022). A target PSF $y = Hx + n.$2 is then chosen analytically, for example a circular Gaussian with FWHM at least as large as the largest PSF core among the local estimates. Each neighborhood is apodized by a root-Hann window,

$y = Hx + n.$3

to suppress FFT edge artifacts; the processed patches are apodized again and merged with half-tile shifts. The identity

$y = Hx + n.$4

ensures exact amplitude reconstruction where no regularization is applied (Hughes et al., 2022).

The regularization parameters have a concrete operational meaning. The parameter $y = Hx + n.$5 controls maximum amplification, approximately $y = Hx + n.$6, with typical values in $y = Hx + n.$7, while $y = Hx + n.$8 controls the sharpness of the transition between $y = Hx + n.$9 and $\{K_x\}$0, with typical values in $\{K_x\}$1 (Hughes et al., 2022). Computationally, for an $\{K_x\}$2 image and $\{K_x\}$3 patches with stride $\{K_x\}$4, the total cost is approximately

$\{K_x\}$5

with practical runtime described as a few $\{K_x\}$6 when $\{K_x\}$7–$\{K_x\}$8 (Hughes et al., 2022).

The reported examples emphasize homogenization rather than full inversion. A model starfield with slow PSF variation recovered a uniform $\{K_x\}$9 px-FWHM Gaussian target PSF with negligible artifacts; DASH wide-field lens data increased effective resolution from a $\phi$0 px FWHM to a $\phi$1 px FWHM Gaussian with uniform wings; and PUNCH engineering-model images reduced PSF-wing amplitude by approximately $\phi$2–$\phi$3 across $\phi$4 frames (Hughes et al., 2022). This suggests that PSF regularization can be used not only to stabilize inverse problems but also to create spatially uniform imaging properties for mosaics and downstream analysis.

3. Sparse and partial deconvolution frameworks

A second major meaning of PSF regularization arises in deconvolution methods that explicitly regularize the latent image and, in some cases, the PSF model. STARRED formulates the deconvolution problem as

$\phi$5

where $\phi$6 is the starlet analysis operator and the $\phi$7 term enforces sparsity of isotropic undecimated wavelet coefficients (Michalewicz et al., 2023). The approach is “two-channel”: the reconstruction is split as

$\phi$8

with $\phi$9 representing point sources and $K_x \otimes T_x = \phi,$0 representing extended sources on a pixel grid. The combined objective includes a starlet penalty on the extended channel and an optional penalty on the point-source amplitudes (Michalewicz et al., 2023).

This decomposition is motivated by the observation that point-source peaks carry very high spatial frequencies that are poorly represented by starlets alone. Explicit modeling of point sources prevents ringing artifacts in the extended component (Michalewicz et al., 2023). Optimization proceeds by a forward-backward proximal scheme: gradient descent on the data term, followed by soft-thresholding in the starlet domain for the extended component, and an optional proximal or projection step such as non-negativity for the point component (Michalewicz et al., 2023).

STARRED also regularizes PSF estimation itself. Given bright, unsaturated star stamps $K_x \otimes T_x = \phi,$1, the narrow PSF is modeled as

$K_x \otimes T_x = \phi,$2

where $K_x \otimes T_x = \phi,$3 is a circular Gaussian or Moffat core and $K_x \otimes T_x = \phi,$4 is a free pixel-grid residual regularized in the starlet domain. The optimization problem is

$K_x \otimes T_x = \phi,$5

Thus, PSF regularization here is not a separate post-processing stage but an internal constraint on the flexible residual component of the PSF model (Michalewicz et al., 2023).

An important feature is partial deconvolution. Rather than fully removing the PSF, STARRED chooses a target narrow PSF $K_x \otimes T_x = \phi,$6, specified as a Gaussian of FWHM $K_x \otimes T_x = \phi,$7 pixels, and solves for an image convolved by $K_x \otimes T_x = \phi,$8. Equivalently, the forward kernel is replaced by a partial kernel $K_x \otimes T_x = \phi,$9 that has no zero-crossings in the Fourier domain up to the desired cutoff (Michalewicz et al., 2023). The stated rationale is to remain within the Nyquist limit of the up-sampled grid and avoid artifacts.

Parameter selection is described pragmatically. Regularization weights are set by the local noise standard deviation $T_x$0, with a rule-of-thumb $T_x$1; the number of starlet scales is often $T_x$2–$T_x$3; convergence may be declared when the relative change in objective is below $T_x$4 or the iteration count reaches approximately $T_x$5; and a typical $T_x$6 stamp with $T_x$7 scales runs in less than or equal to $T_x$8 s on modern hardware (Michalewicz et al., 2023).

4. Kernel regularization for PSF matching and image subtraction

In image subtraction and coaddition, PSF regularization often concerns the matching kernel rather than the scene or the optical PSF directly. The basic model treats the science image as a reference image convolved with a kernel plus noise,

$T_x$9

or, in the regularized delta-basis formulation,

$$\mathcal{K}_x(\omega)\,\mathcal{T}_x(\omega)=\Phi(\omega).$$0

where the vector $$\mathcal{K}_x(\omega)\,\mathcal{T}_x(\omega)=\Phi(\omega).$$1 contains the kernel coefficients (Becker et al., 2012).

The central problem is that a delta-function basis is maximally flexible: a $$\mathcal{K}_x(\omega)\,\mathcal{T}_x(\omega)=\Phi(\omega).$$2 kernel has $$\mathcal{K}_x(\omega)\,\mathcal{T}_x(\omega)=\Phi(\omega).$$3 free coefficients. This flexibility allows the kernel to represent arbitrary shapes, including off-center PSF differences, astrometric shifts, and optical distortions, but it also “fits the noise,” producing kernels with high spatial variance and residual variance in the difference image below unity, identified as a sign of over-fitting (Becker et al., 2012).

Regularization is introduced through a penalized least-squares criterion,

$$\mathcal{K}_x(\omega)\,\mathcal{T}_x(\omega)=\Phi(\omega).$$4

where $$\mathcal{K}_x(\omega)\,\mathcal{T}_x(\omega)=\Phi(\omega).$$5 is constructed from a finite-difference approximation to the 2D Laplacian. Concretely, each row of $$\mathcal{K}_x(\omega)\,\mathcal{T}_x(\omega)=\Phi(\omega).$$6 corresponds to the five-point stencil

$$\mathcal{K}_x(\omega)\,\mathcal{T}_x(\omega)=\Phi(\omega).$$7

so the penalty $$\mathcal{K}_x(\omega)\,\mathcal{T}_x(\omega)=\Phi(\omega).$$8 suppresses rough kernels (Becker et al., 2012). The normal equations yield the closed-form estimator

$$\mathcal{K}_x(\omega)\,\mathcal{T}_x(\omega)=\Phi(\omega).$$9

The regularization strength mediates a bias-variance trade-off. As $\mathcal{K}_x(\omega)$0, one recovers the unregularized, overfitting solution; as $\mathcal{K}_x(\omega)$1, the kernel is forced toward the smoothest shape but underfits, increasing residual variance in the difference image (Becker et al., 2012). The paper discusses both Stein’s unbiased risk estimator and generalized cross-validation:

$\mathcal{K}_x(\omega)$2

with

$\mathcal{K}_x(\omega)$3

and

$\mathcal{K}_x(\omega)$4

Both criteria are reported to suggest moderate regularization, and practical experiments show a broad “sweet-spot” of $\mathcal{K}_x(\omega)$5–$\mathcal{K}_x(\omega)$6 (Becker et al., 2012).

This formulation is distinct from image deconvolution, because it regularizes a matching kernel whose role is to bring images to a common PSF for subtraction or coaddition. Nevertheless, the underlying principle is the same: a highly expressive kernel basis requires explicit control of roughness to achieve stable and spatially interpolable solutions (Becker et al., 2012).

5. Joint estimation of PSF and regularization parameters

Several papers treat PSF regularization as part of a coupled parameter-estimation problem in which blur parameters and regularization strengths are optimized simultaneously. In the SURE-based blind-deconvolution framework of Sanders, the deconvolved image for a candidate PSF $\mathcal{K}_x(\omega)$7 is

$\mathcal{K}_x(\omega)$8

where $\mathcal{K}_x(\omega)$9 is a finite-difference or similar regularizer (Sanders, 2022). The quantity to be minimized is Stein’s unbiased risk estimator,

$$\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr) = [\mathcal{K}_x(\omega)]^* \cdot |\mathcal{K}_x(\omega)|^{\alpha-1} \big/ \left(|\mathcal{K}_x(\omega)|^{\alpha+1} + [\epsilon\,|\Phi(\omega)|]^{\alpha+1}\right),$$0

with $$\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr) = [\mathcal{K}_x(\omega)]^* \cdot |\mathcal{K}_x(\omega)|^{\alpha-1} \big/ \left(|\mathcal{K}_x(\omega)|^{\alpha+1} + [\epsilon\,|\Phi(\omega)|]^{\alpha+1}\right),$$1 (Sanders, 2022).

The paper derives fixed-point updates for both PSF-shape parameters and the regularization parameter. For $$\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr) = [\mathcal{K}_x(\omega)]^* \cdot |\mathcal{K}_x(\omega)|^{\alpha-1} \big/ \left(|\mathcal{K}_x(\omega)|^{\alpha+1} + [\epsilon\,|\Phi(\omega)|]^{\alpha+1}\right),$$2, the exact root condition gives

$$\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr) = [\mathcal{K}_x(\omega)]^* \cdot |\mathcal{K}_x(\omega)|^{\alpha-1} \big/ \left(|\mathcal{K}_x(\omega)|^{\alpha+1} + [\epsilon\,|\Phi(\omega)|]^{\alpha+1}\right),$$3

which induces a fixed-point iteration (Sanders, 2022). For Gaussian PSF parameters, empirical tuning recommends an exponent $$\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr) = [\mathcal{K}_x(\omega)]^* \cdot |\mathcal{K}_x(\omega)|^{\alpha-1} \big/ \left(|\mathcal{K}_x(\omega)|^{\alpha+1} + [\epsilon\,|\Phi(\omega)|]^{\alpha+1}\right),$$4 in the fixed-point map; nearly $$\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr) = [\mathcal{K}_x(\omega)]^* \cdot |\mathcal{K}_x(\omega)|^{\alpha-1} \big/ \left(|\mathcal{K}_x(\omega)|^{\alpha+1} + [\epsilon\,|\Phi(\omega)|]^{\alpha+1}\right),$$5 of Gaussian-PSF trials are reported to converge in fewer than $$\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr) = [\mathcal{K}_x(\omega)]^* \cdot |\mathcal{K}_x(\omega)|^{\alpha-1} \big/ \left(|\mathcal{K}_x(\omega)|^{\alpha+1} + [\epsilon\,|\Phi(\omega)|]^{\alpha+1}\right),$$6 iterations, typically $$\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr) = [\mathcal{K}_x(\omega)]^* \cdot |\mathcal{K}_x(\omega)|^{\alpha-1} \big/ \left(|\mathcal{K}_x(\omega)|^{\alpha+1} + [\epsilon\,|\Phi(\omega)|]^{\alpha+1}\right),$$7–$$\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr) = [\mathcal{K}_x(\omega)]^* \cdot |\mathcal{K}_x(\omega)|^{\alpha-1} \big/ \left(|\mathcal{K}_x(\omega)|^{\alpha+1} + [\epsilon\,|\Phi(\omega)|]^{\alpha+1}\right),$$8 (Sanders, 2022). Here PSF regularization is not an isolated penalty on the PSF itself; rather, it is the automatic tuning of deconvolution regularization in response to PSF estimates.

A fully Bayesian analogue appears in the Wiener–Hunt framework of Orieux, Giovannelli, and Rodet (Orieux et al., 2010). The observation model is

$$\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr) = [\mathcal{K}_x(\omega)]^* \cdot |\mathcal{K}_x(\omega)|^{\alpha-1} \big/ \left(|\mathcal{K}_x(\omega)|^{\alpha+1} + [\epsilon\,|\Phi(\omega)|]^{\alpha+1}\right),$$9

with a PSF parameter vector $$\mathcal{T}_x(\omega)=\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr)\cdot \Phi(\omega).$$0 and a Gaussian Markov random field prior on the image,

$$\mathcal{T}_x(\omega)=\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr)\cdot \Phi(\omega).$$1

Gamma hyperpriors are placed on the precisions, a uniform prior is placed on $$\mathcal{T}_x(\omega)=\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr)\cdot \Phi(\omega).$$2, and inference proceeds from the joint posterior over image, precisions, and PSF parameters by Gibbs sampling and Metropolis–Hastings updates (Orieux et al., 2010).

The paper emphasizes interaction between the image smoothness precision $$\mathcal{T}_x(\omega)=\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr)\cdot \Phi(\omega).$$3 and the PSF widths $$\mathcal{T}_x(\omega)=\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr)\cdot \Phi(\omega).$$4. If the estimated PSF is looser, the chain tends to push $$\mathcal{T}_x(\omega)=\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr)\cdot \Phi(\omega).$$5 downward; if the PSF is estimated narrower than truth, the sampler increases $$\mathcal{T}_x(\omega)=\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr)\cdot \Phi(\omega).$$6 to enforce stronger smoothing and avoid amplifying noise (Orieux et al., 2010). Joint histograms show a clear negative correlation between $$\mathcal{T}_x(\omega)=\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr)\cdot \Phi(\omega).$$7 and PSF widths. In the reported $$\mathcal{T}_x(\omega)=\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr)\cdot \Phi(\omega).$$8 simulation, posterior-mean estimates recovered $$\mathcal{T}_x(\omega)=\mathcal{R}_{\alpha,\epsilon,P}\bigl(\mathcal{K}_x(\omega)\bigr)\cdot \Phi(\omega).$$9, $y = Hx + n.$00, $y = Hx + n.$01, $y = Hx + n.$02, and $y = Hx + n.$03 for true values $y = Hx + n.$04, $y = Hx + n.$05, $y = Hx + n.$06, $y = Hx + n.$07, and $y = Hx + n.$08 (Orieux et al., 2010).

These approaches illustrate a recurrent theme: the amount of regularization cannot be meaningfully separated from the inferred PSF. A plausible implication is that “PSF regularization” in myopic or blind deconvolution is best understood as a coupled estimation problem rather than a single penalty term.

6. Surface-area, null-space, and geometry-aware regularization

A different line of work constructs PSF regularization from null-space structure and geometric smoothness. In Bunyak et al., the blurred image is modeled as

$y = Hx + n.$09

and a shift-invariant autoregressive operator $y = Hx + n.$10 is assumed to annihilate both the original and blurred images. In discrete form,

$y = Hx + n.$11

From this, the lexicographically arranged PSF vector $y = Hx + n.$12 must satisfy

$y = Hx + n.$13

so the PSF lies in the left, or conjugated, null-space of the AR operator matrix (Bunyak et al., 2015). A symmetric operator $y = Hx + n.$14 is then diagonalized, and the PSF and inverse PSF are expanded in the eigenfunctions corresponding to near-zero eigenvalues (Bunyak et al., 2015).

Regularization is imposed through the surface-area functional

$y = Hx + n.$15

or, for the inverse PSF,

$y = Hx + n.$16

This penalty is combined with data fidelity in variational objectives for PSF and IPSF optimization (Bunyak et al., 2015, Bunyak et al., 2012). In the IPSF case, the Euler–Lagrange condition is

$y = Hx + n.$17

subject to nonnegativity, compact support, and normalization constraints on the PSF or IPSF as stated in the 2012 formulation (Bunyak et al., 2012).

The same framework extends to deconvolution of the image estimate. One update has the form

$y = Hx + n.$18

with a dynamic regularization coefficient $y = Hx + n.$19 chosen to balance smoothing and sharpening variations (Bunyak et al., 2012). The practical convergence criteria include

$y = Hx + n.$20

and

$y = Hx + n.$21

The dynamic update is intended to prevent overshoot, with $y = Hx + n.$22 as $y = Hx + n.$23 (Bunyak et al., 2012).

A further extension treats the image estimate as a 2D manifold in $y = Hx + n.$24 with metric

$y = Hx + n.$25

or equivalently

$y = Hx + n.$26

The fidelity functional is weighted by this area element, yielding a curved-space regularization scheme in which the regularization weight becomes spatially varying through a factor $y = Hx + n.$27 (Bunyak et al., 2015, Bunyak et al., 2012). This differs from Tikhonov or TV by coupling fidelity and geometry directly.

7. Learned priors, posterior sampling, and uncertainty propagation

Recent work extends PSF regularization beyond explicit smoothness penalties and basis constraints to learned generative priors. In “Pixellated Posterior Sampling of Point Spread Functions in Astronomical Images,” the unknown supersampled PSF image $y = Hx + n.$28 is inferred from noisy star cutouts through a Gaussian likelihood combined with a diffusion-model prior trained on a library of HST ePSF templates (Stone et al., 24 Nov 2025). The likelihood is

$y = Hx + n.$29

where $y = Hx + n.$30 renders the supersampled PSF at fractional pixel centers and $y = Hx + n.$31 is the diagonal covariance from the ERR map (Stone et al., 24 Nov 2025).

The regularizing prior is defined implicitly by a score-based diffusion model. A forward noising process gradually maps realistic PSFs to isotropic Gaussian noise, and a U-Net score estimator learns $y = Hx + n.$32 at intermediate noise levels (Stone et al., 24 Nov 2025). The posterior combines the likelihood, the diffusion-model prior on $y = Hx + n.$33, and broad priors on shifts, flux, and sky background. Sampling is then performed through a posterior SDE that balances the learned prior score with the likelihood gradient,

$y = Hx + n.$34

integrated backward from $y = Hx + n.$35 to $y = Hx + n.$36 using a Heun scheme with Langevin corrector steps (Stone et al., 24 Nov 2025).

In this framework, regularization is explicitly described as data-driven. The diffusion prior replaces “classical smoothness or parametric priors,” enforces that sampled PSFs lie on the learned manifold of realistic HST ePSFs, and adapts continuously to the signal-to-noise ratio: for faint or heavily masked sources, the posterior broadens, with control shifting from the data-informed core to the prior in the outskirts (Stone et al., 24 Nov 2025). Reported evaluation on HST ACS/WFC F814W cutouts finds posterior samples with orders-of-magnitude higher likelihood than three classical baselines, residual images statistically indistinguishable from white noise with $y = Hx + n.$37 and KS-test $y = Hx + n.$38-values around $y = Hx + n.$39, and photometric bias of classical models around $y = Hx + n.$40 versus an unbiased posterior mean flux with a faithful credible interval (Stone et al., 24 Nov 2025).

This development changes the interpretation of PSF regularization. Rather than selecting a single penalty coefficient that trades smoothness against variance, one regularizes by restricting inference to a high-dimensional learned distribution. This suggests a shift from deterministic stabilization toward posterior inference with explicit morphological uncertainty.

8. Comparative themes and recurring misconceptions

Several themes recur across these otherwise different methods. First, regularization does not necessarily imply full deblurring. STARRED deliberately avoids complete inversion and instead “brings the image to a higher resolution” through partial deconvolution (Michalewicz et al., 2023). Likewise, field-wide PSF regularization in Hughes et al. targets a homogeneous effective PSF rather than the diffraction limit (Hughes et al., 2022).

Second, regularization is not equivalent to arbitrary smoothing. In PSF-matching, the parameter $y = Hx + n.$41 exchanges variance in the difference image with variance in the kernel itself (Becker et al., 2012). In SURE- and Bayesian methods, the optimal regularization level is coupled to the inferred PSF parameters (Sanders, 2022, Orieux et al., 2010). In null-space and surface-area formulations, the regularizer is explicitly geometrical rather than purely quadratic (Bunyak et al., 2015, Bunyak et al., 2012). In diffusion-prior methods, the constraint is learned from a training distribution and is intended to prevent both over-smoothing and over-fitting (Stone et al., 24 Nov 2025).

Third, PSF regularization is often conflated with PSF estimation alone. The literature separates several operations: estimation of a local PSF field from stars; stabilization of that estimate by parametric, sparse, or learned priors; computation of transfer kernels or inverse filters; and regularized reconstruction of the latent scene. Different papers place the regularization at different stages (Hughes et al., 2022, Michalewicz et al., 2023, Becker et al., 2012, Stone et al., 24 Nov 2025).

A concise comparison is useful.

Setting	Regularized object	Characteristic mechanism
Field homogenization (Hughes et al., 2022)	Transfer PSF $y = Hx + n.$42	Regularized reciprocal in Fourier space
Two-channel deconvolution (Michalewicz et al., 2023)	Extended image and PSF residual grid	Starlet $y = Hx + n.$43 sparsity and partial deconvolution
Image subtraction (Becker et al., 2012)	PSF-matching kernel $y = Hx + n.$44	Laplacian quadratic penalty
Blind/myopic deconvolution (Sanders, 2022, Orieux et al., 2010)	Image prior strength and PSF parameters jointly	SURE fixed points or Bayesian posterior inference
Null-space/IPSF methods (Bunyak et al., 2015, Bunyak et al., 2012)	PSF, IPSF, and image estimate	Surface-area and curved-space regularization
Posterior PSF sampling (Stone et al., 24 Nov 2025)	Pixelized PSF morphology	Diffusion-model prior and posterior SDE

Taken together, these works define PSF regularization as a broad methodological class for stabilizing PSF-related inverse problems. The specific mathematical form may be a regularized reciprocal, an $y = Hx + n.$45 sparse prior, a Laplacian quadratic penalty, a surface-area functional, a hierarchical Bayesian prior, or a score-based diffusion model. What unifies them is the attempt to preserve scientifically relevant optical structure while preventing the instability that follows from direct inversion, excessive basis flexibility, or underconstrained pixelized PSF models.