Least-Squares Method for Inverse Medium Problems

We present a two-stage least-squares method to inverse medium problems of reconstructing multiple unknown coefficients simultaneously from noisy data. A direct sampling method is applied to detect the location of the inhomogeneity in the first stage, while a total least-squares method with mixed regularization is used to recover the medium profile in the second stage. The total least-squares method is designed to minimize the residual of the model equation and the data fitting, along with an appropriate regularization, in an attempt to significantly improve the accuracy of the approximation obtained from the first stage. We shall also present an analysis on the well-posedness and convergence of this algorithm. Numerical experiments are carried out to verify the accuracies and robustness of this novel two-stage least-squares algorithm, with great tolerance of noise.



There are no comments yet.


page 14

page 15

page 16

page 17

page 18


Quantum-inspired algorithm for truncated total least squares solution

Total least squares (TLS) methods have been widely used in data fitting....

Convergence analysis of a Crank-Nicolson Galerkin method for an inverse source problem for parabolic systems with boundary observations

This work is devoted to an inverse problem of identifying a source term ...

MINRES for second-order PDEs with singular data

Minimum residual methods such as the least-squares finite element method...

Gradient Descent-based D-optimal Design for the Least-Squares Polynomial Approximation

In this work, we propose a novel sampling method for Design of Experimen...

Recovering Multiple Fractional Orders in Time-Fractional Diffusion in an Unknown Medium

In this work, we investigate an inverse problem of recovering multiple o...

Determining kernels in linear viscoelasticity

In this work, we investigate the inverse problem of determining the kern...

Adaptive spectral decompositions for inverse medium problems

Inverse medium problems involve the reconstruction of a spatially varyin...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In this work, we propose a total least-squares formulation for recovering multiple medium coefficients for a class of inverse medium problems that are governed by the forward model of the form:


where is a bilinear operator on , is the state variable, while represents one or multiple unknown coefficients in the model that are to be recovered, under some measurement data of . Here , and are three Hilbert spaces, and is an observation map from to .

In many applications, it is often required to recover multiple coefficients simultaneously. For instance, the diffusive optical tomography (DOT) aims at recovering the diffusion and absorption coefficients and from the governing equation [1, 2]:


using the Cauchy data collected at the boundary of :


Another example would be the inverse electromagnetic medium problem to recover the unknown magnetic and electric coefficients and in the Maxwell’s system [3, 4]:

using some electric or magnetic measurement data or .

The inverse reconstruction of multiple medium coefficients is generally much more technical and difficult than the single coefficient case. We shall proposed a total least-squares formulation with appropriate regularization to transform the inverse problem into an optimization problem. The total least-squares philosophy is not uncommon. One conventional approach for inverse medium problems is to seek an optimal parameter from a feasible set such that it minimizes an output least-squares functional of the form


where solves the forward model (1.1) when is given, is a regularization term and is a regularization parameter. We refer the readers to [2, 5, 6] for more details about this traditional approach. A relaxed variational method of the least-squares form was proposed and studied in [7, 8] for the impedance computed tomography. The least-squares functional consists of a residual term associated with the governing equation while the measurement data is enforced at the feasible set. Different from the aforementioned approaches, we shall follow some basic principle of a total least-squares approach from [9] and treat the governing equation and the data fitting separately, along with a regularization. That is, we look for an optimal parameter from and a state variable from together such that they minimize an extended functional of the form


This functional combines the residual of model equation, data fitting and constraints on parameters in the least-squares criterion. The combination in (1.5) results in a regularization effect to treat the model equation and allows a robustness as a quasi-reversibility method [10]. Compared with the conventional approaches, the domain of is much more regular and the semi-norm defined by the formulation is much stronger. More precisely, the total least-squares formulation aims to find a solution pair simultaneously in a smooth class and is less sensitive to the noise and uncertainty in the inverse model. Another important feature of this formulation is that the functional is quadratical and convex with respect to each variable and , if the regularization is chosen to be quadratical and convex, while the traditional one in (1.4) is highly nonlinear and non-convex in general. This special feature facilitates us naturally to minimize the functional effectively by the alternating direction iterative (ADI) method [11, 12] so that only two quadratical and convex suboptimizations are required in terms of the variables and respectively at each iteration.

In addition to the functional (1.5) that uses the residual of the forward model (1.1), we will also address another least-squares functional that makes use of the equivalent first-order system of the forward model (1.1) and replaces the first term in (1.5) by the residuals of the corresponding first-order system. Using the first-order system has been the fundamental idea in modern least-squares methods in solving second-order PDEs [13, 14, 15]. The advantages of using first-order formulations are much more significant to the numerical solutions of inverse problems, especially when we aim at simultaneously reconstructing multiple coefficients as we do in this work. First, the multiple coefficients appear in separated first-order equations, hence are naturally decoupled. This would greatly reduce the nonlinearity and enhance the convexity of the resulting optimization systems. Second, the first-order formulation relaxes the regularity requirement of the solutions in the resulting analysis.

A crucial step to an effective reconstruction of multiple coefficients is to seek some reasonable initial approximations to capture some general (possibly rather inaccurate) geometric and physical profiles of all the unknown multiple coefficients. This is a rather technical and difficult task in numerical solutions. For this purpose, we shall propose to adopt the direct sampling-type method (DSM) that we have been developing in recent years (cf. [17, 18, 26, 27]). Using the index functions provided by DSM, we shall determine a computational domain that is often much smaller than the original physical domain, then the restricted index functions on the computational domain serve as the initial guesses of the unknown coefficients. In this work, we will apply a newly developed DSM [19], where two groups of probing and index functions are constructed to identify and decouple the multiple inhomogeneous inclusions of different physical nature, which is different from the classical DSMs targeting the inhomogeneous inclusions of one single physical nature. As we shall see, DSMs turn out to be very effective and fast solvers to provide some reasonable initial approximations.

The rest of the paper is structured as follows. In Section 2, we justify the well-posedness of the least-squares formulation for the general inverse medium problems. In Section 3, we propose an alternating direction iterative method for solving the minimization problem and prove the convergence of the ADI method. We illustrate in Section 4 how this total least-squares method applies to a concrete inverse problem, by taking the DOT problem as a benchmark problem. We present numerical results in Section 5 for a couple of different types of inhomogeneous coefficients for the DOT problem to demonstrate the stability and effectiveness of this proposed method. Throughout the paper, , and denote generic constants which may differ at each occurrence.

2 Well-posedness of the least-squares formulation for inverse medium problem

Recall that to solve the inverse medium problems modeled by (1.1), we propose the following least-squares formulation


This section is devoted to the well-posedness of the total least-squares formulation (2.1), namely, the existence of a solution to (2.1) and the conditional stability of the reconstruction with respect to the measurement. To provide a rigorous justification of the well-posedneness, we present several assumptions on the least-squares formulation, which are minimal for the proof. We will verify these checkable assumptions in Section 4 for a concrete example of such inverse medium problems.

Let us first introduce several notations. For simplicity, for a given (resp. ), we will write (resp. ) as


We denote the subdifferential of the regularization term at by , and denote the inner products of the Hilbert spaces , and by , and respectively.

2.1 Existence of a minimizer

We present the following assumptions on the regularization term and operators and in the forward model:

Assumption 1.

The regularization term is strictly convex and weakly lower semicontinuous. Furthermore, is also coercive [6], i.e., .

This assumption implies that the level set defines a bounded set in .

Assumption 2.

Given a constant , for in the level set , , where with a specific side constraint that is not in the least-squares formulation (2.1), is a closed linear operator and is uniformly coercive, i.e., the graph norm satisfies uniformly in for some constant , and thus has a unique solution in .

Under Assumption 2, we can define the inverse operator , which is uniformly bounded by the coercivity of . We also need the following assumption on the sequentially closedness of operators and .

Assumption 3.

The operators and are weakly sequentially closed, i.e., if a sequence converges to weakly in , then the sequence converges to weakly in and the sequence converges to weakly in .

Then we can verify the existence of the minimizers to the least-squares formulation (2.1).

Theorem 1.

Under Assumptions 13, there exists a minimizer in of the least-squares formulation (2.1).


Since and are nonempty, there exists a minimizing sequence in such that


By Assumptions 12, is a coercive functional and the graph norm is uniformly coercive, thus it follows (1.5) that the sequence is uniformly bounded. Then there exists a subsequence of , still denoted by , and some such that converges to weakly in and converges to weakly in . As and are weakly sequentially closed by Assumption 3, there hold


From the weak lower semicontinuity of the norms and , we have

Together with the lower semicontinuity of the regularization term , we can deduce that

Hence it follows (2.3) that is indeed a minimizer of the functional in . ∎

Remark 1.

While the weakly sequentially closedness in Assumption 3 is nontrivial to verify for most nonlinear inverse medium problems, one can reach the weak convergence (2.4) with the compactness of and the concrete formula of , as shown in Section 4.

2.2 Conditional stability

In this subsection, we present some conditional stability estimates of the total least-squares formulation (

2.1) for the general inverse medium problems. First we introduce two pairs and that satisfy


Letting be the unique minimizer of (2.1) in a neighborhood of , we study the approximation error to illustrate the stability of the least-squares formulation (2.1) with respect to the measurement and also the term in the governing equation (1.1). Denote the residual of the governing equation by . As is the local minimizer of functional in (1.5), we have . Therefore, by the definition of we have the inequality


which directly leads to the following observation on :


If is coercive, (2.7) provides a rough estimate of the reconstruction with respect to the data noise.

We can further derive an estimate of the approximation error , under the following assumption on the operator .

Assumption 4.

There exists a norm on such that for any ,

Assumption 4 holds when belongs to a finite rank subspace of and for all non-zero . We can now deduce the following result of the approximation error .

Lemma 1.

Under Assumptions 14, the approximation error is bounded in -norm by the data noise and the regularization term, i.e.,


Using the bilinear property of , one can rewrite the difference as


By Assumption 2, admits an inverse operator from to , which, together with (2.5), (2.9) and the definition of , implies


Plugging (2.10) into (2.6) leads to an inequality:


It follows Assumption 4 that there exists a norm such that


Then we can deduce from (2.11), the triangle inequality, the boundedness of and (2.12) that

where is a constant, which completes the proof. ∎

The rest of this section is devoted to verifying the consistency of the least-squares formulation (2.1) as the noise level of measurement goes to zero, which is an essential property of a regularization scheme. If we choose an appropriate regularization parameter according to the noise level of the data, we can deduce the convergence result of the reconstructed coefficients associated with the regularization parameter . More precisely, given a set of exact data , we consider a parametric family such that . In the rest of this section, we denote the functional in (1.5) with and by , and the minimizer of by . Then we justify the consistency of the least-squares formulation (2.1) by proving the convergence of the sequence of minimizers to the minimum norm solution [6] of the system (2.5) as .

Theorem 2.

Let be a sequence converging to zero, and be the corresponding sequence of minimizers of . Then under Assumptions 13, has a subsequence that converges weakly to a minimum norm solution of the system (2.5), i.e.,


As is the minimizer of , there holds


By definition, , and thus is uniformly bounded. Following the similar argument in the proof of Theorem 1, there exist a subsequence of , still denoted as , and some such that converges to weakly in . By Assumption 3, we have

From (2.13) one can also derive that


which implies

Therefore, as , will converge to satisfying


Recall that one has an estimate of from (2.14) that

which leads to

as . Using the lower semicontinuity of functional , one obtain that

Together with (2.15), we conclude that is a minimum norm solution of (2.5). ∎

3 ADI method and convergence analysis

An important feature of the least-squares formulation (2.1) is that the functional is quadratical and convex with respect to each variable and . This unique feature facilitates us naturally to minimize the functional effectively by the alternating direction iterative (ADI) method [11, 12] so that only two quadratical and convex suboptimizations of one variable are required at each iteration. We shall carry out the convergence analysis of the ADI method in this section for general inverse medium problems.

Alternating direction iterative method for the minimization of (1.5).

Given an initial pair , find a sequence of pairs for as below:

  • Given , find by solving

  • Given , find by solving


We shall establish the convergence of the sequence generated by the ADI method, under Assumptions 13 on the least-squares formulation (2.1). For this purpose, we would like to introduce the Bregman distance [21] with respect to ,


which is always nonnegative for convex function . Now we are ready to present the following lemma on convergence of the sequence generated by (3.1)–(3.2).

Lemma 2.

Under Assumptions 13, the sequence generated by the ADI method (3.1)–(3.2) converges to a pair that satisfies the optimality condition of (2.1):


Using the optimality condition satisfied by the minimizer of (3.1), we deduce


Similarly, from the optimality condition satisfied by the minimizer of (3.2), we obtain

Taking , and in (3.3), we can derive that


The following equality hold for the first term in (3.6):


Plugging (3.7) into (3.6), it follows that


As the sequence is generated by ADI method (3.1)–(3.2), in each iteration the updated (resp. ) minimizes the functional (resp. ), which would lead to

for all . Then we further derive from (3.5) and (3.8) that for any , satisfies


which implies that is bounded. Then we can conclude using Assumption 2 that forms a Cauchy sequence and thus converges to some . Since is uniformly bounded for all , we can derive that converges to some from the strict convexity of . As the sequence is monotone decreasing, there exists a limit . Following the similar argument in the proof of Theorem 1, we conclude that by Assumption 3. This completes the proof of the convergence. ∎

Remark 2.

If (3.4) has a unique solution in a neighborhood of initial guess , then the solution is a local minimizer of the least-squares formulation (2.1), and we can apply the ADI method to generate a sequence that converges to this local minimizer as a plausible approximation of the exact coefficients.

4 Diffusive optical tomography

We take the diffusive optical tomography (DOT) as an example to illustrate the total least-squares approach for the concrete inverse medium problem in this section. We will introduce the mixed regularization term, present the least-squares formulation of the first-order system of DOT, and then verify the assumptions in Section 2 for the proposed formulation. We shall use the standard notations for Sobolev spaces. The objective of the DOT problem is to determine the unknown diffusion and absorption coefficients and simultaneously in a Lipschitz domain () from the model equation


with a pair of Cauchy data on the boundary , i.e.,

Throughout this section, we shall use the notation in the total least-squares formulation, where denotes the Neumann to source map. To define the Neumann to source map, we first introduce the boundary restriction mapping on , i.e., denotes the boundary value of . Then we use to denote the trace operator [22], which is formally defined to be the unique linear continuous extension of as an operator from onto . Using Riesz representation theorem, there exists a function in , denoted by , such that for any ,

which will be denoted by in the least-squares formulation.

4.1 Mixed regularization

In this subsection, we present the mixed regularization term for the DOT problem. As the regularization term in the least-squares formulation (2.1) shall encode the priori information, e.g., sparsity, continuity, lower or upper bound and other properties of the unknown coefficients, it is essential to choose an appropriate regularization term for a concrete inverse problem to ensure satisfactory reconstructions. In this work, we introduce a mixed regularization term for a coefficient :


In this revision we change the semi-norm of into the full norm in the regularization to ensure the strict convexity of . where is given by

and and are the lower and upper bounds of the coefficient respectively. The first term is the regularization term, the second term corresponds to the regularization, and the third term enforces the reconstruction to meet the constrains of the coefficient.

In practice, the regularization enhances the reconstruction and helps find geometrically sharp inclusions, but might cause spiky results. The regularization generates the reconstructions with overall clear structures, while the retrieved background may be blurry. Compared with other more conventional regularization methods, this mixed regularization technique in (4.2) combines two penalty terms and effectively promotes multiple features of the target solution, that is, it enhances the sparsity of the solution while it still preserves the overall structure at the same time. The scalar parameters and need to be chosen carefully to compromise the two regularization terms.

4.2 First-order formulation of DOT

As we have emphasized in the Introduction, it may have some advantages to make use of the residuals of the first-order system of the model equation (4.1), instead of the residual of the original equation in the formulation (2.1), when we aim at recovering two unknown coefficients and simultaneously. Similarly to the formulation (2.1), we now have , and the operator is given by


where we have introduced an auxiliary vector flux

and each entry of is of a first-order form such that two coefficients are separated naturally. Clearly, is still bilinear with respect to the state variables and coefficients . Using the first-order system, we can then come to the following total least-squares functional:


where is the trace operator, and and are the corresponding mixed regularization terms of and defined as in (4.2), , are lower and upper bounds of , and , are lower and upper bounds of . We shall minimize (4.4) over , that is, , .

We will apply the ADI method to solve the least-squares formulation of and :


Given an initial guess , we find a sequence for as below:

  • Given , , find , by solving

  • Given , , find , by solving

4.3 Well-posedness of the least-squares formulation for DOT

Recall that we have proved the well-posedness of the least-squares formulation in Section 2 for general inverse medium problems. This part is devoted to the verification of assumptions in Section 2 for the formulation (4.5) to ensure its well-posedness.

Firstly we consider Assumption 1 on the regularization terms. It is observed from the formula (4.2) that each term of and