Regularized Maximum Likelihood Imaging

• Regularized maximum likelihood imaging is a family of inverse-problem methods that infers unknown images by optimizing a likelihood function matched to measurement physics while incorporating structural or algorithmic regularization.
• It employs noise-specific data-fidelity terms, such as Gaussian least squares or Poisson Kullback–Leibler, combined with explicit penalties, learned priors, or stopping rules to stabilize the reconstruction.
• These techniques have been successfully applied in modalities like holographic phase retrieval, low-dose STEM, solar DEM inversion, and radio interferometry to improve reconstruction accuracy and robustness.

Regularized maximum likelihood imaging denotes a family of inverse-problem methods in which an unknown image, field, or latent spatial distribution is inferred by optimizing a likelihood matched to the measurement physics while stabilizing the reconstruction with penalties, constraints, structural priors, or algorithmic devices such as early stopping. In its most explicit form, the estimator is a penalized likelihood or MAP reconstruction; in a broader and historically important sense, the same designation also covers likelihood-driven reconstructions whose effective regularization comes from positivity, physical forward models, known references, symmetry constraints, latent template structure, or statistically calibrated stopping rules rather than from an additive penalty alone (Tan et al., 2024, Benvenuto et al., 2012, Barmherzig et al., 2021, Kramberger et al., 2018).

1. Statistical formulation and data-fidelity structure

The canonical RML problem is a likelihood-based variational optimization. In the Bayesian formulation used for unsupervised regularizer learning, reconstruction is written as

$$\hat x \in \arg\min_x \bigl\{D(Ax,y)+R_\theta(x)\bigr\},$$

with the data term identified as a negative log-likelihood and the regularizer as a negative log-prior (Tan et al., 2024). A closely related abstract formulation appears in the general theory of maximum regularized likelihood estimators: $$\widehat{\Lambda} \in \arg\min_{\Lambda \in \mathcal{L}} \left\{ -\log f_\Lambda(X) + r\,u(\Lambda) \right\},$$ where $\mathcal{L}$ is convex and the regularizer $u$ is definite and positively homogeneous (Zhuang et al., 2017).

Across imaging modalities, the data-fidelity term is determined by the acquisition statistics. For Gaussian models, constrained ML reduces to least squares,

$D_{LS}(y,x)=\|Hx-y\|^2,$

whereas for Poisson models it becomes the Kullback–Leibler discrepancy

$$D_{KL}(y,x)=\sum_{i=1}^N y_i \log\frac{y_i}{(Hx)_i} + (Hx)_i - y_i.$$

This distinction is central in photon-limited imaging, because Gaussian surrogates are mismatched when counts are low (Benvenuto et al., 2012). In holographic coherent diffraction imaging, the Poisson negative log-likelihood is applied directly to the intensity predicted by a reference-assisted Fourier model,

$$l(X)=\frac{1}{2}\sum_{(i,j)\in M}\left(\mathcal{B}_{ij}\,|\mathcal{F}(X)_{ij}+B_{ij}|^2 - Y_{ij}\log |\mathcal{F}(X)_{ij}+B_{ij}|^2\right),$$

with beamstop masking embedded in the likelihood rather than treated by interpolation (Barmherzig et al., 2021).

Interferometric RML uses the same pattern with Gaussian visibility noise. For ALMA and exoALMA imaging, the negative log-likelihood is

$$L_{\mathrm{nll}}(\mathbf{I})=\frac12\chi^2(\mathbf{D}\mid \mathbf{I}), \qquad \chi^2(\mathbf{D}\mid \mathbf{I}) = \sum_i \left|\frac{D_i-V_i(\mathbf{I})}{\sigma_i}\right|^2,$$

where $V_i(\mathbf I)$ is the model visibility from the image $\mathbf I$ at sampled $$\widehat{\Lambda} \in \arg\min_{\Lambda \in \mathcal{L}} \left\{ -\log f_\Lambda(X) + r\,u(\Lambda) \right\},$$0 coordinates (Zawadzki et al., 2022, Zawadzki et al., 27 Apr 2025). In solar DEM inversion, the same Poisson logic yields

$$\widehat{\Lambda} \in \arg\min_{\Lambda \in \mathcal{L}} \left\{ -\log f_\Lambda(X) + r\,u(\Lambda) \right\},$$1

which is a literal regularized maximum-likelihood objective over a thermal distribution rather than a conventional image (Massa et al., 2023).

2. What counts as regularization

A central feature of the literature is that “regularization” is not limited to an explicit additive penalty. Some reconstructions are explicitly penalized. Others are pure ML in objective form but effectively stabilized by physics, constraints, or structured latent models. The distinction matters, because several of the strongest empirical gains arise before any hand-crafted penalty is added.

Mechanism	Representative structure	Example
Explicit penalty	$$\widehat{\Lambda} \in \arg\min_{\Lambda \in \mathcal{L}} \left\{ -\log f_\Lambda(X) + r\,u(\Lambda) \right\},$$2	DEM inversion with $$\widehat{\Lambda} \in \arg\min_{\Lambda \in \mathcal{L}} \left\{ -\log f_\Lambda(X) + r\,u(\Lambda) \right\},$$3
Structural prior	known reference, masking, finite parameterization	holographic CDI
Latent model restriction	templates, symmetry averaging, model splitting	low-dose graphene STEM
Algorithmic regularization	stopping rule or inward hedge	CBR stopping; hedged tomography
Learned prior	convex neural $$\widehat{\Lambda} \in \arg\min_{\Lambda \in \mathcal{L}} \left\{ -\log f_\Lambda(X) + r\,u(\Lambda) \right\},$$4	unsupervised marginal-likelihood training

The solar DEM method is an explicit RML construction: it adds the linear temperature-weighted penalty

$$\widehat{\Lambda} \in \arg\min_{\Lambda \in \mathcal{L}} \left\{ -\log f_\Lambda(X) + r\,u(\Lambda) \right\},$$5

to a Poisson likelihood, thereby suppressing spurious hot tails while preserving positivity through multiplicative updates (Massa et al., 2023). Optical phase retrieval with sparse regularization of amplitude and phase is also explicit: the objective combines a Gaussian intensity-fit term with $$\widehat{\Lambda} \in \arg\min_{\Lambda \in \mathcal{L}} \left\{ -\log f_\Lambda(X) + r\,u(\Lambda) \right\},$$6 penalties on BM3D-frame coefficients,

$$\widehat{\Lambda} \in \arg\min_{\Lambda \in \mathcal{L}} \left\{ -\log f_\Lambda(X) + r\,u(\Lambda) \right\},$$7

under propagation and synthesis-analysis constraints (Migukin et al., 2011). In interferometric imaging, MPoL and related workflows combine the Gaussian visibility likelihood with entropy, sparsity, total variation, or total squared variation; exoALMA used entropy, $$\widehat{\Lambda} \in \arg\min_{\Lambda \in \mathcal{L}} \left\{ -\log f_\Lambda(X) + r\,u(\Lambda) \right\},$$8 sparsity, and TSV, while the ALMA continuum study found entropy + sparsity + TSV to be the preferred combination for ringed disks (Zawadzki et al., 27 Apr 2025, Zawadzki et al., 2022).

By contrast, holographic phase retrieval in the low-photon regime is explicitly characterized as pure maximum likelihood rather than MAP. There is no additive penalty such as TV, sparsity, smoothness, or nonnegativity in the HoloML objective. Its effective stabilization comes from the known reference object, the physically informed forward model $$\widehat{\Lambda} \in \arg\min_{\Lambda \in \mathcal{L}} \left\{ -\log f_\Lambda(X) + r\,u(\Lambda) \right\},$$9, beamstop masking inside the likelihood, finite-dimensional parameterization, limited iteration count, and reference design (Barmherzig et al., 2021). Low-dose STEM defect reconstruction is similar: there is no TV or quadratic penalty, but the reconstruction is heavily constrained by a finite mixture of recurring template images, a registered periodic graphene lattice, fixed supercell size, exact symmetry handling on a hexagonal grid, iterative model cloning, and a stopping rule based on whether new meaningful structures continue to emerge (Kramberger et al., 2018).

Algorithmic regularization is another recurrent form. For constrained ML with Gaussian or Poisson noise, the constrained backprojected residual

$\mathcal{L}$0

defines statistically justified stopping rules that regularize ISRA and EM without altering the objective itself (Benvenuto et al., 2012). In quantum state tomography, an opposite lesson appears: plain likelihood optimization is often too willing to hit the boundary of the feasible set. The proposed cure is “hedging” away from pure, rank-deficient states, for example by moving a boundary estimate slightly toward the interior; this is a form of regularization aimed at reducing overconfidence rather than enforcing spatial smoothness (Ferrie et al., 2018).

3. Optimization methodologies

Because RML objectives inherit nonlinearity from the forward model and non-smoothness from the regularizer, optimization is usually modality-specific. Yet several recurring patterns appear: multiplicative EM-type schemes for Poisson models, proximal-point and trust-region reinterpretations of likelihood ascent, augmented Lagrangian splitting, direct pixel optimization in visibility space, and stochastic approximation for empirical-Bayes training.

For Poisson inverse problems, the classical constrained ML update is EM / Richardson–Lucy,

$\mathcal{L}$1

while the Gaussian analogue is ISRA,

$\mathcal{L}$2

What the stopping-rule work adds is not a new objective, but a proof that these iterative solvers become genuine regularization methods when stopped by the CBR criterion rather than by discrepancy rules that can fail when $\mathcal{L}$3 (Benvenuto et al., 2012). Chretien and Hero reinterpret EM itself as a KL-proximal point algorithm: $\mathcal{L}$4 with standard EM recovered at $\mathcal{L}$5. When $\mathcal{L}$6, the proximal damping is progressively relaxed and the convergence becomes superlinear under the paper’s assumptions; the practical implementation uses a trust-region strategy in which the relaxation parameter is the dual variable for a KL-defined trust region (Chrétien et al., 2012).

Low-photon holographic phase retrieval remains nonconvex because of the map $\mathcal{L}$7, but the paper derives explicit Wirtinger gradients for complex and real-valued objects and reports that conjugate gradient and trust-region methods converge to almost identical high-quality solutions. The two implementations are HoloML-CG and HoloML-TR; no EM, ADMM, proximal splitting, or majorization–minimization scheme is introduced there (Barmherzig et al., 2021). Optical multi-plane phase retrieval with sparse regularization uses a decoupled augmented Lagrangian, splitting the data/propagation part from the sparse coefficient part, with soft-thresholding of BM3D-frame coefficients,

$\mathcal{L}$8

followed by synthesis of a filtered object proxy and propagation updates (Migukin et al., 2011).

In solar DEM inversion, the Poisson likelihood and linear temperature penalty produce a multiplicative RML update

$\mathcal{L}$9

which preserves nonnegativity by construction (Massa et al., 2023). Interferometric MPoL workflows instead parameterize the image directly in pixel space and minimize the differentiable loss with gradient-based optimization on GPU-accelerated PyTorch pipelines; exoALMA similarly uses direct pixel optimization in MPoL after gridding visibilities to Fourier cells (Zawadzki et al., 2022, Zawadzki et al., 27 Apr 2025).

A distinct optimization layer arises when the regularizer itself is learned. Unsupervised convex-regularizer training by marginal likelihood uses the gradient identity

$u$0

estimated from prior and posterior Langevin chains. The resulting SAPG-ULA scheme updates the regularizer parameters without ground-truth images and later uses the trained $u$1 in a conventional MAP reconstruction (Tan et al., 2024).

4. Parameter selection, stopping, and validation

RML imaging depends not only on the objective and optimizer, but also on how regularization strength, stopping time, and reconstruction validity are chosen. The literature spans discrepancy principles, cross-validation, marginal likelihood, statistical stopping rules, and ex post tests of global optimality.

In the DEM setting, the regularization weight $u$2 is selected by a Morozov discrepancy principle: start from a deliberately large $u$3, reconstruct a DEM, compute reduced $u$4, decrease $u$5 by a fixed factor—typically $u$6—and stop when reduced $u$7. The paper emphasizes that $u$8 should vary by pixel and that a global constant degrades performance substantially (Massa et al., 2023). In ALMA continuum RML, model selection is performed by $u$9-fold cross-validation on visibility cells,

$D_{LS}(y,x)=\|Hx-y\|^2,$0

with both random-cell and dartboard partitions studied. The principal conclusion is not that CV identifies a unique optimum, but that it reliably identifies a broad family of models with comparably strong predictive power (Zawadzki et al., 2022). ExoALMA extends this logic to multi-channel spectral cubes and finds that hyperparameters tuned on one representative channel are usually stable across non-adjacent channels, continuum-subtracted and non-continuum-subtracted cubes, and even across multiple lines of the same source, which makes cube-wide RML imaging computationally tractable (Zawadzki et al., 27 Apr 2025).

Empirical Bayes supplies a more explicitly likelihood-based route. For convex imaging problems with prior

$D_{LS}(y,x)=\|Hx-y\|^2,$1

the hyperparameters are estimated by

$D_{LS}(y,x)=\|Hx-y\|^2,$2

and then plugged into the MAP reconstruction. The SAPG methodology estimates one or several regularization parameters directly from the data and was demonstrated for denoising, deconvolution, hyperspectral unmixing, and TGV-like models (Vidal et al., 2019). The unsupervised convex-neural-prior work applies the same empirical-Bayes principle to learn the regularizer itself, again through marginal likelihood (Tan et al., 2024).

Stopping is itself a regularization parameter in iterative ML. The constrained backprojected residual criterion

$D_{LS}(y,x)=\|Hx-y\|^2,$3

aligns the stopping decision with the KKT condition $D_{LS}(y,x)=\|Hx-y\|^2,$4, and it is proved to yield bona fide regularization methods for ISRA and EM even when classical discrepancy criteria never trigger (Benvenuto et al., 2012). At the other end of the pipeline, nonconvex likelihood problems can be subjected to an explicit statistical hypothesis test

$D_{LS}(y,x)=\|Hx-y\|^2,$5

with a reparameterized embedding used to improve power. In the camera-blur application, this global-optimum test is used as a stopping rule for multistart local optimization (LeBlanc et al., 2019).

5. Representative imaging modalities

The breadth of RML imaging is best seen across applications. In holographic coherent diffraction imaging, a known reference object modifies the forward model from $D_{LS}(y,x)=\|Hx-y\|^2,$6 to $D_{LS}(y,x)=\|Hx-y\|^2,$7, allowing direct Poisson ML reconstruction in the low-photon regime. The study considers photon fluxes $D_{LS}(y,x)=\|Hx-y\|^2,$8, reports that HoloML-CG and HoloML-TR outperform inverse and Wiener filtering at $D_{LS}(y,x)=\|Hx-y\|^2,$9, remain robust with a $$D_{KL}(y,x)=\sum_{i=1}^N y_i \log\frac{y_i}{(Hx)_i} + (Hx)_i - y_i.$$0 beamstop, and support reduced oversampling down to about $$D_{KL}(y,x)=\sum_{i=1}^N y_i \log\frac{y_i}{(Hx)_i} + (Hx)_i - y_i.$$1, zero specimen-reference separation, and irregular references such as an annulus; URA references consistently give the best reconstructions, while block references are best among simple geometries (Barmherzig et al., 2021).

Low-dose STEM of graphene provides a different paradigm: the unknown is not a single denoised frame but a mixture of recurring defect templates on a periodic lattice. The dataset contains 1187 MAADF frames of $$D_{KL}(y,x)=\sum_{i=1}^N y_i \log\frac{y_i}{(Hx)_i} + (Hx)_i - y_i.$$2 pixels over a 12 nm field of view, with dwell time $$D_{KL}(y,x)=\sum_{i=1}^N y_i \log\frac{y_i}{(Hx)_i} + (Hx)_i - y_i.$$3, dose $$D_{KL}(y,x)=\sum_{i=1}^N y_i \log\frac{y_i}{(Hx)_i} + (Hx)_i - y_i.$$4 per frame, and post-processed single-frame SNR about $$D_{KL}(y,x)=\sum_{i=1}^N y_i \log\frac{y_i}{(Hx)_i} + (Hx)_i - y_i.$$5 dB. Reconstruction proceeds by likelihood maximization over model images, symmetry states, and weights using an empirical Gamma noise model, eventually recovering four archetypal vacancy defects with reported aggregate weights $$D_{KL}(y,x)=\sum_{i=1}^N y_i \log\frac{y_i}{(Hx)_i} + (Hx)_i - y_i.$$6, $$D_{KL}(y,x)=\sum_{i=1}^N y_i \log\frac{y_i}{(Hx)_i} + (Hx)_i - y_i.$$7, $$D_{KL}(y,x)=\sum_{i=1}^N y_i \log\frac{y_i}{(Hx)_i} + (Hx)_i - y_i.$$8, and $$D_{KL}(y,x)=\sum_{i=1}^N y_i \log\frac{y_i}{(Hx)_i} + (Hx)_i - y_i.$$9 (Kramberger et al., 2018).

Solar differential emission measure inversion is a pixel-wise thermal imaging problem. The RML method uses five AIA EUV channels—94, 131, 171, 193, and 211 Å—while excluding 335 Å, discretizes the DEM into 12 temperature bins centered at $$l(X)=\frac{1}{2}\sum_{(i,j)\in M}\left(\mathcal{B}_{ij}\,|\mathcal{F}(X)_{ij}+B_{ij}|^2 - Y_{ij}\log |\mathcal{F}(X)_{ij}+B_{ij}|^2\right),$$0, and evaluates performance with reduced $$l(X)=\frac{1}{2}\sum_{(i,j)\in M}\left(\mathcal{B}_{ij}\,|\mathcal{F}(X)_{ij}+B_{ij}|^2 - Y_{ij}\log |\mathcal{F}(X)_{ij}+B_{ij}|^2\right),$$1 and NRMSE over 25 Poisson realizations. Its main empirical claim is that the temperature-weighted regularizer suppresses spurious high-temperature tails that unregularized ML frequently invents, while still recovering a true high-temperature component when present; the implementation runs in about 35 s for $$l(X)=\frac{1}{2}\sum_{(i,j)\in M}\left(\mathcal{B}_{ij}\,|\mathcal{F}(X)_{ij}+B_{ij}|^2 - Y_{ij}\log |\mathcal{F}(X)_{ij}+B_{ij}|^2\right),$$2 pixels on an Apple M1, excluding uncertainty estimation (Massa et al., 2023).

Radio interferometric imaging has become a major modern RML domain. For ALMA continuum observations of HD 143006, MPoL reconstructs images directly from uniformly weighted gridded visibilities and reports spatial-resolution improvement by up to a factor of 3 without sacrificing sensitivity; the recommended workflow uses entropy, sparsity, and TSV, with cross-validation over coarse grids of regularization strengths (Zawadzki et al., 2022). For simulated ngVLA stellar radio photospheres, RML in SMILI is compared with multi-scale CLEAN and achieves better goodness-of-fit to the data, lower residual errors, and better recovery of representative structures for most stellar models, especially when the synthesized beam is highly non-Gaussian (Akiyama et al., 2019). In exoALMA, MPoL-based RML imaging of multi-channel $$l(X)=\frac{1}{2}\sum_{(i,j)\in M}\left(\mathcal{B}_{ij}\,|\mathcal{F}(X)_{ij}+B_{ij}|^2 - Y_{ij}\log |\mathcal{F}(X)_{ij}+B_{ij}|^2\right),$$3CO, $$l(X)=\frac{1}{2}\sum_{(i,j)\in M}\left(\mathcal{B}_{ij}\,|\mathcal{F}(X)_{ij}+B_{ij}|^2 - Y_{ij}\log |\mathcal{F}(X)_{ij}+B_{ij}|^2\right),$$4CO, and CS cubes independently reproduces the non-Keplerian features seen in CLEAN images, and is used specifically as an independent verification of marginal structures rather than as a replacement for CLEAN (Zawadzki et al., 27 Apr 2025).

Other modalities reinforce the same pattern. Multi-plane optical phase retrieval with sparse BM3D-frame regularization reports phase RMSE $$l(X)=\frac{1}{2}\sum_{(i,j)\in M}\left(\mathcal{B}_{ij}\,|\mathcal{F}(X)_{ij}+B_{ij}|^2 - Y_{ij}\log |\mathcal{F}(X)_{ij}+B_{ij}|^2\right),$$5 and amplitude RMSE $$l(X)=\frac{1}{2}\sum_{(i,j)\in M}\left(\mathcal{B}_{ij}\,|\mathcal{F}(X)_{ij}+B_{ij}|^2 - Y_{ij}\log |\mathcal{F}(X)_{ij}+B_{ij}|^2\right),$$6, compared with $$l(X)=\frac{1}{2}\sum_{(i,j)\in M}\left(\mathcal{B}_{ij}\,|\mathcal{F}(X)_{ij}+B_{ij}|^2 - Y_{ij}\log |\mathcal{F}(X)_{ij}+B_{ij}|^2\right),$$7 and $$l(X)=\frac{1}{2}\sum_{(i,j)\in M}\left(\mathcal{B}_{ij}\,|\mathcal{F}(X)_{ij}+B_{ij}|^2 - Y_{ij}\log |\mathcal{F}(X)_{ij}+B_{ij}|^2\right),$$8 for the earlier AL algorithm and $$l(X)=\frac{1}{2}\sum_{(i,j)\in M}\left(\mathcal{B}_{ij}\,|\mathcal{F}(X)_{ij}+B_{ij}|^2 - Y_{ij}\log |\mathcal{F}(X)_{ij}+B_{ij}|^2\right),$$9 and $$L_{\mathrm{nll}}(\mathbf{I})=\frac12\chi^2(\mathbf{D}\mid \mathbf{I}), \qquad \chi^2(\mathbf{D}\mid \mathbf{I}) = \sum_i \left|\frac{D_i-V_i(\mathbf{I})}{\sigma_i}\right|^2,$$0 for SBMIR in the reported phase-only simulation (Migukin et al., 2011). Muon tomography from multiple Coulomb scattering formulates a Gaussian-like scattering likelihood over voxel scattering densities and uses an MLS-EM update rather than PoCA’s single-scatter assignment; in simulation, ROC AUC rises from $$L_{\mathrm{nll}}(\mathbf{I})=\frac12\chi^2(\mathbf{D}\mid \mathbf{I}), \qquad \chi^2(\mathbf{D}\mid \mathbf{I}) = \sum_i \left|\frac{D_i-V_i(\mathbf{I})}{\sigma_i}\right|^2,$$1 to $$L_{\mathrm{nll}}(\mathbf{I})=\frac12\chi^2(\mathbf{D}\mid \mathbf{I}), \qquad \chi^2(\mathbf{D}\mid \mathbf{I}) = \sum_i \left|\frac{D_i-V_i(\mathbf{I})}{\sigma_i}\right|^2,$$2 for W vs Air and from $$L_{\mathrm{nll}}(\mathbf{I})=\frac12\chi^2(\mathbf{D}\mid \mathbf{I}), \qquad \chi^2(\mathbf{D}\mid \mathbf{I}) = \sum_i \left|\frac{D_i-V_i(\mathbf{I})}{\sigma_i}\right|^2,$$3 to $$L_{\mathrm{nll}}(\mathbf{I})=\frac12\chi^2(\mathbf{D}\mid \mathbf{I}), \qquad \chi^2(\mathbf{D}\mid \mathbf{I}) = \sum_i \left|\frac{D_i-V_i(\mathbf{I})}{\sigma_i}\right|^2,$$4 for W vs Fe, while the abstract summarizes the discrimination gain as about 15% over PoCA (Wang et al., 2018).

6. Theory, caveats, and recurrent misconceptions

Several recurrent misconceptions are explicitly contradicted by the literature. First, plain maximum likelihood is not automatically optimal. In quantum state tomography, MLE is proved inadmissible for fidelity, squared Hilbert–Schmidt distance, and relative entropy, and the same mechanism extends to constrained least squares and nuclear-norm-based estimators whenever they can output pure or rank-deficient states. The practical message is that boundary solutions may be statistically overconfident, and small hedges toward the interior can uniformly improve risk (Ferrie et al., 2018).

Second, “regularized” does not always mean “explicitly penalized.” Holographic low-photon phase retrieval and low-dose STEM defect reconstruction are both central to the subject precisely because they show how much stabilization can come from the correct likelihood, acquisition design, and latent structure before TV, sparsity, or quadratic penalties are introduced (Barmherzig et al., 2021, Kramberger et al., 2018). A plausible implication is that likelihood mismatch can dominate penalty choice in photon-limited or structurally constrained problems.

Third, theoretical guarantees are often predictive rather than pixelwise. The general MRLE theory yields finite-sample oracle inequalities in KL divergence under convex parametrization and positively homogeneous regularization, without restricted eigenvalue or similar conditions. What it controls is

$$L_{\mathrm{nll}}(\mathbf{I})=\frac12\chi^2(\mathbf{D}\mid \mathbf{I}), \qquad \chi^2(\mathbf{D}\mid \mathbf{I}) = \sum_i \left|\frac{D_i-V_i(\mathbf{I})}{\sigma_i}\right|^2,$$5

that is, prediction accuracy in model space. The paper is explicit that such guarantees do not automatically imply image-space recovery in $$L_{\mathrm{nll}}(\mathbf{I})=\frac12\chi^2(\mathbf{D}\mid \mathbf{I}), \qquad \chi^2(\mathbf{D}\mid \mathbf{I}) = \sum_i \left|\frac{D_i-V_i(\mathbf{I})}{\sigma_i}\right|^2,$$6, PSNR, or perceptual quality (Zhuang et al., 2017).

Fourth, low residuals can be misleading. In exoALMA, the very low RMS of RML images is attributed to sparsity suppressing unsupported background pixels, not to an actual improvement in thermal sensitivity, and the native RML image has no unique restoring beam or straightforward effective resolution (Zawadzki et al., 27 Apr 2025). The ALMA continuum study makes a related point by showing that cross-validation minima are broad rather than sharply localized; model families, not single hyperparameter values, are often the stable object of inference (Zawadzki et al., 2022).

Finally, classical stopping criteria and local optimization heuristics are not universally reliable. For constrained ML, discrepancy rules such as Morozov, Pearson, or Poisson discrepancy can fail entirely when data lie outside the nonnegative range of the forward operator, whereas the CBR rule remains valid because it is tied to the KKT residual rather than to unattainable zero discrepancy (Benvenuto et al., 2012). For nonconvex likelihood landscapes, local convergence alone does not certify global optimality; the reparameterized-embedding test shows that this question can itself be turned into a statistical inference problem (LeBlanc et al., 2019).

Taken together, these results define regularized maximum likelihood imaging less as a single algorithm than as a design principle: specify the measurement law and forward physics accurately, choose the weakest regularization that resolves ill-posedness in the relevant geometry, optimize with algorithms matched to the likelihood structure, and validate both the hyperparameters and the attained solution statistically rather than heuristically.