Fast linear inversion of highly overdetermined linearized inverse scattering problems, introduced by Markel, Levinson and Schotland (Inverse Problems 35, 124002, 2019), exploits the separable structure of the forward matrix to solve diffuse optical tomography problems with up to about 10810^8 data points without ever assembling the full matrix. Its central data-reduction step, however, truncates a diagonal-scaling preconditioner at a cutoff ε\varepsilon chosen by inspecting singular spectra, with no systematic criterion. The open problem is to build the missing theory: quantitative or adaptive rules for the truncation cutoffs that balance information retention, the intrinsic ill-posedness of the inverse problem, and finite-precision round-off; rigorous error-propagation estimates under realistic noise models beyond additive shot noise, including the paper's modified unscaled diagonal operator; and the analytic mapping between these geometry-independent algebraic methods and the symmetry-exploiting analytic methods for slab geometries, which the authors deferred.

Background

In diffuse optical tomography (DOT) one reconstructs the spatially varying absorption coefficient of a highly light-scattering medium, such as biological tissue, from measurements of diffuse light collected at the boundary under many source illuminations; light propagation is modeled by the diffusion approximation to the radiative transport equation, and the reconstruction task is a nonlinear inverse problem. Under the Born (linearizing) approximation the unknown perturbation VV of the absorption is assumed small, so the data depend linearly on VV. After discretization, the linearized forward map is a dense matrix KK whose rows are indexed by source-detector pairs (i,j)(i,j) — M=NdNsM=N_dN_s rows, with NdN_d detectors and NsN_s sources — and whose columns are indexed by the N=NvN=N_v voxels. When the forward model uses Green's functions of the background medium, this matrix has the separable (Khatri-Rao) structure K(ij),n=AinBnj,K_{(ij),n}=A_{in}B_{nj}, where AA contains detector-to-voxel and BB voxel-to-source Green's-function values. In modern noncontact experiments MM can reach 10710^7-10810^8 while N∼104N\sim10^4-10510^5: the system is highly overdetermined, and generic dense linear algebra — a singular value decomposition (SVD, the factorization of a matrix as unitary times diagonal times unitary) at O(MN2)O(MN^2) cost and O(MN)O(MN) storage — is infeasible.

Markel, Levinson and Schotland (2019) observed that the separable structure allows the normal matrix and right-hand side to be computed without forming KK: with Φ\Phi the Nd×NsN_d\times N_s data matrix, (K∗K)nm=(A∗A)nm(BB∗)mn,cn=(A∗ΦB∗)nn,(K^*K)_{nm}=(A^*A)_{nm}(BB^*)_{mn},\qquad c_n=(A^*\Phi B^*)_{nn}, where ∗^* denotes conjugate transpose and the first product is the Hadamard (entrywise) product. This yields the regularized normal equations (K∗K+λ2I)x=c(K^*K+\lambda^2I)x=c (Tikhonov regularization with parameter λ\lambda) at a cost independent of the large factor MM. They further proposed a diagonal-scaling preconditioner: from the SVDs A=∑μσμAfμAgμA∗A=\sum_\mu\sigma_\mu^A f_\mu^Ag_\mu^{A*} and B=∑νσνBfνBgνB∗B=\sum_\nu\sigma_\nu^B f_\nu^Bg_\nu^{B*}, transform the system by the unitary map built from the singular vectors and rescale each retained mode-pair equation (μ,ν)(\mu,\nu) by 1/(σμAσνB)1/(\sigma_\mu^A\sigma_\nu^B), discarding equations whose singular values fall below a cutoff ε\varepsilon. This reduces the effective number of equations from NdNsN_dN_s to MAMBM_AM_B (retained mode pairs) and enabled direct inversions of DOT problems with up to M=101,727,396M=101{,}727{,}396 data points and Nv=68,921N_v=68{,}921 unknowns on a supercomputer (hours of wall-clock time, about 76 GB at the 2Nv22N_v^2-word memory bottleneck), with roughly six-fold memory reduction from the cutoff in the largest example. The efficiency matters because nonlinear reconstruction schemes, such as Newton-type iterations or data-compatible T-matrix completion — recovering the scattering transition matrix, which relates incident fields to scattered fields, from incomplete data (Levinson and Markel 2016) — require a linearized inversion at each outer step. An immediate predecessor line (Machida, Panasyuk, Wang, Markel and Schotland 2016) treated large DOT datasets in the radiative transport regime with mode-expansion techniques.

The method's weak point is that the cutoff ε\varepsilon — defined absolutely on the singular values in eq. (13) of the paper — is chosen by inspecting the singular spectra: for a 21321^3-voxel cubic sample, ε=10−3max⁡μσμA\varepsilon=10^{-3}\max_\mu\sigma_\mu^A was placed inside an observed gap in the singular values of AA, while for a 41341^3 cube, ε=10−2max⁡μσμA\varepsilon=10^{-2}\max_\mu\sigma_\mu^A was, per the paper's Section 5.4, well above that gap (reducing it to ∼10−4\sim10^{-4}, inside the gap, gave no improvement of image quality) — in neither case with any criterion linking the cutoff to noise level, conditioning, or reconstruction fidelity. In the additive Poisson shot-noise simulations, the bare cutoff value ε=0.1\varepsilon=0.1 reported in Section 6 of the paper retained only 28×28=78428\times28=784 equations on the 21321^3 cube against 9,2619{,}261 unknowns and preserved the outer shell of a two-shell target while losing the inner one. The 1/σ1/\sigma rescaling amplifies noise in weak modes, and a modified diagonal operator (eq. (23) of the paper) that discards weak equations without rescaling was proposed to moderate this but never analyzed quantitatively. The simulations used additive shot noise only. For measurement geometries with translational symmetry, such as the slab, the linearized forward operator can be diagonalized analytically: Markel, Mital and Schotland (JOSA A 20, 890, 2003) derived closed-form inversion formulas and the singular-value decomposition of the slab-geometry diffusion tomography operator via Fourier transformation in the invariant directions. These symmetry-exploiting analytic methods motivated the algebraic, geometry-independent construction of 2019; in its Section 7 the source notes that the fast inversion methods of its refs. [28,29] (Markel and Schotland 2002, 2004), which require translational symmetry, are related to the present method when the SVD decompositions of AA and BB can be computed analytically, and explicitly defers the mathematical connection between the two approaches ("will be explained elsewhere"). A related rigorous framework appeared later: Dölz, Egger and Schlottbom (2021) cast inverse problems with the same factorized structure T(c)=V0D(c)UT(c)=V_0D(c)U as operator-valued data problems in Hilbert spaces and constructed certified reduced-order models of quasi-optimal rank with complete error analysis under additive Hilbert-Schmidt-norm noise; stochastic row-sampling algorithms address generic overdetermined inverse problems with many measurements (Roosta-Khorasani, van den Doel and Ascher 2014). None of these analyzes the diagonal-scaling preconditioner of the 2019 method itself.

Problem Statement

The problem is to develop the systematic theory that converts the diagonal-scaling data-reduction method of Markel, Levinson and Schotland (2019) for highly overdetermined linearized inverse scattering problems from an empirically tuned algorithm into a principled method with predictable accuracy. The common objective uniting the requested results is a complete quantitative understanding of what the truncation step preserves, amplifies, and destroys; specifically:

(1) Truncation criteria. Derive systematic criteria or adaptive rules for the diagonal-scaling cutoff ε\varepsilon — and for the cutoff of the modified diagonal operator (eq. (23) of the source) that discards equations without the 1/σ1/\sigma rescaling — that account simultaneously for information retention (which singular mode pairs of the factor matrices AA and BB carry recoverable signal), the intrinsic ill-posedness of the inverse problem (stability and regularization behavior of the reduced normal operator W=(A+A)∘(BB+)W=(A^+A)\circ(BB^+), with ∘\circ the Hadamard product and A+A^+, B+B^+ the truncated pseudoinverses), and finite-precision round-off.

(2) Noise propagation. Provide rigorous error-propagation estimates for the truncated and preconditioned reconstructions under noisy data, extending beyond the additive Poisson shot-noise model used in the source's simulations to at least one structurally different noise model (e.g., multiplicative noise), and quantify the noise-amplification trade-off between the rescaling diagonal operator and the unscaled modified operator of eq. (23).

(3) Analytic connection. Establish the mathematical connection, deferred by the authors, between the geometry-independent algebraic construction (SVD-based diagonal scaling of the factor matrices) and the symmetry-exploiting analytic methods for slab/translational-invariance geometries (the fast inversion methods of Markel and Schotland 2002, 2004, cited as refs. [28,29] of the source), in the regime, identified by the source, where the SVD decompositions of AA and BB can be computed analytically.

These are independently identified open directions of the source rather than a jointly mandated program; a numerical implementation (for instance, integration of the eq.-(23) variant into a full reconstruction pipeline) may accompany the theory, but programming-level memory-reduction strategies alone are not part of the scientific completion boundary.

References

  1. Markel, V. A., Levinson, H., and Schotland, J. C. (2019). Fast linear inversion for highly overdetermined inverse scattering problems. Inverse Problems 35, 124002. DOI: 10.1088/1361-6420/ab44e7. https://doi.org/10.1088/1361-6420/ab44e7
  2. Markel, V. A., Mital, V., and Schotland, J. C. (2003). Inverse problem in optical diffusion tomography III: inversion formulas and singular-value decomposition. Journal of the Optical Society of America A 20, 890. DOI: 10.1364/josaa.20.000890. https://doi.org/10.1364/josaa.20.000890
  3. Machida, M., Panasyuk, G. Y., Wang, Z.-M., Markel, V. A., and Schotland, J. C. (2016). Radiative transport and optical tomography with large datasets. Journal of the Optical Society of America A 33, 551. DOI: 10.1364/josaa.33.000551. https://doi.org/10.1364/josaa.33.000551
  4. Dölz, J., Egger, H., and Schlottbom, M. (2021). A model reduction approach for inverse problems with operator valued data. Numerische Mathematik 148, 889-917. DOI: 10.1007/s00211-021-01224-5. Preprint arXiv:2004.11827. https://arxiv.org/abs/2004.11827
  5. Markel, V. A. (2020). Investigation of the effect of super-resolution in nonlinear inverse scattering. Physical Review E 102, 053313. DOI: 10.1103/PhysRevE.102.053313. https://doi.org/10.1103/PhysRevE.102.053313
  6. Wang, Z., Hauptmann, A., Lu, L., and Schotland, J. C. (2025). Learning the detector in optical tomography. arXiv:2510.05341. https://arxiv.org/abs/2510.05341
  7. Khan, A., and Venkatapathi, M. (2021). Approximations for modeling light scattering by spheres with uncertainty in physical parameters. arXiv:2112.13786. https://arxiv.org/abs/2112.13786
  8. Roosta-Khorasani, F., van den Doel, K., and Ascher, U. (2014). Stochastic algorithms for inverse problems involving PDEs and many measurements. SIAM Journal on Scientific Computing 36, S3-S22. DOI: 10.1137/130922756. https://doi.org/10.1137/130922756
  9. Levinson, H. W., and Markel, V. A. (2016). Solution of the nonlinear inverse scattering problem by T-matrix completion. I. Theory. Physical Review E 94, 043317. DOI: 10.1103/PhysRevE.94.043317. https://doi.org/10.1103/PhysRevE.94.043317
  10. Markel, V. A., and Schotland, J. C. (2002). Effects of sampling and limited data in optical tomography. Applied Physics Letters 81, 1180-1182. DOI: 10.1063/1.1495543. https://doi.org/10.1063/1.1495543 (Reference [28] of the source.)
  11. Markel, V. A., and Schotland, J. C. (2004). Symmetries, inversion formulas, and image reconstruction for optical tomography. Physical Review E 70, 056616. DOI: 10.1103/PhysRevE.70.056616. https://doi.org/10.1103/PhysRevE.70.056616 (Reference [29] of the source.)