---
title: PSF Regularization in Astronomical Imaging
url: https://www.emergentmind.com/topics/psf-regularization
type: topic
---

# PSF Regularization in Astronomical Imaging

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 [2212.02594]; sparse-regularized deconvolution and PSF reconstruction using starlet priors and partial deconvolution [2305.18526]; regularized PSF-matching kernels for image subtraction and coaddition [1202.2902]; 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 [2202.13242], [1004.5538], [1502.07781]. 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 [2212.02594], [2305.18526].

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 [2305.18526]. 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 [1202.2902].

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 [2202.13242], [2305.18526], [1502.07781]. Others regularize the PSF representation, the inverse PSF, or the matching kernel directly through smoothness penalties, basis restrictions, or learned priors [1202.2902], [1502.07781], [2511.19594]. 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” [2212.02594]. 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 [2212.02594].

The implementation is tile-based. The full image is partitioned into overlapping $N\times N$ neighborhoods with overlap $N/2$ 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 [2212.02594]. A target PSF $\phi$ 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,
$$
w(x,y) = \sin\!\left[\frac{(x+1/2)\pi}{N}\right]\sin\!\left[\frac{(y+1/2)\pi}{N}\right],
$$
to suppress FFT edge artifacts; the processed patches are apodized again and merged with half-tile shifts. The identity
$$
w^2(x,y)+w^2(x+N/2,y)+\dots=1
$$
ensures exact amplitude reconstruction where no regularization is applied [2212.02594].

The regularization parameters have a concrete operational meaning. The parameter $\epsilon$ controls maximum amplification, approximately $1/\epsilon$, with typical values in $[0.1,0.3]$, while $\alpha$ controls the sharpness of the transition between $1/\mathcal{K}$ and $0$, with typical values in $[2,10]$ [2212.02594]. Computationally, for an $L\times L$ image and $N\times N$ patches with stride $N/2$, the total cost is approximately
$$
O(L^2\log N),
$$
with practical runtime described as a few $\times L^2\log L$ when $N\approx L/8$–$L/16$ [2212.02594].

The reported examples emphasize homogenization rather than full inversion. A model starfield with slow PSF variation recovered a uniform $3$ px-FWHM Gaussian target PSF with negligible artifacts; DASH wide-field lens data increased effective resolution from a $7.5$ px FWHM to a $3.5$ px FWHM Gaussian with uniform wings; and PUNCH engineering-model images reduced PSF-wing amplitude by approximately $30$–$100\times$ across $540$ frames [2212.02594]. 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
$$
\min_{x\in\mathbb R^N}\;
\frac12\|Hx-y\|_2^2 + \lambda \|Wx\|_1,
$$
where $W$ is the starlet analysis operator and the $\ell_1$ term enforces sparsity of isotropic undecimated wavelet coefficients [2305.18526]. The approach is “two-channel”: the reconstruction is split as
$$
x = x_{\rm pt} + x_{\rm ext},
$$
with $x_{\rm pt}$ representing point sources and $x_{\rm ext}$ 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 [2305.18526].

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 [2305.18526]. 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 [2305.18526].

STARRED also regularizes PSF estimation itself. Given bright, unsaturated star stamps $\{s_i\}_{i=1}^P$, the narrow PSF is modeled as
$$
k = k_{\rm anal} + d,
$$
where $k_{\rm anal}$ is a circular Gaussian or Moffat core and $d$ is a free pixel-grid residual regularized in the starlet domain. The optimization problem is
$$
\min_d \sum_{i=1}^P \frac12\|k*s_i^{\rm ref}-s_i^{\rm obs}\|_2^2 + \mu \|Wd\|_1.
$$
Thus, PSF regularization here is not a separate post-processing stage but an internal constraint on the flexible residual component of the PSF model [2305.18526].

An important feature is **partial deconvolution**. Rather than fully removing the PSF, STARRED chooses a target narrow PSF $r$, specified as a Gaussian of FWHM $=2$ pixels, and solves for an image convolved by $r$. Equivalently, the forward kernel is replaced by a partial kernel $K=r^{-1}*H$ that has no zero-crossings in the Fourier domain up to the desired cutoff [2305.18526]. 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 $\sigma$, with a rule-of-thumb $\lambda\approx \sigma\sqrt{2\ln N}$; the number of starlet scales is often $J=4$–$6$; convergence may be declared when the relative change in objective is below $10^{-5}$ or the iteration count reaches approximately $200$; and a typical $64\times64$ stamp with $6$ scales runs in less than or equal to $1$ s on modern hardware [2305.18526].

## 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,
$$
\mathbf{S} = \mathbf{C}\mathbf{a} + \boldsymbol\epsilon,
$$
or, in the regularized delta-basis formulation,
$$
\mathbf{D}=A\mathbf{k}+\boldsymbol\epsilon,
$$
where the vector $\mathbf{k}$ contains the kernel coefficients [1202.2902].

The central problem is that a delta-function basis is maximally flexible: a $19\times19$ kernel has $361$ 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 [1202.2902].

Regularization is introduced through a penalized least-squares criterion,
$$
J(\mathbf{k}) = \|\mathbf{D}-A\mathbf{k}\|^2 + \lambda\,\mathbf{k}^T L\,\mathbf{k},
$$
where $L=R^TR$ is constructed from a finite-difference approximation to the 2D Laplacian. Concretely, each row of $R$ corresponds to the five-point stencil
$$
k_{i+1,j}+k_{i-1,j}+k_{i,j+1}+k_{i,j-1}-4k_{i,j},
$$
so the penalty $\mathbf{k}^T L\mathbf{k} = \|R\mathbf{k}\|^2$ suppresses rough kernels [1202.2902]. The normal equations yield the closed-form estimator
$$
\hat{\mathbf{k}}
=
(A^TA+\lambda L)^{-1}A^T\mathbf{D}.
$$

The regularization strength mediates a bias-variance trade-off. As $\lambda\to0$, one recovers the unregularized, overfitting solution; as $\lambda\to\infty$, the kernel is forced toward the smoothest shape but underfits, increasing residual variance in the difference image [1202.2902]. The paper discusses both Stein’s unbiased risk estimator and generalized cross-validation:
$$
\hat R(\lambda)=\frac{1}{n}\|\mathbf{D}-A\hat{\mathbf{k}}\|^2
+\frac{2\sigma^2}{n}\operatorname{tr}(S(\lambda))-\sigma^2,
$$
with
$$
S(\lambda)=A(A^TA+\lambda L)^{-1}A^T,
$$
and
$$
\operatorname{GCV}(\lambda)
=
\frac{\|\mathbf{D}-A\hat{\mathbf{k}}\|^2}{[n-\operatorname{tr}S(\lambda)]^2}.
$$
Both criteria are reported to suggest moderate regularization, and practical experiments show a broad “sweet-spot” of $\lambda\approx0.1$–$1.0$ [1202.2902].

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 [1202.2902].

## 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 $H$ is
$$
u(\lambda)=\arg\min_v \|Hv-b\|_2^2 + \lambda \|Tv\|_2^2
=
(H^TH+\lambda T^TT)^{-1}H^Tb,
$$
where $T$ is a finite-difference or similar regularizer [2202.13242]. The quantity to be minimized is Stein’s unbiased risk estimator,
$$
\operatorname{SURE}(u)
=
-N\sigma^2+\|Hu-b\|_2^2+2\sigma^2\operatorname{trace}(HMH^T),
$$
with $M=(H^TH+\lambda T^TT)^{-1}$ [2202.13242].

The paper derives fixed-point updates for both PSF-shape parameters and the regularization parameter. For $\lambda$, the exact root condition gives
$$
\lambda
=
\frac{\sigma^2\,\operatorname{trace}(HM T^TT MH^T)}
{b^THM T^TT M T^TT MH^Tb},
$$
which induces a fixed-point iteration [2202.13242]. For Gaussian PSF parameters, empirical tuning recommends an exponent $p=0.25$ in the fixed-point map; nearly $100\%$ of Gaussian-PSF trials are reported to converge in fewer than $200$ iterations, typically $50$–$100$ [2202.13242]. 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 [1004.5538]. The observation model is
$$
y = H_w x + \varepsilon,
\qquad
\varepsilon\sim\mathcal N(0,\gamma_\varepsilon^{-1}I),
$$
with a PSF parameter vector $w$ and a Gaussian Markov random field prior on the image,
$$
\Lambda_x=\gamma_0 e_0 e_0^T + \gamma_1 D^TD.
$$
Gamma hyperpriors are placed on the precisions, a uniform prior is placed on $w$, and inference proceeds from the joint posterior over image, precisions, and PSF parameters by Gibbs sampling and Metropolis–Hastings updates [1004.5538].

The paper emphasizes interaction between the image smoothness precision $\gamma_1$ and the PSF widths $w$. If the estimated PSF is looser, the chain tends to push $\gamma_1$ downward; if the PSF is estimated narrower than truth, the sampler increases $\gamma_1$ to enforce stronger smoothing and avoid amplifying noise [1004.5538]. Joint histograms show a clear negative correlation between $\gamma_1$ and PSF widths. In the reported $128\times128$ simulation, posterior-mean estimates recovered $\hat\gamma_\varepsilon=0.49$, $\hat\gamma_1=1.65$, $\hat w_a=20.07$, $\hat w_b=7.19$, and $\hat\varphi=1.03$ for true values $\gamma_\varepsilon=0.5$, $\gamma_1=2$, $w_a=20$, $w_b=7$, and $\varphi=\pi/3$ [1004.5538].

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
$$
X(x,y)=[H*S](x,y),
$$
and a shift-invariant autoregressive operator $\mathcal A$ is assumed to annihilate both the original and blurred images. In discrete form,
$$
\sum_{i=0}^{P-1}\sum_{k=0}^{Q-1} a_{i,k}\,X_{n-i,m-k}=0.
$$
From this, the lexicographically arranged PSF vector $\mathbf h$ must satisfy
$$
\mathbf h^T \mathbf A = \mathbf 0,
$$
so the PSF lies in the left, or conjugated, null-space of the AR operator matrix [1502.07781]. A symmetric operator $\mathbf B=\mathbf A\mathbf A^T$ is then diagonalized, and the PSF and inverse PSF are expanded in the eigenfunctions corresponding to near-zero eigenvalues [1502.07781].

Regularization is imposed through the surface-area functional
$$
R(H)=\iint_\Omega \sqrt{1+H_x^2+H_y^2}\,dx\,dy,
$$
or, for the inverse PSF,
$$
E_G[G]=\iint_\Xi \sqrt{1+G_x^2+G_y^2}\,dx\,dy.
$$
This penalty is combined with data fidelity in variational objectives for PSF and IPSF optimization [1502.07781], [1206.3594]. In the IPSF case, the Euler–Lagrange condition is
$$
Y^\top*(G*Y-x)
-
\lambda \nabla\!\cdot\!\left(\frac{\nabla G}{\sqrt{1+|\nabla G|^2}}\right)
=0,
$$
subject to nonnegativity, compact support, and normalization constraints on the PSF or IPSF as stated in the 2012 formulation [1206.3594].

The same framework extends to deconvolution of the image estimate. One update has the form
$$
S^{(k+1)}
=
S^{(k)}
+
\tau\Bigl[
X-H*S^{(k)}
+
\alpha^{(k)}\,G*L(S^{(k)})
\Bigr],
$$
with a dynamic regularization coefficient $\alpha^{(k)}$ chosen to balance smoothing and sharpening variations [1206.3594]. The practical convergence criteria include
$$
\|S^{(k+1)}-S^{(k)}\|_2 \le \delta \|S^{(k)}-S^{(k-1)}\|_2,
\qquad
\delta\approx1.0001,
$$
and
$$
\|S^{(k+1)}-S^{(k)}\|_2 < \varepsilon,
\qquad
\varepsilon\approx 10^{-8}.
$$
The dynamic update is intended to prevent overshoot, with $\alpha^{(k)}\to0$ as $\|\Delta S\|\to0$ [1206.3594].

A further extension treats the image estimate as a 2D manifold in $\mathbb R^3$ with metric
$$
g_{ij}=\delta_{ij}+S_{,i}S_{,j},
$$
or equivalently
$$
\sqrt{\det \mathbf g}=\sqrt{1+S_x^2+S_y^2}.
$$
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 $A(H,S)$ [1502.07781], [1206.3594]. 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 $M$ is inferred from noisy star cutouts through a Gaussian likelihood combined with a diffusion-model prior trained on a library of HST ePSF templates [2511.19594]. The likelihood is
$$
\mathcal L(\mathbf d\mid M,x,y,f,s)
=
\frac{1}{\sqrt{(2\pi)^N\det\Sigma}\,
\exp\!\left[
-\frac12(\mathbf d-f\,\mathbf F(M,x,y)-s)^\top
\Sigma^{-1}
(\mathbf d-f\,\mathbf F(M,x,y)-s)
\right]},
$$
where $\mathbf F(M,x,y)$ renders the supersampled PSF at fractional pixel centers and $\Sigma$ is the diagonal covariance from the ERR map [2511.19594].

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 $\nabla_{M_t}\log q_t(M_t)$ at intermediate noise levels [2511.19594]. The posterior combines the likelihood, the diffusion-model prior on $M$, 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,
$$
dM_t
=
\Bigl[
f(t)\,\mathbf s_\theta(M_t,t)
+
g(t)^2 \nabla_{M_t}\log \mathcal L(\mathbf d\mid M_t)
\Bigr]dt
+
g(t)\,d\mathbf W_t,
$$
integrated backward from $t=T$ to $0$ using a Heun scheme with Langevin corrector steps [2511.19594].

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 [2511.19594]. 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 $\chi^2/\mathrm{dof}\approx1$ and KS-test $p$-values around $0.5$, and photometric bias of classical models around $5\%$ versus an unbiased posterior mean flux with a faithful credible interval [2511.19594].

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 [2305.18526]. Likewise, field-wide PSF regularization in Hughes et al. targets a homogeneous effective PSF rather than the diffraction limit [2212.02594].

Second, regularization is not equivalent to arbitrary smoothing. In PSF-matching, the parameter $\lambda$ exchanges variance in the difference image with variance in the kernel itself [1202.2902]. In SURE- and Bayesian methods, the optimal regularization level is coupled to the inferred PSF parameters [2202.13242], [1004.5538]. In null-space and surface-area formulations, the regularizer is explicitly geometrical rather than purely quadratic [1502.07781], [1206.3594]. In diffusion-prior methods, the constraint is learned from a training distribution and is intended to prevent both over-smoothing and over-fitting [2511.19594].

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 [2212.02594], [2305.18526], [1202.2902], [2511.19594].

A concise comparison is useful.

| Setting | Regularized object | Characteristic mechanism |
|---|---|---|
| Field homogenization [2212.02594] | Transfer PSF $T_x$ | Regularized reciprocal in Fourier space |
| Two-channel deconvolution [2305.18526] | Extended image and PSF residual grid | Starlet $\ell_1$ sparsity and partial deconvolution |
| Image subtraction [1202.2902] | PSF-matching kernel $\mathbf{k}$ | Laplacian quadratic penalty |
| Blind/myopic deconvolution [2202.13242], [1004.5538] | Image prior strength and PSF parameters jointly | SURE fixed points or Bayesian posterior inference |
| Null-space/IPSF methods [1502.07781], [1206.3594] | PSF, IPSF, and image estimate | Surface-area and curved-space regularization |
| Posterior PSF sampling [2511.19594] | 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 $\ell_1$ 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.

Source: https://www.emergentmind.com/topics/psf-regularization