---
title: Regularized Maximum Likelihood Imaging
url: https://www.emergentmind.com/topics/regularized-maximum-likelihood-imaging
type: topic
---

# Regularized Maximum Likelihood Imaging

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 [2404.05445] [1212.3258] [2105.11512] [1805.01712].

## 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 [2404.05445]. 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 [1710.02950].

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 [1212.3258]. 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 [2105.11512].

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 \((u,v)\) coordinates [2209.11813] [2504.19111]. In solar DEM inversion, the same Poisson logic yields
\[
\Phi(\boldsymbol{\xi})=\sum_i \left[(\mathbf{K}\boldsymbol{\xi})_i-g_i\log(\mathbf{K}\boldsymbol{\xi})_i\right]+\lambda\sum_j T_j\,\xi_j\,\Delta T_j,
\qquad \boldsymbol{\xi}\ge 0,
\]
which is a literal regularized maximum-likelihood objective over a thermal distribution rather than a conventional image [2301.04688].

## 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 | \(D+\lambda R\) | DEM inversion with \(\sum_j T_j\xi_j\Delta T_j\) |
| 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 \(g_\theta(x)\) | unsupervised marginal-likelihood training |

The solar DEM method is an explicit RML construction: it adds the linear temperature-weighted penalty
\[
R(\boldsymbol{\xi})=\sum_j T_j\,\xi_j\,\Delta T_j
\]
to a Poisson likelihood, thereby suppressing spurious hot tails while preserving positivity through multiplicative updates [2301.04688]. Optical phase retrieval with sparse regularization of amplitude and phase is also explicit: the objective combines a Gaussian intensity-fit term with \(\ell_1\) penalties on BM3D-frame coefficients,
\[
\sum_{r=1}^K \frac{1}{2\sigma_r^2}\|o_r-|u_r|^2\|_2^2+\tau_a\|\theta_a\|_1+\tau_\varphi\|\theta_\varphi\|_1,
\]
under propagation and synthesis-analysis constraints [1108.3251]. In interferometric imaging, MPoL and related workflows combine the Gaussian visibility likelihood with entropy, sparsity, total variation, or total squared variation; exoALMA used entropy, \(L_1\) sparsity, and TSV, while the ALMA continuum study found entropy + sparsity + TSV to be the preferred combination for ringed disks [2504.19111] [2209.11813].

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 \( |\mathcal{F}(X)+B|^2 \), beamstop masking inside the likelihood, finite-dimensional parameterization, limited iteration count, and reference design [2105.11512]. 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 [1805.01712].

Algorithmic regularization is another recurrent form. For constrained ML with Gaussian or Poisson noise, the constrained backprojected residual
\[
r(x,y)=\|x\,\nabla L_y(x)\|^2
\]
defines statistically justified stopping rules that regularize ISRA and EM without altering the objective itself [1212.3258]. 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 [1808.01072].

## 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,
\[
x^{(k+1)}=\frac{x^{(k)}}{H^T1}\,H^T\frac{y}{Hx^{(k)}},
\]
while the Gaussian analogue is ISRA,
\[
x^{(k+1)}=x^{(k)}\frac{H^Ty}{H^THx^{(k)}}.
\]
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 \(y\notin H(\mathcal C)\) [1212.3258]. Chretien and Hero reinterpret EM itself as a KL-proximal point algorithm:
\[
\theta^{k+1}=\arg\max_\theta\left\{l_y(\theta)-\beta_k I_y(\theta^k,\theta)\right\},
\]
with standard EM recovered at \(\beta_k=1\). When \(\beta_k\to 0\), 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 [1201.5907].

Low-photon holographic phase retrieval remains nonconvex because of the map \(X\mapsto |\mathcal{F}(X)+B|^2\), 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 [2105.11512]. 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,
\[
\theta=\mathfrak{Th}_\tau(u)=\operatorname{sign}(u)\circ(|u|-\tau)_+,
\]
followed by synthesis of a filtered object proxy and propagation updates [1108.3251].

In solar DEM inversion, the Poisson likelihood and linear temperature penalty produce a multiplicative RML update
\[
\xi^{(l+1)}=\frac{\xi^{(l)}}{\mathbf{K}^T\mathbf{1}+\lambda\mathbf{T}}\,\mathbf{K}^T\!\left(\frac{\mathbf g}{\mathbf K\xi^{(l)}}\right),
\qquad \xi^{(0)}=\mathbf 1,
\]
which preserves nonnegativity by construction [2301.04688]. 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 [2209.11813] [2504.19111].

A distinct optimization layer arises when the regularizer itself is learned. Unsupervised convex-regularizer training by marginal likelihood uses the gradient identity
\[
\nabla_\theta \log p(y\mid \theta)
=
\mathbb E_{x\mid \theta}[\nabla_\theta g_\theta(x)]
-
\mathbb E_{x\mid y,\theta}[\nabla_\theta g_\theta(x)],
\]
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 \(g_\theta\) in a conventional MAP reconstruction [2404.05445].

## 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 \(\lambda\) is selected by a Morozov discrepancy principle: start from a deliberately large \(\lambda\), reconstruct a DEM, compute reduced \(\chi^2\), decrease \(\lambda\) by a fixed factor—typically \(2/3\)—and stop when reduced \(\chi^2<1\). The paper emphasizes that \(\lambda\) should vary by pixel and that a global constant degrades performance substantially [2301.04688]. In ALMA continuum RML, model selection is performed by \(K\)-fold cross-validation on visibility cells,
\[
\mathrm{CV}=\sum_{k=1}^K \chi^2(\mathbf D_k\mid \mathbf I_{\mathrm{train}}),
\]
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 [2209.11813]. 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 [2504.19111].

Empirical Bayes supplies a more explicitly likelihood-based route. For convex imaging problems with prior
\[
p(x\mid \theta)=e^{-\theta^T g(x)}/Z(\theta),
\]
the hyperparameters are estimated by
\[
\hat\theta\in \arg\max_{\theta\in\Theta} p(y\mid \theta),
\qquad
p(y\mid \theta)=\int p(y\mid x)p(x\mid \theta)\,dx,
\]
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 [1911.11709]. The unsupervised convex-neural-prior work applies the same empirical-Bayes principle to learn the regularizer itself, again through marginal likelihood [2404.05445].

Stopping is itself a regularization parameter in iterative ML. The constrained backprojected residual criterion
\[
r(x,y)=\|x\,\nabla L_y(x)\|^2
\]
aligns the stopping decision with the KKT condition \(x\,\nabla L_y(x)=0\), and it is proved to yield bona fide regularization methods for ISRA and EM even when classical discrepancy criteria never trigger [1212.3258]. At the other end of the pipeline, nonconvex likelihood problems can be subjected to an explicit statistical hypothesis test
\[
H_0:\ \hat\theta=\hat\theta_{\mathrm{Global}}
\qquad\text{vs.}\qquad
H_1:\ \hat\theta\neq \hat\theta_{\mathrm{Global}},
\]
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 [1906.00101].

## 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 \( |\mathcal F(X)|^2 \) to \( |\mathcal F(X)+B|^2 \), allowing direct Poisson ML reconstruction in the low-photon regime. The study considers photon fluxes \(N_p=1000,100,10,1,0.1\), reports that HoloML-CG and HoloML-TR outperform inverse and Wiener filtering at \(N_p=1\), remain robust with a \(25\times 25\) beamstop, and support reduced oversampling down to about \(1.25\), 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 [2105.11512].

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 \(2048\times 2048\) pixels over a 12 nm field of view, with dwell time \(0.5\,\mu\mathrm{s}\), dose \(1.8\cdot 10^4\, e^- / \text{\AA}^2\) per frame, and post-processed single-frame SNR about \(-6.9\) 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 \(10.4\%\), \(0.5\%\), \(1.8\%\), and \(1.7\%\) [1805.01712].

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 \(\log_{10}T(\mathrm K)=5.85,6.00,6.15,\dots,7.50\), and evaluates performance with reduced \(\chi^2\) 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 \(500\times 500\) pixels on an Apple M1, excluding uncertainty estimation [2301.04688].

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 [2209.11813]. 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 [1910.00013]. In exoALMA, MPoL-based RML imaging of multi-channel \(^{12}\)CO, \(^{13}\)CO, 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 [2504.19111].

Other modalities reinforce the same pattern. Multi-plane optical phase retrieval with sparse BM3D-frame regularization reports phase RMSE \(0.036\) and amplitude RMSE \(0.026\), compared with \(0.26\) and \(0.23\) for the earlier AL algorithm and \(0.58\) and \(0.35\) for SBMIR in the reported phase-only simulation [1108.3251]. 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 \(0.8708\) to \(0.9964\) for W vs Air and from \(0.8366\) to \(0.9764\) for W vs Fe, while the abstract summarizes the discrimination gain as about 15% over PoCA [1806.01161].

## 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 [1808.01072].

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 [2105.11512] [1805.01712]. 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
\[
d(\widehat\Lambda)=E_{\Lambda^*}\log\!\left(\frac{f_{\Lambda^*}(X)}{f_{\widehat\Lambda}(X)}\right),
\]
that is, prediction accuracy in model space. The paper is explicit that such guarantees do not automatically imply image-space recovery in \(\ell_2\), PSNR, or perceptual quality [1710.02950].

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 [2504.19111]. 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 [2209.11813].

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 [1212.3258]. 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 [1906.00101].

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.

Source: https://www.emergentmind.com/topics/regularized-maximum-likelihood-imaging