PaperlayerView exact revision on arXiv

Source-linked AI summary

Uncertainty Quantification for Landweber Iteration with Randomized Normal-Operator Approximation

Anuj Abhishek, Sean Holman

arXiv:2610.00882v1math.NAcs.NAmath.STstat.TH

TL;DR

Uncertainty quantification for reconstructions from noisy, ill-posed inverse problems has received less attention than convergence analysis. This paper propagates observational noise through finite-step Landweber iteration and combines it with a bias bound; simulations report increased mean interior coverage at the lowest noise level when randomized SVD rank increases.

  • Problem

    Quantifying how random observation noise affects iterative reconstructions in linear ill-posed inverse problems has received less attention than studying convergence.

    Source: §1 · pp. 1–2

  • Method

    The method propagates observational noise through a finite-step Landweber map and combines its uncertainty with a regularization-bias bound to construct confidence regions.

    Source: §1 · p. 2, §3.3 · p. 7, §3.4 · p. 7

  • Results

    At σ = 0.1, mean interior coverage increased from 0.8740 at rank 500 to 0.9311 at rank 1300.

    Source: §4.2 · p. 16

  • Takeaways & Limitations

    Randomized SVD rank should depend on observational noise: lower noise calls for sufficient rank to limit truncation-induced undercoverage, while higher noise may need a smaller effective rank for coverage near 0.95.

    Source: §4.2 · p. 19

    The appropriate randomized SVD rank varies with noise level, and the results motivate developing an adaptive rank choice.

    Source: §4.2 · p. 19

Abstract

from arXiv · show

Iterative methods are extremely popular for solving linear ill-posed inverse problems. While deterministic convergence, or more precisely, semi-convergence of such methods is widely studied, the problem of quantifying uncertainty in the resulting reconstructions due to random noise in the observed data has received far less attention. In this work, we study uncertainty quantification for the Landweber iteration in linear inverse problems using the Radon transform as a motivating example. We interpret the Landweber reconstruction statistically by analyzing how uncertainty in the measured data is propagated through the reconstruction map. This perspective is closely related to generalized fiducial inference, where uncertainty about the parameter is induced by inverting the relation between the observed data and the unknown (fixed) quantity. We combine the resulting stochastic uncertainty with a theoretical bound on the regularization bias to construct confidence intervals for the reconstructed solution. The evaluation of the resulting uncertainty estimates require repeated operations with the forward operator that can become expensive in large-scale problems. To reduce this cost, we use randomized singular value decomposition to obtain a low-rank approximation of the normal operator associated with the Radon transform and incorporate this approximation in quantifying uncertainty. {{In particular, this requires that the additional randomization error introduced due to the use of such randomized techniques be included in our approach for uncertainty quantification.}} Numerical results show that the proposed method provides accurate reconstructions and reliable uncertainty estimates at substantially reduced computational cost.

1 Introduction

The paper studies uncertainty in reconstructions from noisy, ill-posed inverse problems, where exact inversion is unstable. It propagates observational noise through a finite-step Landweber map and combines the resulting uncertainty with a bias bound to construct confidence regions.

1 Introduction

  • The method is motivated by uncertainty quantification for linear ill-posed inverse problems, with computerized tomography and the Radon transform as an example.The observed data contain Gaussian noise whose distribution is known or can be estimated reliably.Source: §1 · p. 1
  • The approach is closely related to generalized fiducial inference, transferring data randomness to the parameter space without a Bayesian prior.The paper distinguishes this perspective from both frequentist estimators and Bayesian posterior distributions.Source: §1 · p. 2
  • The finite-step Landweber map produces a data-dependent distribution at each iteration, replacing unstable exact inversion with regularized inversion for uncertainty quantification.The stopping index controls both regularization and propagated uncertainty; early stopping limits noise amplification but introduces bias that must be estimated.Source: §1 · p. 2
  • Randomized SVD approximates the normal operator to reduce computational overhead and add regularization, requiring revised bias and variance estimates for confidence regions.The paper describes a new stopping rule based on randomized-SVD bias and variance approximations and provides lower bounds on frequentist coverage.Source: §1 · p. 3
  • The article presents the mathematical setup and method, then uses numerical simulations to check the proposed approach.Auxiliary analysis and discussion are also included.Source: §1 · p. 3

2 Observation model

The observation model represents noisy Radon-transform measurements of an unknown image on a fixed grid. The image and measurements are discretized into vectors so computations can use the forward and adjoint operators without explicitly constructing the system matrix.

2 Observation model

  • The model begins with an unknown image function observed through noisy samples of its Radon transform on fixed angle and detector grids.Source: §2 · p. 3
  • The measured sinogram is modeled as Radon-transform samples corrupted by independent, identically distributed Gaussian noise.The noise has mean zero and variance σ^2.Source: §2 · p. 3
  • The unknown image is discretized on a finite-dimensional basis, such as pixel values on a Cartesian grid, and flattened into a vector.The discrete Radon transform is a weighted sum of pixel values, with weights encoding ray-pixel intersection lengths.Source: §2 · p. 3
  • The flattened observations and image define a matrix A, but the simulations apply MATLAB's radon and iradon operators matrix-free.The noise variance is assumed given and may be estimated from observed data using a variance-based estimator.Source: §2 · pp. 3–4

3.1 Landweber iteration in Discrete Setting

The section presents classical Landweber iteration as an iterative method for linear inverse problems. For noisy data, its iterates require an appropriate stopping rule because the stopping index acts as a regularization parameter.

3.1 Landweber iteration in Discrete Setting

  • For noisy data, Landweber iteration exhibits semi-convergence, so the iteration must be stopped appropriately.With a suitable relaxation parameter, iterates for exact data converge to the generalized solution; appropriately chosen stopping rules can yield regularization methods under suitable smoothness assumptions [18].Source: §3.1 · p. 4
  • The iteration is rewritten using G = I − ωA^T A and f = ωA^T y; the relaxation parameter satisfies 0 < ω < 2/∥A^T A∥.The stated choice makes G a contraction, with spectral radius below one when A has full column rank.Source: §3.1 · p. 4
  • The paper uses the classical Landweber iteration to develop generalized fiducial inference for iterative methods.It defines T_K = ω Σ from s = 0 to K−1 of G^s A^T, then starts the iteration at x_0 = 0.Source: §3.1 · p. 4

3.2 Fiducial inference for the inverse problem

Exact fiducial inversion is unstable for ill-posed problems, so the paper substitutes the finite-step Landweber map and propagates auxiliary randomness through it. The resulting distribution has the classical iterate as its mean and quantifies propagated noise through its covariance, while early stopping introduces bias.

3.2 Fiducial inference for the inverse problem

  • In the linear data-generating model, the auxiliary variable U is Gaussian and the unknown parameter is fixed; fiducial inference inverts the model while retaining U's randomness.The stated inversion map assumes A has full column rank.Source: §3.2 · p. 5
  • Because exact inversion is unstable for ill-posed problems, the paper replaces it with a regularized inverse generated by finitely many Landweber iterations.The exact inversion's mean and covariance involve (A^T A)^−1 for an ill-conditioned matrix.Source: §3.2 · p. 5
  • At iteration K, the regularized fiducial map pushes the auxiliary randomness forward to produce a data-dependent distribution on x.This distribution does not depend on an initial belief distribution, unlike the cited probabilistic iterative method formulation.Source: §3.2 · p. 5
  • The distribution's mean is the classical K-th Landweber iterate, and its covariance Γ_K quantifies propagated observational noise.The iteration index controls both reconstruction regularization and the amount of uncertainty transferred from data to parameter space.Source: §3.2 · p. 5
  • Early stopping introduces bias in the regularized distribution, so uncertainty quantification requires estimating that bias.Source: §3.2 · p. 5

3.3 Distribution of the error in the reconstruction

The reconstruction error separates into deterministic regularization bias and stochastic error from observational noise. Its stochastic distribution is characterized by the covariance matrix ΓK.

3.3 Distribution of the error in the reconstruction

  • The mean of the regularized fiducial distribution is the usual K-th Landweber iterate, xK = TKy.Confidence-region construction therefore requires a quantitative bound on the bias alongside the stochastic error.Source: §3.3 · pp. 6–7
  • The error x† − xK follows a normal distribution with mean given by the bias term and covariance ΓK.The distribution is derived for Gaussian observational noise with x0 = 0.Source: §3.3 · p. 6
  • The reconstruction error has a deterministic bias component and a stochastic component induced by observational noise.The stochastic component is characterized by ΓK, while the bias results from terminating Landweber iteration after finitely many steps.Source: §3.3 · p. 7

3.4 Source condition and regularization bias

A source condition restricts admissible solutions by their regularity relative to the forward operator. Under this condition, the paper bounds the finite-iteration Landweber bias, with its decay rate determined by ν.

3.4 Source condition and regularization bias

  • The parameter ν quantifies solution smoothness relative to A^T A and determines the rate at which regularization error can decay.The condition is used to obtain a quantitative bound on the deterministic bias from terminating Landweber iteration.Source: §3.4 · p. 7
  • The source condition represents x† as (A^T A)^νw for a vector w with bounded norm.It restricts the class of admissible exact solutions by prescribing their regularity relative to the forward operator.Source: §3.4 · p. 7
  • Under the source condition, the bias has an upper bound that decays algebraically with K at a rate determined by ν.The stated estimate is bK(ν,A) ≤ Cν M K^-ν, with Cν depending on ν and ∥A∥.Source: §3.4 · p. 7

3.5 Construction of a guaranteed confidence region

The section constructs confidence regions by combining propagated stochastic uncertainty with a bound on regularization bias. It also develops pixelwise intervals and motivates a randomized low-rank covariance approximation for practical computation.

3.5 Construction of a guaranteed confidence region

  • Under the stated assumptions, CK(Y) is a (1 − α)-level confidence region for x†.The construction accounts for the deterministic bias bound as well as the stochastic reconstruction error.Source: §3.5 · p. 8
  • For imaging applications, the method also considers marginal, pixelwise confidence regions based on each coordinate of the reconstruction error.These complement the joint ellipsoidal confidence region with separate uncertainty estimates for individual pixels.Source: §3.5 · p. 9
  • The stopping index is chosen by balancing the regularization bias against the variance term capturing stochastic uncertainty.The stopping index depends on the noise level and forward problem, but not on the observed-data realization.Source: §3.5 · p. 9
  • A Monte Carlo method estimates the covariance diagonal entries needed for stopping and pixelwise intervals; the estimates are unbiased.The construction uses independent Rademacher vectors.Source: §3.5 · p. 9
  • Exact propagated covariance is impractical at large scale, and small-eigenvalue directions can inflate interval widths.The proposed response is a randomized rank-r approximation for the principal covariance contribution, with the omitted spectral component treated separately.Source: §3.5 · p. 10

3.6 Randomized SVD for practical UQ

The method uses randomized low-rank approximations of the normal operator to make uncertainty quantification more practical while retaining the standard Landweber reconstruction. It balances approximations to regularization bias and propagated uncertainty, then uses them to construct pixelwise confidence intervals.

3.6 Randomized SVD for practical UQ

  • The standard Landweber iteration still supplies the point estimate; randomized approximation of the normal operator is used only to evaluate propagated uncertainty.This avoids confining the reconstruction to the retained rank-r subspace, which could exclude components of the truth.Source: §3.6 · p. 10

3.6.1 Bias and Covariance approximation for TSVD

  • For large-scale problems, deterministic truncated SVD is computationally infeasible, motivating a randomized rank-r approximation for the stopping rule and uncertainty estimates.The truncated-SVD bias analysis notes that omitted-tail bias can be smaller when the truth has little projection onto omitted eigendirections, motivating a practical rule that ignores this tail contribution.Source: §3.6.1 · p. 12

3.6.2 Stopping rule and confidence regions from rSVD

  • The randomized stopping rule balances an approximate bias against pixelwise uncertainty from covariance diagonals, estimating the residual covariance diagonal rather than the full covariance diagonal.The residual diagonal estimate can be more accurate when the low-rank approximation captures the dominant covariance structure well.Source: §3.6.2 · p. 13
  • Under the stated assumptions, the pixelwise confidence intervals have a coverage guarantee whose lower bound can approach the desired level, while simulations indicate coverage will usually exceed it.The guarantee assumes the randomized SVD and residual-diagonal probes are independent of the observational noise.Source: §3.6.2 · pp. 13–14

4 Numerical Simulation Results

The simulations evaluate randomized-rank uncertainty estimates for a noisy Shepp-Logan CT reconstruction. Empirical coverage reflects a rank–noise trade-off: low ranks can under-cover at low noise, while larger ranks can yield conservative coverage at higher noise.

4.1 Simulation set-up

  • The simulations use a 128 × 128 Shepp-Logan phantom, Radon data at 180 angles, Gaussian noise, and matrix-free Landweber operations.They assess pixelwise intervals and empirical frequentist coverage over repeated noisy-data trials.Source: §4.1 · pp. 15–16, §4 · p. 15

4.2 Simulation results and interpretation

  • At fixed rank, empirical coverage increases with noise; once dominant spectral components are captured, adding rank can inflate intervals and make coverage increasingly conservative.At rank 700, mean interior coverage rises from 0.8968 for σ = 0.1 to 0.9432 for σ = 0.2 and 0.9599 for σ = 0.3.Source: §4.2 · p. 16
  • Higher noise leads the stopping rule to stop earlier, whereas increasing retained rank allows more Landweber iterations, consistent with complementary regularization from stopping and spectral truncation.The study concludes that rank should depend on noise level and identifies adaptive rank selection as future work.Source: §4.2 · pp. 16, 19
  • At σ = 0.1, mean interior coverage rises from 0.8740 at rank 500 to 0.9311 at rank 1300, indicating low-rank undercoverage relative to the nominal 0.95.Table 1 reports coverage for different randomized SVD ranks and noise levels; interior coverage is computed inside the phantom.Source: §4.2 · p. 16, Table 1 · p. 19

5 Appendix

The appendix analyzes how randomized SVD approximation of the normal operator affects the propagated covariance. It provides a probabilistic covariance-error bound and discusses rank selection and computational alternatives.

5.1 Covariance matrix from SVD

  • For p ≥ 4, the randomized-SVD covariance approximation satisfies a probabilistic error bound with probability at least 1 − 3e^-p.The error threshold is denoted Δ_K,r,p.Source: §5.1 · pp. 20–21
  • The randomized approximation uses a range finder, a reduced eigen-decomposition, and truncation to obtain a rank-r approximation of the positive semidefinite normal operator.The construction uses an oversampling parameter p alongside target rank r.Source: §5.1 · pp. 19, 21
  • Rank can be chosen adaptively to meet a covariance-approximation tolerance with high probability, at extra one-time computational cost.The appendix also notes that the relevant bound is expected to decay as rank increases, independently of discretization.Source: §5.1 · p. 21
Loading 2610.00882v1…