Inverse Problems

Partial differential equation (PDE)-constrained inverse problems are central to many scientific research areas and engineering applications, including medical imaging, seismic imaging, and electromagnetic exploration. Full waveform inversion (FWI) is a typical PDE-constrained optimization problem widely used in seismic exploration and has achieved notable real-world successes. However, FWI is susceptible to local minima, making it difficult to obtain accurate inversion results without a sufficiently good initial model. Wavefield reconstruction inversion (WRI) has emerged as a promising approach to address the local minima issues inherent in FWI. Despite this potential, WRI suffers from the high computational cost associated with solving augmented systems for optimal data-fitting wavefields, limiting its feasibility in realistic 3D scenarios. This study presents a GPU-accelerated 3D source free adaptive WRI (GPU-SF-AWRI) method that overcomes these computational limitations by adaptively controlling the accuracy of wavefield simulations and optimizing GPU utilization, thereby making WRI more applicable to 3D problems. The inclusion of an on-the-fly source estimation technique further boosts its performance on realistic problems. Numerical experiments demonstrate that the proposed GPU-accelerated method achieves a 195-fold speedup over CPU-based implementations. By incorporating adaptive accuracy adjustments and total variation regularization, the method achieves a twofold speedup while preserving inversion accuracy. We applied the GPU-AWRI method to both synthetic and real physical modeling data, demonstrating its effectiveness in handling real data challenges and mitigating the local minima issues associated with conventional FWI.

Online optimisation studies the convergence of optimisation methods as the data embedded in the problem changes. Based on this idea, we propose a primal dual online method for nonlinear time-discrete inverse problems. We analyse the method through regret theory and demonstrate its performance in real-time monitoring of moving bodies in a fluid with electrical impedance tomography. To do so, we also prove the second-order differentiability of the complete electrode model solution operator on $L^\infty$.

In this paper we study the simultaneous reconstruction of two coefficients in a reaction–subdiffusion equation, namely a nonlinearity and a space dependent factor. The fact that these are coupled in a multiplicative manner makes the reconstruction particularly challenging. Several situations of overposed data are considered: boundary observations over a time interval, interior observations at final time, as well as a combination thereof. We devise fixed point schemes and also describe application of a frozen Newton method. In the final time data case we prove convergence of the fixed point scheme as well as uniqueness of both coefficients. Numerical experiments illustrate performance of the reconstruction methods, in particular dependence on the differentiation order in the subdiffusion equation.

In this work, we present and analyze a numerical method for recovering a piecewise constant conductivity with multiple phases in inverse conductivity problems. Specifically, we consider two types of inverse conductivity problems: problems with boundary measurements or with internal measurements. The conductivity is assumed to be constant in each phase, and a Voronoi diagram generated by a set of sites is used to model the phases. An optimization problem with respect to the position of the sites is described to approximate the solution of the inverse problem. Combining techniques from non-smooth shape calculus and the sensitivity of Voronoi diagrams, we prove shape differentiability and compute the gradient of the cost function. Two different formulas for the gradient, a volumetric one and an interface one, are provided. The dependence of the reconstruction on the problem parameters, such as noise, number of sites, and initialization, is investigated through several numerical experiments.

We are concerned with the inverse boundary problem of determining anomalies associated with a semilinear elliptic equation of the form $-\Delta u+a(\mathbf{x}, u) = 0$, where $a(\mathbf{x}, u)$ is a general nonlinear term that belongs to a Hölder class. It is assumed that the inhomogeneity of $a(\mathbf{x}, u)$ is contained in a bounded domain D in the sense that outside D, $a(\mathbf{x}, u) = \lambda u$ with $\lambda\in\mathbb{C}$. We establish novel unique identifiability results in several general scenarios of practical interest. These include determining the support of the inclusion (i.e. D) independent of its content (i.e. $a(\mathbf{x}, u)$ in D) by a single boundary measurement; and determining both D and $a(\mathbf{x}, u)|_D$ by M boundary measurements, where $M\in\mathbb{N}$ signifies the number of unknown coefficients in $a(\mathbf{x}, u)$. The mathematical argument is based on microlocally characterizing the singularities in the solution u induced by the geometric singularities of D, and does not rely on any linearization technique.

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).

This article investigates reduced-order models (ROMs) for efficiently solving geometric inverse source problems in parabolic equations. To reconstruct source supports in diffusion processes, a reduced-order approach combining proper orthogonal decomposition (POD) and incremental singular value decomposition (ISVD) is proposed. This method significantly reduces the computational complexity and storage requirements typically associated with numerical shape and topology optimization. Numerical experiments are conducted to validate the effectiveness and efficiency of the proposed methodology.

Traction force microscopy is a method widely used in biophysics and cell biology to determine forces that biological cells apply to their environment. In the experiment, the cells adhere to a soft elastic substrate, which is then deformed in response to cellular traction forces. The inverse problem consists in computing the traction stress applied by the cell from microscopy measurements of the substrate deformations. In this work, we consider a linear model, in which 3D forces are applied at a 2D interface, called 2.5D traction force microscopy, and a nonlinear pure 2D model, from which we directly obtain a linear pure 2D model. All models lead to a linear resp. nonlinear parameter identification problem for a boundary value problem of elasticity. We analyze the respective forward operators and conclude with some numerical experiments for simulated and experimental data.

This article studies the inverse problem of determining the discontinuity jump of the displacement field, subject to a steady linear elasticity equation, on an interface from boundary measurements of displacements and tractions on boundary subdomains located on different sides of the interface. As a partial result, a conditional stability result of the discontinuity jump through internal measurements of the displacements is established. A conditional stability result for the discontinuity jump from boundary measurements is derived by a domain extension. Finally, an optimal control problem is presented to recover a tangential discontinuity jump, together with a numerical experiment to illustrate the theoretical findings.

In this paper, we study the affine phase retrieval problem, which aims to recover signals from the magnitudes of affine measurements. We develop second-order optimization methods based on Newton and Gauss-Newton iterations and establish that, under specific a priori conditions, the problem exhibits strong convexity. Theoretically, we prove that the Newton method with resampling achieves global quadratic convergence in the noiseless setting for both Gaussian measurements and admissible coded diffraction patterns (CDPs). Furthermore, we demonstrate that the same theoretical framework naturally extends to the Gauss-Newton method, implying its quadratic convergence. To validate our theoretical findings, we conduct extensive numerical experiments. The results confirm the quadratic convergence of second-order methods, while their computational efficiency remains comparable to that of first-order methods. Additionally, our experiments demonstrate that second-order methods achieve exact recovery with relatively few measurements, highlighting their practical feasibility and robustness.

The inverse problem of identifying an unknown space-dependent potential coefficient in the parabolic equation is considered from the additional observation at the terminal time in this work. A novel conditional stability estimate is established for a large terminal time $T$ with suitable assumptions on the input data. Then the potential coefficient and solution of the parabolic equation are parameterized by separate deep neural networks (DNNs), and a new loss function is proposed to reconstruct the unknown potential coefficient. The DNN approximations of the potential coefficient for both continuous and empirical loss functions are analyzed rigorously via utilizing analogous arguments for the conditional stability. Meanwhile, the error estimates are expressed explicitly by the noise level and neural network architectural parameters, which yields a prior rule for determining the number of observations and choosing the size of neural networks. Some numerical experiments are provided to illustrate the robustness of the approach against various levels of measured observation and the accuracy of the numerical solutions.

Traction force microscopy is a method widely used in biophysics and cell biology to determine forces that biological cells apply to their environment. In the experiment, the cells adhere to a soft elastic substrate, which is then deformed in response to cellular traction forces. The inverse problem consists in computing the traction stress applied by the cell from microscopy measurements of the substrate deformations. In this work, we consider a linear model, in which 3D forces are applied at a 2D interface, called 2.5D traction force microscopy, and a nonlinear pure 2D model, from which we directly obtain a linear pure 2D model. All models lead to a linear resp. nonlinear parameter identification problem for a boundary value problem of elasticity. We analyze the respective forward operators and conclude with some numerical experiments for simulated and experimental data.

In this paper we apply the stochastic variance reduced gradient (SVRG) method, which is a popular variance reduction method in optimization for accelerating the stochastic gradient method, to solve large scale linear ill-posed systems in Hilbert spaces. Under {\it a priori} choices of stopping indices, we derive a convergence rate result when the sought solution satisfies a benchmark source 
condition and establish a convergence result without using any source condition. To terminate the method in an {\it a posteriori} manner, we consider the discrepancy principle and show that it terminates the method in finitely many iteration steps almost surely. Various numerical results are reported to test the performance of the method.

We present a mathematical framework for viscoelastic full waveform inversion (FWI) in vertically transverse isotropic media. FWI can be formulated as the nonlinear inverse problem of identifying parameters in the underlying attenuating anisotropic wave equation given partial wave field measurements (seismograms). From a mathematical point of view, one has to solve an operator equation for the full waveform forward operator, which is the corresponding parameter-to-state map. We give a rigorous definition of this operator, show its Fr'echet differentiability, and explicitly characterize the adjoint operator of its Fr'echet derivative. Thus, we provide the main ingredients to implement Newton-type/gradient-based regularization schemes for FWI. Our approach can be directly applied to other concepts of anisotropy.

The Monotonicity Principle (MP), stating a monotonic relationship between a material property and a proper corresponding boundary operator, is attracting great interest in the field of inverse problems, because of its fundamental role in developing real time imaging methods. Moreover, under quite general assumptions, a MP for elliptic PDEs with nonlinear coefficients has been established. This MP provided the basis for introducing a new imaging method to deal with the inverse obstacle problem, in the presence of nonlinear anomalies. This constitutes a relevant novelty because there is a general lack of quantitative and physic based imaging method, when nonlinearities are present. The introduction of a MP based imaging method poses a set of fundamental questions regarding the performance of the method in the presence of noise. The main contribution of this work is focused on theoretical aspects and consists in proving that (i) the imaging method is stable and robust with respect to the noise, (ii) the reconstruction approaches monotonically to a well-defined limit, as the noise level approaches to zero, and that (iii) the limit contains the unknown set and is contained in the outer boundary of the unknown set. Results (i) and (ii) come directly from the MP, while result (iii) requires to prove the so-called Converse of the MP, a theoretical result of fundamental relevance to evaluate the ideal (noise-free) performances of the imaging method. The results are provided in a quite general setting for Calderón problem, and proved for three wide classes where the nonlinearity of the anomaly can be either bounded from infinity and zero, or bounded from zero only, or bounded by infinity only. These classes of constitutive relationships cover the wide majority of cases encountered in applications.

We introduce subgradient-based Lavrentiev regularisation of the form $\mathcal{A}\left(u\right) + \alpha \partial \mathcal{R}\left(u\right) \ni f^{\,{\delta}}$ for linear and nonlinear ill-posed problems with monotone operators ${\mathcal{A} }$ and general regularisation functionals ${\mathcal{R} }$. In contrast to Tikhonov regularisation, this approach perturbs the equation itself and avoids the use of the adjoint of the derivative of ${\mathcal{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.

We are concerned with the well-posedness of an inverse problem for determining the wedge-boundary and associated two-dimensional steady supersonic Euler flow past the wedge, provided that the pressure distribution on the boundary-surface of the wedge and the incoming state of the flow are given. We first establish the existence of wedge-boundaries and associated entropy solutions of the inverse problem, when the pressure on the wedge-boundary is larger than that of the incoming flow but less than a critical value, and the total variation of the incoming flow and the pressure distribution is sufficiently small. This is achieved by a careful construction of suitable approximate solutions and corresponding approximate boundaries via developing a wave-front tracking algorithm and the rigorous proof of their strong convergence subsequentially to a global entropy solution and a wedge-boundary, respectively. Then we establish the $ L^{\infty}$--stability of the wedge-boundaries, by introducing a modified Lyapunov functional for two different solutions with two distinct boundaries, each of which may contain a strong shock-front. The modified Lyapunov functional is carefully designed to control the distance between the two boundaries and is proved to be Lipschitz continuous with respect to the differences of the incoming flow and the pressure on the wedge, which leads to the existence of the Lipschitz semigroup as a converging limit of the approximate solutions and boundaries. Finally, when the pressure distribution on the wedge-boundary is sufficiently close to that of the incoming flow, using this semigroup, we compare two solutions of the inverse problem in the respective supersonic full Euler flow and potential flow and prove that, at $x>0$, the distance between the two boundaries and the difference of the two solutions are of the same order of $x$ multiplied by the cube of the perturbations of the initial boundary data in $L^\infty\cap BV$.

Landweber-type methods are prominent for solving ill-posed inverse problems in Banach spaces and their convergence has been well-understood. However, how to derive their convergence rates remains a challenging open question. In this paper, we tackle the challenge of deriving convergence rates for Landweber-type methods applied to ill-posed inverse problems, where forward operators map from a Banach space to a Hilbert space. Under a benchmark source condition, we introduce a novel strategy to derive convergence rates when the method is terminated by either an a priori stopping rule or the discrepancy principle. Our results offer substantial flexibility regarding step sizes, by allowing the use of variable step sizes. By extending the strategy to deal with the stochastic mirror descent method for solving nonlinear ill-posed systems with exact data, under a benchmark source condition we also obtain an almost sure convergence rate in terms of the number of iterations.

Computational inverse problems utilize a finite number of measurements to infer a discrete approximation of the unknown parameter function. With motivation from the setting of PDE-based optimization, we study the unique reconstruction of discretized inverse problems by examining the positivity of the Hessian matrix. What is the reconstruction power of a fixed number of data observations? How many parameters can one reconstruct? Here we describe a probabilistic approach, and spell out the interplay of the observation size (r) and the number of parameters to be uniquely identified (m). The technical pillar here is the random sketching strategy, in which the matrix concentration inequality and sampling theory are largely employed. By analyzing a randomly subsampled Hessian matrix, we attain a well-conditioned reconstruction problem with high probability. Our main theory is validated in numerical experiments, using an elliptic inverse problem as an example.

We study continuous-domain linear inverse problems that involve a general data-fidelity term and a regularisation term. We consider a regularisation that is formed by the sparsity-promoting total-variation norm, pre-composed with a differential operator that specifies some underlying dictionary of atoms. It has been previously shown that such problems have sparse spline solutions with adaptive knots. These knots are part of the parameterization of the solution and their estimation is itself a difficult non-convex problem. To alleviate this difficulty, we rely on an exact discretization of the optimization problem, where the spline knots are chosen on a dense regular grid. We then follow a multiresolution strategy to refine this grid. In this work, we investigate the convergence of the discretization to the original continuous-domain problem when the grid goes from coarse to fine. We provide an in-depth study of this convergence, concluding that its strength depends on the regularity of the Green's function of the differential operator. We show that uniform convergence holds in very general settings. We carry a numerical analysis to illustrate our theoretical results.

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 cannot 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.