Inverse Problems

In this work, the inverse problem of quantitative thermoacoustic tomography is studied. In quantitative thermoacoustic tomography, dielectric parameters of an imaged target are estimated from an absorbed energy density induced by an externally introduced electromagnetic excitation. In this work, simultaneous estimation of electrical conductivity and permittivity is considered. We approach this problem in the framework of Bayesian inverse problems. The dielectric parameters are estimated by computing maximum a posteriori estimates, and the reliability of the estimates is studied using the Laplace's approximation. The forward model to describe electromagnetic wave propagation is based on a vectorial Maxwell's equation, that is numerically approximated using a finite element method with edge elements. The proposed methodology is evaluated using numerical simulations utilising one or two electromagnetic excitations at multiple excitation frequencies. The results show that the electrical conductivity and permittivity can be simultaneously estimated in quantitative thermoacoustic tomography. However, the problem can suffer from non-uniqueness, which could be overcome using multiple electromagnetic excitations.

Compressive sensing magnetic resonance imaging (CS-MRI) accelerates data acquisition by reconstructing high-quality images from a limited set of k-space samples. To solve this ill-posed inverse problem, the plug-and-play (PnP) framework integrates image priors using convolutional neural network (CNN) denoisers. However, CNN denoisers often prioritize local details and may neglect broader degradation effects, leading to visually plausible but structurally inaccurate artifacts. Additionally, the theoretical convergence of PnP methods remains a significant challenge. In this work, we propose a novel method, Plug-And-pLAy 3D MRI recoNstruction, to bridge the gap between denoising and MRI reconstruction. Our model employs the tensor tubal nuclear norm (TNN) to capture intrinsic correlations in 3D MRI data. It also incorporates two implicit regularizers. The first leverages CNN denoisers to exploit image priors. The second, introduced here for the first time, is formulated as a CS-MRI reconstruction subproblem and solved using a deep learning-based method to preserve global spatial structure. We solve the proposed model using the alternating direction method of multipliers. We extend existing theoretical results to prove the algorithm's convergence to a fixed point under reasonable assumptions. Experiments on two datasets with three sampling masks show that our method outperforms state-of-the-art MRI reconstruction methods. Ablation studies confirm that the TNN and the two implicit regularizers work together to improve reconstruction quality.

This paper concerns the inverse shape problem of recovering an unknown clamped cavity embedded in a thin infinite plate. The model problem is assumed to be governed by the two-dimensional biharmonic wave equation in the frequency domain. Based on the far-field data, a resolution analysis is conducted for cavity recovery via the direct sampling method. The Funk–Hecke integral identity is employed to analyze the performance of two imaging functions. Our analysis demonstrates that the same imaging functions commonly used for acoustic inverse shape problems are applicable to the biharmonic wave context. This work presents the first extension of direct sampling methods to biharmonic waves using far-field data. Numerical examples are provided to illustrate the effectiveness of these imaging functions in recovering a clamped cavity.

Consider the linear ill-posed problems of the form $\sum_{i = 1}^{b} A_i x_i = y$, where, for each i, A_i is a bounded linear operator between two Hilbert spaces X_i and ${\mathcal Y}$. When b is huge, solving the problem by an iterative method using the full gradient at each iteration step is both time-consuming and memory insufficient. Although randomized block coordinate descent (RBCD) method has been shown to be an efficient method for well-posed large-scale optimization problems with a small amount of memory, there still lacks a convergence analysis on the RBCD method for solving ill-posed problems. In this paper, we investigate the convergence property of the RBCD method with noisy data under either a priori or a posteriori stopping rules. We prove that the RBCD method combined with an a priori stopping rule yields a sequence that converges weakly to a solution of the problem almost surely. We also consider the early stopping of the RBCD method and demonstrate that the discrepancy principle can terminate the iteration after finite many steps almost surely. For a class of ill-posed problems with special tensor product form, we obtain strong convergence results on the RBCD method. Furthermore, we consider incorporating the convex regularization terms into the RBCD method to enhance the detection of solution features. To illustrate the theory and the performance of the method, numerical simulations from the imaging modalities in computed tomography and compressive temporal imaging are reported.

We propose and analyze a stochastic-gradient type method for solving systems of nonlinear ill-posed equations. The method considered here extends the stochastic gradient descent (SGD) type iteration introduced in Rabelo et al (2022 Inverse Problems38 025003) for solving linear ill-posed systems. A distinctive feature of our method resides in the adaptive choice of the stepsize, which promotes a relaxed orthogonal projection of the current iterate onto a conveniently chosen convex set. This characteristic distinguish our method from other SGD type methods in the literature (where the stepsize is typically chosen a priori) and accounts for the faster convergence observed in the numerical experiments conducted in this manuscript. The convergence analysis discussed here includes: monotonicity and mean square convergence of the iteration error (exact data case), stability and semi-convergence (noisy data case). In the later case, our method is coupled with an a priori stopping rule. Numerical experiments are presented for two large scale nonlinear inverse problems in machine learning (both with real data): (i) we address, using neural networks, the big data problem of CO-concentration prediction considered in the above cited article; (ii) we tackle the classification problem for the MNIST database (http://yann.lecun.com/exdb/mnist/). Additionally, a parameter identification problem in a 3D elliptic PDE system is considered.

Recent years have witnessed a growth in mathematics for deep learning—which seeks a deeper understanding of the concepts of deep learning with mathematics and explores how to make it more robust—and deep learning for mathematics, where deep learning algorithms are used to solve problems in mathematics. The latter has popularised the field of scientific machine learning where deep learning is applied to problems in scientific computing. Specifically, more and more neural network (NN) architectures have been developed to solve specific classes of partial differential equations (PDEs). Such methods exploit properties that are inherent to PDEs and thus solve the PDEs better than standard feed-forward NNs, recurrent NNs, or convolutional neural networks. This has had a great impact in the area of mathematical modelling where parametric PDEs are widely used to model most natural and physical processes arising in science and engineering. In this work, we review such methods as well as their extensions for parametric studies and for solving the related inverse problems. We also show their relevance in various industrial applications.

We present a review of methods for optimal experimental design (OED) for Bayesian inverse problems governed by partial differential equations with infinite-dimensional parameters. The focus is on problems where one seeks to optimize the placement of measurement points, at which data are collected, such that the uncertainty in the estimated parameters is minimized. We present the mathematical foundations of OED in this context and survey the computational methods for the class of OED problems under study. We also outline some directions for future research in this area.

Over the last decades, the total variation (TV) has evolved to be one of the most broadly-used regularisation functionals for inverse problems, in particular for imaging applications. When first introduced as a regulariser, higher-order generalisations of TV were soon proposed and studied with increasing interest, which led to a variety of different approaches being available today. We review several of these approaches, discussing aspects ranging from functional-analytic foundations to regularisation theory for linear inverse problems in Banach space, and provide a unified framework concerning well-posedness and convergence for vanishing noise level for respective Tikhonov regularisation. This includes general higher orders of TV, additive and infimal-convolution multi-order total variation, total generalised variation, and beyond. Further, numerical optimisation algorithms are developed and discussed that are suitable for solving the Tikhonov minimisation problem for all presented models. Focus is laid in particular on covering the whole pipeline starting at the discretisation of the problem and ending at concrete, implementable iterative procedures. A major part of this review is finally concerned with presenting examples and applications where higher-order TV approaches turned out to be beneficial. These applications range from classical inverse problems in imaging such as denoising, deconvolution, compressed sensing, optical-flow estimation and decompression, to image reconstruction in medical imaging and beyond, including magnetic resonance imaging, computed tomography, magnetic-resonance positron emission tomography, and electron tomography.

Electrical impedance tomography (EIT) is an imaging modality where a patient or object is probed using harmless electric currents. The currents are fed through electrodes placed on the surface of the target, and the data consists of voltages measured at the electrodes resulting from a linearly independent set of current injection patterns. EIT aims to recover the internal distribution of electrical conductivity inside the target. The inverse problem underlying the EIT image formation task is nonlinear and severely ill-posed, and hence sensitive to modeling errors and measurement noise. Therefore, the inversion process needs to be regularized. However, traditional variational regularization methods, based on optimization, often suffer from local minima because of nonlinearity. This is what makes regularized direct (non-iterative) methods attractive for EIT. The most developed direct EIT algorithm is the D-bar method, based on complex geometric optics solutions and a nonlinear Fourier transform. Variants and recent developments of D-bar methods are reviewed, and their practical numerical implementation is explained.

Magnetic particle imaging (MPI) is a relatively new imaging modality. The nonlinear magnetization behavior of nanoparticles in an applied magnetic field is exploited to reconstruct an image of the concentration of nanoparticles. Finding a sufficiently accurate model to reflect the behavior of large numbers of particles for MPI remains an open problem. As such, reconstruction is still computed using a measured forward operator obtained in a time-consuming calibration process. The model commonly used to illustrate the imaging methodology and obtain first model-based reconstructions relies on substantial model simplifications. By neglecting particle–particle interactions, the forward operator can be expressed by a Fredholm integral operator of the first kind when describing the inverse problem. Here, we review previously proposed models derived from single-particle behavior in the MPI context and consider future research on linear and nonlinear problems beyond concentration reconstruction applications. This survey complements a recent topical review on MPI (Knopp et al 2017 Phys. Med. Biol. 62 R124).

The functional linear regression model has been widely studied and utilized for dealing with functional predictors. In this paper, we study the Nyström subsampling method, a strategy used to tackle the computational complexities inherent in big data analytics, especially within the domain of functional linear regression model in the framework of reproducing kernel Hilbert space. 

By adopting a Nyström subsampling strategy, our aim is to mitigate the computational overhead associated with kernel methods, which often struggle to scale gracefully with dataset size. Specifically, we investigate a regularization-based approach combined with Nyström subsampling for functional linear regression model, effectively reducing the computational complexity from $O(n^3)$ to $O(m^2 n)$, where $n$ represents the size of the observed empirical dataset and $m$ is the size of subsampled dataset. Notably, we establish that these methodologies will achieve optimal convergence rates, provided that the subsampling level is appropriately selected. We have also demonstrated a numerical example of Nyström subsampling in the RKHS framework for the functional linear regression model.

We introduce subgradient-based Lavrentiev regularisation of the form A(u) + α ∂R(u) ∋ f^δ for linear and nonlinear ill-posed problems with monotone operators A and general regularisation functionals R. In contrast to Tikhonov regularisation, this approach perturbs the equation itself and avoids the use of the adjoint of the derivative of A. It is therefore especially suitable for time-causal problems that only depend on information in the past and allows for real-time computation of regularised solutions. We establish a general well-posedness theory in Banach spaces and prove convergence-rate results with variational source conditions. Furthermore, we demonstrate its application in total-variation denoising in linear Volterra integral operators of the first kind and parameter-identification problems in semilinear parabolic PDEs.

In this paper, we investigate an inverse Cauchy problem for a stochastic hyperbolic equation. A Lipschitz type observability estimate is established using a pointwise Car- leman identity. By minimizing the constructed Tikhonov-type functional, we obtain a regularized approximation to the problem. The properties of the approximation are studied by means of the Carleman estimate and Riesz representation theorem. Leveraging kernel-based learning theory, we simulate numerical algorithms based on the proposed regularization method. These reconstruction algorithms are implemented and validated through several numerical experiments, demonstrating their feasibility and accuracy.

We study the problem of finding the resistors in a resistor network from measurements of the power dissipated by the resistors under different loads. We give sufficient conditions for local uniqueness, i.e. conditions that guarantee that the linearization of this non-linear inverse problem admits a unique solution. Our method is inspired by a method to study local uniqueness of inverse problems with internal functionals in the continuum, where the inverse problem is reformulated as a redundant system of differential equations. We use our method to derive local uniqueness conditions for other discrete inverse problems with internal functionals including a discrete analogue of the inverse Schrödinger problem and problems where the resistors are replaced by impedances and dissipated power at the zero and a positive frequency are available. Moreover, we show that the dissipated power measurements can be obtained from measurements of thermal noise induced currents.

This paper presents an accelerated preconditioned proximal gradient algorithm (APPGA) for effectively solving a class of positron emission tomography (PET) image reconstruction models with differentiable regularizers. We establish the convergence of APPGA with the generalized Nesterov (GN) momentum scheme, demonstrating its ability to converge to a minimizer of the objective function with rates of o(1/k^{2ω}) and o(1/k^ω) in terms of the function value and the distance between consecutive iterates, respectively, where ω∈(0, 1] is the power parameter of the GN momentum. To achieve an efficient algorithm with high-order convergence rate for the higher-order isotropic total variation (ITV) regularized PET image reconstruction model, we replace the ITV term by its smoothed version and subsequently apply APPGA to solve the smoothed model. Numerical results presented in this work indicate that as ω∈(0, 1] increase, APPGA converges at a progressively faster rate. Furthermore, APPGA exhibits superior performance compared to the preconditioned proximal gradient algorithm and the preconditioned Krasnoselskii-Mann algorithm. The extension of the GN momentum technique for solving a more complex optimization model with multiple nondifferentiable terms is also discussed.

We introduce subgradient-based Lavrentiev regularisation of the form A(u) + α ∂R(u) ∋ f^δ for linear and nonlinear ill-posed problems with monotone operators A and general regularisation functionals R. In contrast to Tikhonov regularisation, this approach perturbs the equation itself and avoids the use of the adjoint of the derivative of A. It is therefore especially suitable for time-causal problems that only depend on information in the past and allows for real-time computation of regularised solutions. We establish a general well-posedness theory in Banach spaces and prove convergence-rate results with variational source conditions. Furthermore, we demonstrate its application in total-variation denoising in linear Volterra integral operators of the first kind and parameter-identification problems in semilinear parabolic PDEs.

In this work, the inverse problem of quantitative thermoacoustic tomography is studied. In quantitative thermoacoustic tomography, dielectric parameters of an imaged target are estimated from an absorbed energy density induced by an externally introduced electromagnetic excitation. In this work, simultaneous estimation of electrical conductivity and permittivity is considered. We approach this problem in the framework of Bayesian inverse problems. The dielectric parameters are estimated by computing maximum a posteriori estimates, and the reliability of the estimates is studied using the Laplace's approximation. The forward model to describe electromagnetic wave propagation is based on a vectorial Maxwell's equation, that is numerically approximated using a finite element method with edge elements. The proposed methodology is evaluated using numerical simulations utilising one or two electromagnetic excitations at multiple excitation frequencies. The results show that the electrical conductivity and permittivity can be simultaneously estimated in quantitative thermoacoustic tomography. However, the problem can suffer from non-uniqueness, which could be overcome using multiple electromagnetic excitations.

This paper concerns the inverse shape problem of recovering an unknown clamped cavity embedded in a thin infinite plate. The model problem is assumed to be governed by the two-dimensional biharmonic wave equation in the frequency domain. Based on the far-field data, a resolution analysis is conducted for cavity recovery via the direct sampling method. The Funk–Hecke integral identity is employed to analyze the performance of two imaging functions. Our analysis demonstrates that the same imaging functions commonly used for acoustic inverse shape problems are applicable to the biharmonic wave context. This work presents the first extension of direct sampling methods to biharmonic waves using far-field data. Numerical examples are provided to illustrate the effectiveness of these imaging functions in recovering a clamped cavity.

Due to their uncertainty quantification, Bayesian solutions to inverse problems are the framework of choice in applications that are risk averse. These benefits come at the cost of computations that are in general, intractable. New advances in machine learning and variational inference (VI) have lowered this computational barrier by leveraging data-driven learning. Two VI paradigms have emerged that represent different tradeoffs: amortized and non-amortized. Amortized VI can produce fast results but due to generalizing to many observed datasets it produces suboptimal inference results. Non-amortized VI is slower at inference but finds better posterior approximations since it is specialized towards a single observed dataset. Current amortized VI techniques run into a sub-optimality wall that can not be improved without more expressive neural networks or extra training data. We present a solution that enables iterative improvement of amortized posteriors that uses the same networks architectures and training data. The benefits of our method requires extra computations but these remain frugal since they are based on physics-hybrid methods and summary statistics. Importantly, these computations remain mostly offline thus our method maintains cheap and reusable online evaluation while bridging the optimality gap between these two paradigms. We denote our proposed method ASPIRE - Amortized posteriors with Summaries that are Physics-based and Iteratively REfined. We first validate our method on a stylized problem with a known posterior then demonstrate its practical use on a high-dimensional and nonlinear transcranial medical imaging problem with ultrasound. Compared with the baseline and previous methods in the literature, ASPIRE stands out as an computationally efficient and high-fidelity method for posterior inference.

Iterative regularization is a classic idea in regularization theory, that has recently become popular in machine learning. On the one hand, it allows to design efficient algorithms controlling at the same time numerical and statistical accuracy. On the other hand it allows to shed light on the learning curves observed while training neural networks. In this paper, we focus on iterative regularization in the context of classification. After contrasting this setting with that of linear inverse problems, we develop an iterative regularization approach based on the use of the hinge loss function. More precisely we consider a diagonal approach for a family of algorithms for which we prove convergence as well as rates of convergence and stability results for a suitable classification noise model. Our approach compares favorably with other alternatives, as confirmed by numerical simulations.

This paper studies inverse problems in quantitative photoacoustic tomography with additional optical current data supplemented from diffuse optical tomography. We propose a three-stage image reconstruction method for the simultaneous recovery of the absorption, diffusion, and Grüneisen coefficients. We demonstrate, through numerical simulations, that: (i) when the Grüneisen coefficient is known, the addition of the optical measurements allows a more accurate reconstruction of the scattering and absorption coefficients; and (ii) when the Grüneisen coefficient is not known, the addition of optical current measurements allows us to reconstruct uniquely the Grüneisen, the scattering and absorption coefficients. Numerical simulations based on synthetic data are presented to demonstrate the effectiveness of the proposed idea.

The forward problem arising in several hybrid imaging modalities can be modeled by the Cauchy problem for the free space wave equation. Solution to this problems describes propagation of a pressure wave, generated by a source supported inside unit sphere S. The data g represent the time-dependent values of the pressure on the observation surface S. Finding initial pressure f from the known values of g consitutes the inverse problem. The latter is also frequently formulated in terms of the spherical means of f with centers on S. Here we consider a problem of range description of the wave operator mapping f into g. Such a problem was considered before, with data g known on time interval at least $[0,2]$ (assuming the unit speed of sound). Range conditions were also found in terms of spherical means, with radii of integration spheres lying in the range $[0,2]$. However, such data are redundant. We present necessary and sufficient conditions for function g to be in the range of the wave operator, for g given on a half-time interval $[0,1]$. This also implies range conditions on spherical means measured for the radii in the range $[0,1]$.

Poisson distributed measurements in inverse problems often stem from Poisson point processes that are observed through discretized or finite-resolution detectors, one of the most prominent examples being positron emission tomography (PET). These inverse problems are typically reconstructed via Bayesian methods. A natural question then is whether and how the reconstruction converges as the signal-to-noise ratio tends to infinity and how this convergence interacts with other parameters such as the detector size. In this article we carry out a corresponding variational analysis for the exemplary Bayesian reconstruction functional from Mauritz et al (2024 SIAM J. Math. Anal.56 5840–80); Schmitzer et al (2020 IEEE Trans. Med. Imaging39 1626–35), which considers dynamic PET imaging (i.e. the object to be reconstructed changes over time) and uses an optimal transport regularization.

The inverse radiation therapy planning problem is crucial in achieving tumoricidal doses while minimizing radiation-induced normal-tissue toxicity. Compared to conventional radiation therapy (with conventional dose rates), flash proton radiation therapy (with ultra-high dose rates) can provide additional normal tissue sparing. However, this groundbreaking advancement in radiation oncology introduces a challenging nonconvex and nonsmooth optimization problem. In this paper, we propose a stochastic three-operator splitting (STOS) algorithm to address the flash proton radiation therapy problem. We establish the convergence and convergence rates of the STOS algorithm under the nonconvex framework for both unbiased gradient estimators and variance-reduced gradient estimators. These stochastic gradient estimators include the most popular ones, such as SGD, SAGA, SARAH, and SAG, among others. Experimental results demonstrate that the flash proton radiation therapy plans obtained by the STOS algorithm can effectively kill tumors while better protecting normal tissues.

We introduce parametric level-set enhanced to improve reconstruction (PaLEnTIR), a significantly enhanced parametric level-set (PaLS) method addressing the restoration and reconstruction of piecewise constant objects. Our key contribution involves a unique PaLS formulation utilizing a single level-set function to restore scenes containing multi-contrast piecewise constant objects without requiring knowledge of the number of objects or their contrasts. Unlike standard PaLS methods employing radial basis functions (RBFs), our model integrates anisotropic basis functions (ABFs), thereby expanding its capacity to represent a wider class of shapes. Furthermore, PaLEnTIR streamlines the model by reducing redundancy and indeterminacy in the parameterization, resulting in improved numerical performance. We compare PaLEnTIR's performance to state-of-the art alternatives via a diverse collection of experiments encompassing denoising, deconvolution, sparse and limited angle of view x-ray computed tomography (2D and 3D), and nonlinear diffuse optical tomography tasks using both real and simulated data sets.