Non-Linear Iterative Richardson-Lucy (NIRL)

• The NIRL algorithm is an advanced deconvolution method that iteratively incorporates non-linear corrections to recover true signals in complex data.
• It updates the forward operator at each iteration using effective kernel corrections and regularization to ensure convergence and physical positivity.
• Empirical tests in CMB delensing and image restoration demonstrate reduced reconstruction error and superior feature preservation compared to standard RL methods.

The Non-Linear Iterative Richardson-Lucy (NIRL) algorithm is an extension of the classic Richardson-Lucy (RL) deconvolution method, designed to address cases where the measurement process introduces non-linear or non-local effects, such as non-linear convolutions, weak lensing, or Bayesian spectral priors. Unlike standard RL algorithms, which assume a linear convolution operator, NIRL incorporates non-linear corrections iteratively, updating the forward operator or introducing additional regularization or prior terms based on the current estimate of the solution. This approach has found particular utility in astrophysical data analysis (e.g., Cosmic Microwave Background delensing), image restoration, and applications where higher fidelity recovery of structure is required without explicit prior assumptions.

1. Mathematical Foundations and Model Structure

In the standard RL framework, the observed blurred data $g$ is modeled as a linear convolution of the true signal $f$ via a known, positive-definite kernel or operator $H$: $g = Hf$. The RL update is multiplicative, updating $f$ via

$$f^{(k+1)} = f^{(k)} \odot H^T \left(\frac{g}{H f^{(k)}}\right),$$

where $\odot$ is element-wise multiplication and the ratio is taken pointwise. This iteration finds maximum-likelihood Poisson solutions while enforcing non-negativity.

NIRL generalizes this update to cases where the measurement process is non-linear, and the forward operator depends on the current estimate. This structure arises, for example, in CMB delensing, where the observed lensed spectrum $\widetilde{C}_\ell^{TT}$ is not a linear convolution of the primordial power spectrum $P(k)$, but involves higher-order corrections:

$$\widetilde{C}_\ell^{TT} = \sum_k P(k) \left[ G_{\ell k}^{TT} + \Delta G_{\ell k}^{TT}(P) \right],$$

where $\Delta G_{\ell k}^{TT}$ is a small, but $P$-dependent, non-linear kernel correction representing lensing effects (Chandra et al., 2021). Analogously, NIRL can implement prior-based or spectral corrections, as in Dirichlet-prior Bayesian MAP estimation for image deblurring (Hendrix et al., 2024).

2. NIRL Iterative Scheme and Algorithmic Steps

The NIRL approach is defined by an iteration in which, at each step, the forward operator or correction term is updated using the most recent estimate of the latent field:

1. Kernel Evaluation: Compute the effective, possibly non-linear kernel $\widetilde{G}^{(i)}$ or forward operator at iteration $i$ based on $P^{(i)}$ or $f^{(i)}$.
2. Forward Projection: Calculate the forward-model spectrum or blurred image.
3. Residual Calculation: Form the ratio between observed data and the current modeled prediction.
4. Update Rule: Multiply the current estimate by an RL-type correction using the current effective kernel.
5. Regularization or Priors: Optionally, incorporate regularization or spectral prior terms multiplicatively or additively, yielding a fully non-linear iteration.
6. Convergence Assessment: Check stopping criteria based on fractional changes, $\chi^2$ plateaus, or relative error metrics.

In CMB applications, the effective kernel update is given by

$$\widetilde{G}_{\ell k}^{TT(i)} = G_{\ell k}^{TT} \left(1 + \sum_{k_1} P^{(i)}(k_1) \sum_L G_{L k_1}^{\phi \phi} S^{(b)}_{\ell L} \right) + \sum_{k_1} P^{(i)}(k_1) \sum_{L,\ell'} G_{L k_1}^{\phi \phi} G_{\ell' k}^{TT} S^{(a)}_{\ell L \ell'},$$

and the primary NIRL update is

$$P^{(i+1)}(k) = P^{(i)}(k) \sum_\ell \widetilde{G}_{\ell k}^{TT(i)} \frac{\hat{C}_\ell^{TT}}{C_\ell^{TT(i)}},$$

with $$C_\ell^{TT(i)} = \sum_k \widetilde{G}_{\ell k}^{TT(i)} P^{(i)}(k)$$ (Chandra et al., 2021).

For Bayesian MAP NIRL (Hendrix et al., 2024), the update is

$$x^{(k+1)} = x^{(k)} \circ \frac{H^T \big(y / (H x^{(k)}) \big) + \lambda F^{-1}\left(\frac{\alpha-1}{F x^{(k)}}\right)}{1 + \lambda F^{-1}(\mathbf{1})},$$

where $F$ and $F^{-1}$ are the Fourier transform and its inverse and $\alpha$ specifies the spectral Dirichlet prior.

3. Theoretical Properties and Convergence

NIRL maintains the positivity and multiplicative structure characteristic of RL schemes, which is crucial for interpreting recovered spectra and images physically. Under small non-linearity assumptions for the perturbative correction (e.g., weak lensing), convergence to the true solution is guaranteed to the same order as in the standard RL by Neumann series expansion (Chandra et al., 2021). The multiplicative correction avoids the imposition of smoothness priors, preserving the ability to reconstruct sharp features provided they are above the noise or cosmic-variance threshold.

Stopping criteria include:

• Fixed number of iterations,
• Maximal fractional change threshold in the estimate,
• Plateauing of fit residuals or $\chi^2$.

Formal guarantees are preserved when the non-linear correction is sufficiently small and the underlying forward operator remains positive-definite.

4. Performance and Comparison to Alternative Deconvolution Approaches

Empirical studies in CMB PPS reconstruction demonstrate that:

• Linear RL on lensed data produces spurious high-frequency oscillations ("saw-tooth artifacts").
• Template-based delensing approaches with a fixed power-law prior exactly recover power-law PPS, but incur residual errors for non-power-law inputs.
• NIRL, making no prior assumptions on $P(k)$, achieves robust, high-fidelity reconstructions. Mean absolute error is reduced by $\sim$2–10% relative to template-based RL for Gaussian bump and wavepacket feature test cases, saturating the ideal performance (Chandra et al., 2021).

For Bayesian imaging, NIRL-MAP outperforms both classic RL and RL+TV:

• On test patterns: RL PSNR peaks and then degrades with iteration; NIRL-MAP achieves higher, stable PSNR (e.g., 25.2 dB vs. 23.8 dB for Siemens stars at optimal iteration count).
• SSIM and artifact suppression (e.g., removal of ring and staircasing structures) are markedly superior for NIRL-MAP (Hendrix et al., 2024).

5. Algorithmic Implementation and Computational Cost

A simplified pseudocode for the CMB-delensing NIRL algorithm is as follows (Chandra et al., 2021):

input: observed lensed spectrum Ĉ[ℓ], unlensed kernel G_TT[ℓ,k], lensing kernel G_pp[L,k], shape factors S(a)[ℓ,L,ℓ′], S(b)[ℓ,L] initialize P[k] = P₀(k) # e.g., smooth power law for i = 0 … N_iter-1 do: # 1) compute lensing correction ΔG for ℓ, k: ΔG[ℓ,k] = G_TT[ℓ,k] * sum_{k1,L} P[k1]*G_pp[L,k1]*S(b)[ℓ,L] + sum_{k1,L,ℓ′} P[k1]*G_pp[L,k1]*G_TT[ℓ′,k]*S(a)[ℓ,L,ℓ′] G_eff[ℓ,k] = G_TT[ℓ,k] + ΔG[ℓ,k] # 2) forward model for ℓ: C_model[ℓ] = sum_k G_eff[ℓ,k]*P[k] # 3) RL update for k: correction = sum_ℓ G_eff[ℓ,k] * ( Ĉ[ℓ]/C_model[ℓ] ) P_new[k] = P[k] * correction P[:] = P_new[:] end for output: P_rec[k]

Computational overhead per iteration is 2–4x that of template delensing plus linear RL, due to the non-linear kernel update step (Chandra et al., 2021). For Bayesian image NIRL (Hendrix et al., 2024), the cost is dominated by FFTs and remains $\mathcal{O}(N \log N)$ per iteration, similar to linear RL, with only minor additional computational requirements.

6. Extensions, Regularization, and Related Non-Linear Iterative Schemes

NIRL is related to a broader class of non-linear iterative deconvolution methods:

• Edge-preserving or local regularizers are incorporated multiplicatively or through auxiliary fields to control amplification of small singular-value modes, as in RL with edge-preserving regularization (Mamba et al., 2024).
• Non-local acceleration schemes, such as the Dyson–Albertini algorithm (Dyson et al., 2015), introduce non-linear, non-local Laplacian corrections to RL. These methods perform an RL-type back-projection followed by a Laplacian sharpening, with parameters updated via $L^2$ line-search for optimal convergence, leading to improved high-frequency recovery and accelerated convergence in the noise-free setting.

These extensions maintain the Cauchy sequence nature of the RL algorithm, delivering fast and stable recovery of structure across a range of forward operators, and can be adapted to domain-specific requirements such as CMB delensing, physically motivated Bayesian imaging, or non-local image restoration.

7. Applications, Impact, and Empirical Validation

NIRL has been instrumental in making primordial power spectrum reconstruction from lensed CMB spectra robust to weak lensing contamination, without relying on power-law priors and while accurately capturing features that would otherwise be lost or biased (Chandra et al., 2021). Bayesian NIRL variants achieve statistically principled deconvolutions in photon-limited imaging, yielding improved signal-to-noise performance, artifact suppression, and direct uncertainty quantification (Hendrix et al., 2024).

A key impact of the NIRL framework is its ability to accommodate the physical and statistical structure of the measurement process without a priori smoothness or template assumptions, while remaining numerically tractable and convergence-guaranteed. Empirical results across domains report improved reconstruction error, enhanced robustness to non-linearities, and superior feature preservation compared to linear or ad hoc regularized RL.