Least-Squares Method for Inverse Medium Problems

Kazufumi Ito¹¹1 Department of Mathematics and Center for Research in Scientific Computation, North Carolina State University, Raleigh, NC 27695. ([email protected]). Ying Liang²²2Department of Mathematics, Purdue University, West Lafayette, IN 47907. ([email protected]). Jun Zou³³3Department of Mathematics, The Chinese University of Hong Kong, Shatin, N.T., Hong Kong. The work of this author was substantially supported by Hong Kong RGC General Research Fund (projects 14306719 and 14306718). ([email protected]).

Abstract

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.

Keywords: Inverse medium problem, least-squares method, reconstruction algorithm, mixed regularization

Mathematics Subject Classification (MSC2000): 35R30, 65F22, 65N21, 65N06

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:

L(u,q)=g

(1.1)

where $L$ is a bilinear operator on $(u,q)$ , $u\in X$ is the state variable, while $q\in Z$ represents one or multiple unknown coefficients in the model that are to be recovered, under some measurement data of $f:=Cu\in Y$ . Here $X$ , $Y$ and $Z$ are three Hilbert spaces, and $C$ is an observation map from $X$ to $Y$ .

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 $\sigma$ and $\mu$ from the governing equation [1, 2]:

-\nabla\cdot(\sigma(x)\nabla u)+\mu(x)\,u=0\,\quad\mbox{in}~{}~{}\Omega

(1.2)

using the Cauchy data $(f,h)$ collected at the boundary $\Gamma$ of $\Omega$ :

u|_{\Gamma}=f,\quad\frac{\partial u}{\partial\nu}|_{\Gamma}=h\,.

(1.3)

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

\begin{array}[]{l}\nabla\times\vec{H}+i\omega\mu(x)\,\vec{E}=0\,\quad\mbox{in}~{}~{}\Omega\,,\\ \nabla\times\vec{E}-i\omega\lambda(x)\,\vec{H}=0\,\quad\mbox{in}~{}~{}\Omega\,,\end{array}

using some electric or magnetic measurement data $\vec{E}$ or $\vec{H}$ .

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 $q$ from a feasible set $K\subset Z$ such that it minimizes an output least-squares functional $j$ of the form

j(q)=\|Cu(q)-f\|_{Y}^{2}+\alpha\,\psi(q),

(1.4)

where $u(q)$ solves the forward model (1.1) when $q$ is given, $\psi$ is a regularization term and $\alpha>0$ 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 $q$ from $Z$ and a state variable $u$ from $X$ together such that they minimize an extended functional of the form

J(u,q)=\|L(u,q)-g\|_{X}^{2}+\|Cu-f\|_{Y}^{2}+\alpha\,\psi(q)\,.

(1.5)

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 $J(u,q)$ 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 $(u,q)$ 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 $J(u,q)$ is quadratical and convex with respect to each variable $u$ and $q$ , if the regularization $\psi$ is chosen to be quadratical and convex, while the traditional one $j(q)$ in (1.4) is highly nonlinear and non-convex in general. This special feature facilitates us naturally to minimize the functional $J(u,q)$ 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 $u$ and $q$ 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, $c$ , $c_{0}$ and $c_{1}$ 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

\min_{u\in X,q\in Z}J(u,q)=\|L(u,q)-g\|_{X}^{2}+\|Cu-f\|_{Y}^{2}+\alpha\,\psi(q)\,.

(2.1)

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 $q\in Z$ (resp. $u\in X$ ), we will write $L_{q}$ (resp. $\Phi_{u}$ ) as

L_{q}u:=L(u,q)~{}~{}(\mbox{resp.}\,\Phi_{u}q:=L(u,q))\,.

(2.2)

We denote the subdifferential of the regularization term $\psi$ at $q$ by $\partial\psi(q)$ , and denote the inner products of the Hilbert spaces $X$ , $Y$ and $Z$ by $(\cdot,\cdot)_{X}$ , $(\cdot,\cdot)_{Y}$ and $\langle\cdot,\cdot\rangle$ respectively.

2.1 Existence of a minimizer

We present the following assumptions on the regularization term $\psi$ and operators $L$ and $C$ in the forward model:

Assumption 1.

The regularization term $\psi$ is strictly convex and weakly lower semicontinuous. Furthermore, $\psi$ is also coercive [6], i.e., $\psi(q)\geq c\|q\|^{2}_{Z}$ .

This assumption implies that the level set $\{q\in Z:\psi(q)\leq c_{0}\}$ defines a bounded set in $Z$ .

Assumption 2.

Given a constant $c_{0}$ , for $q$ in the level set $\{q\in Z:\psi(q)\leq c_{0}\}$ , $L_{q}:dom(L)\rightarrow X$ , where $dom(L)\subset\{u\in X:L(u,q)\in X\text{ and }Cu\in Y\text{ for all }q\in Z\text{ satisfying }\psi(q)\leq c_{0}\}$ 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 $|u|_{W,q}:=\|L(u,q)\|_{X}$ satisfies $|u|_{W,q}\geq c_{1}\|u\|_{X}$ uniformly in $q$ for some constant $c_{1}>0$ , and thus $L(u,q)=g\in X$ has a unique solution in $dom(L)$ .

Under Assumption 2, we can define the inverse operator $L_{q}^{-1}:X\rightarrow dom(L)$ , which is uniformly bounded by the coercivity of $L_{q}$ . We also need the following assumption on the sequentially closedness of operators $L$ and $C$ .

Assumption 3.

The operators $L$ and $C$ are weakly sequentially closed, i.e., if a sequence $\{(u_{n},q_{n})\}_{n=1}^{\infty}$ converges to $(u,q)$ weakly in $X\times Z$ , then the sequence $\{L(u_{n},q_{n})\}_{n=1}^{\infty}$ converges to $L(u,q)$ weakly in $X$ and the sequence $\{Cu_{n}\}_{n=1}^{\infty}$ converges to $Cu$ weakly in $Y$ .

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

Theorem 1.

Under Assumptions 1–3, there exists a minimizer $(u^{\star},q^{\star})$ in $X\times Z$ of the least-squares formulation (2.1).

Proof.

Since $X$ and $Z$ are nonempty, there exists a minimizing sequence $\{(u_{n},q_{n})\}_{n=1}^{\infty}$ in $X\times Z$ such that

\lim_{n\rightarrow\infty}J(u_{n},q_{n})=\inf_{(u,q)\in X\times Z}J(u,q).

(2.3)

By Assumptions 1–2, $\psi$ is a coercive functional and the graph norm $|\cdot|_{W,q}$ is uniformly coercive, thus it follows (1.5) that the sequence $\{(u_{n},q_{n})\}_{n=1}^{\infty}$ is uniformly bounded. Then there exists a subsequence of $\{(u_{n},q_{n})\}_{n=1}^{\infty}$ , still denoted by $\{(u_{n},q_{n})\}_{n=1}^{\infty}$ , and some $(u^{\star},q^{\star})\in X\times Z$ such that $u_{n}$ converges to $u^{\star}$ weakly in $X$ and $q_{n}$ converges to $q^{\star}$ weakly in $Z$ . As $L$ and $C$ are weakly sequentially closed by Assumption 3, there hold

\displaystyle L(u_{n},q_{n})\text{ converges to }L(u^{\star},q^{\star})\text{ weakly in }X;\quad Cu_{n}\text{ converges to }Cu^{\star}\text{ weakly in }Y.

(2.4)

From the weak lower semicontinuity of the norms $\|\cdot\|_{X}$ and $\|\cdot\|_{Y}$ , we have

\|L(u^{\star},q^{\star})-g\|^{2}_{X}+\|Cu^{\star}-f\|^{2}_{Y}\leq\lim_{n\rightarrow\infty}\inf\left(\|L(u_{n},q_{n})-g\|^{2}_{X}+\|Cu_{n}-f\|^{2}_{Y}\right).

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

J(u^{\star},q^{\star})\leq\lim_{n\rightarrow\infty}\inf J(u_{n},q_{n}).

Hence it follows (2.3) that $(u^{\star},q^{\star})$ is indeed a minimizer of the functional $J$ in $X\times Z$ . ∎

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 $dom(L)$ and the concrete formula of $L$ , 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 $(\bar{u},\bar{q})$ and $(\bar{g},\bar{f})$ that satisfy

L(\bar{u},\bar{q})=\bar{g},\;\;C\bar{u}=\bar{f}.

(2.5)

Letting $(u^{\star},q^{\star})$ be the unique minimizer of (2.1) in a neighborhood of $(\bar{u},\bar{q})$ , we study the approximation error $(\delta u,\delta q):=(u^{\star}-\bar{u},q^{\star}-\bar{q})$ to illustrate the stability of the least-squares formulation (2.1) with respect to the measurement $f$ and also the term $g$ in the governing equation (1.1). Denote the residual of the governing equation by $\epsilon_{1}:=L(u^{\star},q^{\star})-g$ . As $(u^{\star},q^{\star})$ is the local minimizer of functional $J$ in (1.5), we have $J(u^{\star},q^{\star})\leq J(\bar{u},\bar{q})$ . Therefore, by the definition of $J$ we have the inequality

\|\epsilon_{1}\|_{X}^{2}+\|C(u^{\star}-\bar{u})-(f-\bar{f})\|_{Y}^{2}+\alpha\psi(q^{\star})\leq\|g-\bar{g}\|_{X}^{2}+\|f-\bar{f}\|_{Y}^{2}+\alpha\psi(\bar{q}),

(2.6)

which directly leads to the following observation on $\psi(q^{\star})$ :

\alpha\psi(q^{\star})\leq\|g-\bar{g}\|_{X}^{2}+\|f-\bar{f}\|_{Y}^{2}+\alpha\,\psi(\bar{q}).

(2.7)

If $\psi$ is coercive, (2.7) provides a rough estimate of the reconstruction $q^{\star}$ with respect to the data noise.

We can further derive an estimate of the approximation error $\delta q$ , under the following assumption on the operator $L$ .

Assumption 4.

There exists a norm $\|\cdot\|_{W}$ on $Z$ such that for any $q\in Z$ ,

\|CL_{q}^{-1}\Phi_{\bar{u}}(q-\bar{q})\|^{2}_{Y}\geq\|q-\bar{q}\|_{W}^{2}.

Assumption 4 holds when $\delta q$ belongs to a finite rank subspace $S$ of $Z$ and $\|\delta q\|_{W}\neq 0$ for all non-zero $\delta q\in S$ . We can now deduce the following result of the approximation error $\delta q$ .

Lemma 1.

Under Assumptions 1–4, the approximation error $\delta q$ is bounded in $W$ -norm by the data noise and the regularization term, i.e.,

\|\delta q\|_{W}^{2}+\alpha\psi(q^{\star})\leq c_{0}\left(\|g-\bar{g}\|_{X}^{2}+\|f-\bar{f}\|_{Y}^{2}+\alpha\,\psi(\bar{q})\right).

(2.8)

Proof.

Using the bilinear property of $L$ , one can rewrite the difference $L(u^{\star},q^{\star})-L(\bar{u},\bar{q})$ as

L(u^{\star},q^{\star})-L(\bar{u},\bar{q})=L_{q^{\star}}(\delta u)+\Phi_{\bar{u}}(\delta q).

(2.9)

By Assumption 2, $L_{q}$ admits an inverse operator $L_{q}^{-1}$ from $X$ to $dom(L)$ , which, together with (2.5), (2.9) and the definition of $\epsilon_{1}$ , implies

\delta u=L_{q^{\star}}^{-1}(\epsilon_{1}+g-\bar{g}-\Phi_{\bar{u}}(\delta q)).

(2.10)

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

\displaystyle\|CL_{q^{\star}}^{-1}\Phi_{\bar{u}}\delta q-(\epsilon_{1}+L_{q^{\star}}^{-1}(g-\bar{g}))+f-\bar{f}\|_{Y}^{2}+\|\epsilon_{1}\|_{X}^{2}+\alpha\,\psi(q^{\star})\leq\|g-\bar{g}\|_{X}^{2}+\|f-\bar{f}\|_{Y}^{2}+\alpha\,\psi(\bar{q}).

(2.11)

It follows Assumption 4 that there exists a norm $\|\cdot\|_{W}$ such that

\|CL_{q}^{-1}\Phi_{\bar{u}}\delta q\|_{Y}^{2}\geq\|\delta q\|^{2}_{W}.

(2.12)

Then we can deduce from (2.11), the triangle inequality, the boundedness of $L_{q}^{-1}$ and (2.12) that

\displaystyle\|\delta q\|_{W}^{2}+\alpha\,\psi(q^{\star})\leq c_{0}(\|g-\bar{g}\|_{X}^{2}+\|f-\bar{f}\|_{Y}^{2}+\alpha\,\psi(\bar{q})),

where $c_{0}$ 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 $\alpha$ according to the noise level of the data, we can deduce the convergence result of the reconstructed coefficients associated with the regularization parameter $\alpha$ . More precisely, given a set of exact data $(\overline{g},\overline{f})$ , we consider a parametric family $\{(g_{\alpha},f_{\alpha})\}$ such that $\|g_{\alpha}-\overline{g}\|_{X}^{2}+\|f_{\alpha}-\overline{f}\|_{Y}^{2}=o(\alpha)$ . In the rest of this section, we denote the functional $J$ in (1.5) with $g=g_{\alpha}$ and $f=f_{\alpha}$ by $J_{\alpha}$ , and the minimizer of $J_{\alpha}$ by $(u_{\alpha},q_{\alpha})$ . Then we justify the consistency of the least-squares formulation (2.1) by proving the convergence of the sequence of minimizers $\{q_{\alpha}\}$ to the minimum norm solution [6] of the system (2.5) as $\alpha\rightarrow 0$ .

Theorem 2.

Let $\{\alpha_{n}\}_{n=1}^{\infty}\subset\mathbb{R}^{+}$ be a sequence converging to zero, and $\{(u_{\alpha_{n}},q_{\alpha_{n}})\}_{n=1}^{\infty}$ be the corresponding sequence of minimizers of $J_{\alpha_{n}}$ . Then under Assumptions 1–3, $\{(u_{\alpha_{n}},q_{\alpha_{n}})\}_{n=1}^{\infty}$ has a subsequence that converges weakly to a minimum norm solution $(\hat{u},\hat{q})$ of the system (2.5), i.e.,

	$\displaystyle L(\hat{u},\hat{q})=\bar{g},\quad C\hat{u}=\bar{f}.$
	$\displaystyle\psi(\hat{q})\leq\psi(\bar{q})\mbox{ for all $(\bar{u},\bar{q})$ satisfying \eqref{0}}.$

Proof.

As $(u_{\alpha_{n}},q_{\alpha_{n}})$ is the minimizer of $J_{\alpha_{n}}$ , there holds

J_{\alpha_{n}}(u_{\alpha_{n}},q_{\alpha_{n}})\leq J_{\alpha_{n}}(\bar{u},\bar{q}).

(2.13)

By definition, $J_{\alpha_{n}}(\bar{u},\bar{q})=\|L(\bar{u},\bar{q})-g_{\alpha_{n}}\|_{X}^{2}+\|C\bar{u}-f_{\alpha_{n}}\|_{Y}^{2}+\alpha_{n}\,\psi(\bar{q})$ , and thus $\{J_{\alpha_{n}}(u_{\alpha_{n}},q_{\alpha_{n}})\}_{n=1}^{\infty}$ is uniformly bounded. Following the similar argument in the proof of Theorem 1, there exist a subsequence of $\{(u_{\alpha_{n}},q_{\alpha_{n}})\}_{n=1}^{\infty}$ , still denoted as $\{(u_{\alpha_{n}},q_{\alpha_{n}})\}_{n=1}^{\infty}$ , and some $(\hat{u},\hat{q})$ such that $\{(u_{\alpha_{n}},q_{\alpha_{n}})\}_{n=1}^{\infty}$ converges to $(\hat{u},\hat{q})$ weakly in $X\times Z$ . By Assumption 3, we have

\displaystyle L(u_{\alpha_{n}},q_{\alpha_{n}})\text{ converges to }L(\hat{u},\hat{q})\text{ weakly in }X;\quad Cu_{\alpha_{n}}\text{ converges to }C\hat{u}\text{ weakly in }Y.

From (2.13) one can also derive that

\|L(u_{\alpha_{n}},q_{\alpha_{n}})-g_{\alpha_{n}}\|_{X}^{2}+\|Cu_{\alpha_{n}}-f_{\alpha_{n}}\|_{Y}^{2}+\alpha_{n}\,\psi(q_{\alpha_{n}})\leq\|\bar{g}-g_{\alpha_{n}}\|_{X}^{2}+\|\bar{f}-f_{\alpha_{n}}\|_{Y}^{2}+\alpha_{n}\,\psi(\bar{q}),

(2.14)

which implies

\displaystyle L(u_{\alpha_{n}},q_{\alpha_{n}})\text{ converges to }\bar{g}\text{ strongly in }X;\quad Cu_{\alpha_{n}}\text{ converges to }\bar{f}\text{ strongly in }Y.

Therefore, as $\alpha_{n}\rightarrow 0$ , $\{(u_{\alpha_{n}},q_{\alpha_{n}})\}_{n=1}^{\infty}$ will converge to $(\hat{u},\hat{q})$ satisfying

L(\hat{u},\hat{q})=\bar{g},\quad C\hat{u}=\bar{f}.

(2.15)

Recall that one has an estimate of $\psi(q_{\alpha_{n}})$ from (2.14) that

\alpha_{n}\psi(q_{\alpha_{n}})\leq\|g_{\alpha_{n}}-\bar{g}\|_{X}^{2}+\|f_{\alpha_{n}}-\bar{f}\|_{Y}^{2}+\alpha_{n}\,\psi(\bar{q}),

which leads to

\psi(q_{\alpha_{n}})\leq\psi(\bar{q})+o(1),

as $\alpha_{n}\rightarrow 0$ . Using the lower semicontinuity of functional $\psi$ , one obtain that

\psi(\hat{q})\leq\psi(\bar{q})\mbox{ for all $(\bar{u},\bar{q})$ satisfying \eqref{0}}.

Together with (2.15), we conclude that $(\hat{u},\hat{q})$ 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 $J(u,q)$ is quadratical and convex with respect to each variable $u$ and $q$ . This unique feature facilitates us naturally to minimize the functional $J$ 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 $(u_{0},q_{0})$ , find a sequence of pairs $(u_{k},q_{k})$ for $k\geq 1$ as below:

•

Given $q_{k}$ , find $u=u_{k+1}\in X$ by solving

$\min_{u\in X}\|L(u,q_{k})-g\|_{X}^{2}+\|Cu-f\|_{Y}^{2}\,.$ (3.1)
•

Given $u_{k+1}$ , find $q=q_{k+1}\in Z$ by solving

$\min_{q\in Z}\|L(u_{k+1},q)-g\|_{X}^{2}+\alpha\psi(q)\,.$ (3.2)

We shall establish the convergence of the sequence $\{(u_{k},q_{k})\}_{k=1}^{\infty}$ generated by the ADI method, under Assumptions 1–3 on the least-squares formulation (2.1). For this purpose, we would like to introduce the Bregman distance [21] with respect to $\xi\in\partial\psi(p)$ ,

E(q,p):=\psi(q)-\psi(p)-\langle\xi,q-p\rangle~{}~{}\forall\,q,p\in Z\,,

(3.3)

which is always nonnegative for convex function $\psi$ . Now we are ready to present the following lemma on convergence of the sequence $\{(u_{k},q_{k})\}_{k=1}^{\infty}$ generated by (3.1)–(3.2).

Lemma 2.

Under Assumptions 1–3, the sequence $\{(u_{k},q_{k})\}_{k=1}^{\infty}$ generated by the ADI method (3.1)–(3.2) converges to a pair $(u^{\star},q^{\star})$ that satisfies the optimality condition of (2.1):

	$\displaystyle(L_{q})^{}(L_{q}u-g)+C^{}(Cu-f)=0,$		(3.4)
	$\displaystyle-2(\Phi_{u})^{*}(\Phi_{u}q-g)\in\partial\big{(}\alpha\psi(q)\big{)}\,.$		(3.4)

Proof.

Using the optimality condition satisfied by the minimizer $u_{k+1}$ of (3.1), we deduce

$\displaystyle 0$	$\displaystyle=\big{(}L_{q_{k}}^{}(L_{q_{k}}u_{k+1}-g),u_{k+1}-u_{k}\big{)}_{X}+\big{(}C^{}(Cu_{k+1}-f),u_{k+1}-u_{k}\big{)}_{Y}$
	$\displaystyle=\frac{1}{2}\Big{(}\\|L_{q_{k}}u_{k+1}-g\\|_{X}^{2}+\\|Cu_{k+1}-f\\|_{Y}^{2}-\left(\\|L_{q_{k}}u_{k}-g\\|_{X}^{2}+\\|Cu_{k}-f\\|_{Y}^{2}\right)$
	$\displaystyle\quad+\\|L_{q_{k}}(u_{k+1}-u_{k})\\|_{X}^{2}+\\|C(u_{k+1}-u_{k})\\|_{Y}^{2}\Big{)}\,.$	(3.5)

Similarly, from the optimality condition satisfied by the minimizer $q_{k+1}$ of (3.2), we obtain

-2\Phi_{u_{k+1}}^{*}(\Phi_{u_{k+1}}q_{k+1}-g)\in\partial\big{(}\alpha\psi(q_{k+1})\big{)}.

Taking $q=q_{k}$ , $p=q_{k+1}$ and $\xi=-2\Phi_{u_{k+1}}^{*}(\Phi_{u_{k+1}}q_{k+1}-g)$ in (3.3), we can derive that

\langle-2\Phi_{u_{k+1}}^{*}(\Phi_{u_{k+1}}q_{k+1}-g),q_{k}-q_{k+1}\rangle+\alpha\psi(q_{k+1})-\alpha\psi(q_{k})+E(q_{k},q_{k+1})=0.

(3.6)

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

	$\displaystyle\,\,\quad\langle\Phi_{u_{k+1}}^{*}(\Phi_{u_{k+1}}q_{k+1}-g),q_{k+1}-q_{k}\rangle$
	$\displaystyle=\langle\Phi_{u_{k+1}}q_{k+1}-g,\Phi_{u_{k+1}}(q_{k+1}-q_{k})\rangle$
	$\displaystyle=\frac{1}{2}\langle(\Phi_{u_{k+1}}q_{k+1}-g)+(\Phi_{u_{k+1}}q_{k}-g),\Phi_{u_{k+1}}q_{k+1}-g-(\Phi_{u_{k+1}}q_{k}-g)\rangle$
	$\displaystyle\quad+\frac{1}{2}\langle(\Phi_{u_{k+1}}q_{k+1}-g)-(\Phi_{u_{k+1}}q_{k}-g),\Phi_{u_{k+1}}q_{k+1}-g-(\Phi_{u_{k+1}}q_{k}-g)\rangle$
	$\displaystyle=\frac{1}{2}\left(\\|\Phi_{u_{k+1}}q_{k+1}-g\\|_{X}^{2}-\\|\Phi_{u_{k+1}}q_{k}-g\\|_{X}^{2}+\\|\Phi_{u_{k+1}}(q_{k+1}-q_{k})\\|_{X}^{2}\right).$		(3.7)

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

\|\Phi_{u_{k+1}}q_{k+1}-g\|_{X}^{2}+\alpha\psi(q_{k+1})-\left(\|\Phi_{u_{k+1}}q_{k}-g\|_{X}^{2}+\alpha\psi\left(q_{k}\right)\right)+\|\Phi_{u_{k+1}}(q_{k+1}-q_{k})\|_{X}^{2}+E(q_{k},q_{k+1})=0.

(3.8)

As the sequence $\{(u_{k},q_{k})\}_{k=1}^{\infty}$ is generated by ADI method (3.1)–(3.2), in each iteration the updated $u_{k+1}$ (resp. $q_{k+1}$ ) minimizes the functional $J(u,q_{k})$ (resp. $J(u_{k+1},q)$ ), which would lead to

J(u_{k},q_{k})\geq J(u_{k+1},q_{k})\geq J(u_{k+1},q_{k+1})

for all $k\geq 0$ . Then we further derive from (3.5) and (3.8) that for any $m\geq 1$ , $J(u_{m},q_{m})$ satisfies

	$\displaystyle J(u_{m},q_{m})+\sum_{k=0}^{m-1}E(q_{k},q_{k+1})$
	$\displaystyle+\sum_{k=0}^{m-1}\left(\\|L\left(u_{k+1},q_{k}\right)-L\left(u_{k},q_{k}\right)\\|_{X}^{2}+\\|C\left(u_{k+1}-u_{k}\right)\\|_{Y}^{2}+\\|L\left(u_{k+1},q_{k+1}\right)-L\left(u_{k+1},q_{k}\right)\\|_{X}^{2}\right)$
	$\displaystyle\leq J(u_{0},q_{0}),$		(3.9)

which implies that $\sum_{k=0}^{\infty}\|L\left(u_{k+1}-u_{k},q_{k}\right)\|_{X}^{2}$ is bounded. Then we can conclude using Assumption 2 that $\{u_{k}\}_{k=1}^{\infty}$ forms a Cauchy sequence and thus converges to some $u^{\star}\in{dom(L)}$ . Since $\sum_{k=0}^{m-1}E(q_{k},q_{k+1})$ is uniformly bounded for all $m$ , we can derive that $\{q_{k}\}_{k=1}^{\infty}$ converges to some $q^{\star}\in Z$ from the strict convexity of $\psi$ . As the sequence $\{J(u_{k},q_{k})\}_{k=1}^{\infty}$ is monotone decreasing, there exists a limit $J^{\star}$ . Following the similar argument in the proof of Theorem 1, we conclude that $J^{\star}=J(u^{\star},q^{\star})$ 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 $(u_{0},q_{0})$ , 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 $\sigma$ and $\mu$ simultaneously in a Lipschitz domain $\Omega\in\mathbb{R}^{d}$ ( $d=2,3$ ) from the model equation

-\nabla\cdot(\sigma(x)\nabla u)+\mu(x)\,u=0~{}~{}\mbox{in}~{}~{}\Omega

(4.1)

with a pair of Cauchy data $(f,h)$ on the boundary $\partial\Omega$ , i.e.,

u|_{\partial\Omega}=f,\quad\frac{\partial u}{\partial\nu}|_{\partial\Omega}=h.

Throughout this section, we shall use the notation $g:=\delta_{\partial\Omega}h$ in the total least-squares formulation, where $\delta_{\partial\Omega}$ denotes the Neumann to source map. To define the Neumann to source map, we first introduce the boundary restriction mapping $\gamma_{0}$ on $\mathscr{D}(\overline{\Omega})$ , i.e., $\gamma_{0}u$ denotes the boundary value of $u\in\mathscr{D}(\overline{\Omega})$ . Then we use $T$ to denote the trace operator [22], which is formally defined to be the unique linear continuous extension of $\gamma_{0}$ as an operator from $L^{2}(\Omega)$ onto $H^{-1/2}(\partial\Omega)$ . Using Riesz representation theorem, there exists a function in $L^{2}(\Omega)$ , denoted by $\delta_{\partial\Omega}h$ , such that for any $v\in L^{2}(\Omega)$ ,

(\delta_{\partial\Omega}h,v)_{L^{2}(\Omega)}=\langle h,Tv\rangle_{H^{1/2}(\partial\Omega),H^{-1/2}(\partial\Omega)},

which will be denoted by $g$ 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 $\psi$ 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 $L^{1}$ – $H^{1}$ regularization term $\phi$ for a coefficient $q:\Omega\rightarrow\mathbb{R}$ :

\phi(q;\alpha,\beta,q_{{}_{0}},{q_{{}_{1}}})=\int_{\Omega}\frac{\alpha}{2}(|\nabla q|^{2}+q^{2})dx+\int_{\Omega}\beta|q|\,dx+\chi(q;q_{{}_{0}},q_{{}_{1}}),

(4.2)

In this revision we change the semi-norm of $H^{1}$ into the full norm in the regularization to ensure the strict convexity of $\phi$ . where $\chi(q;q_{{}_{0}},q_{{}_{1}})$ is given by

\chi(q;q_{{}_{0}},q_{{}_{1}})=\left\{\begin{array}[]{cc}0&{q_{{}_{0}}}\leq q\leq q_{{}_{1}},\\ \\ \infty&\mbox{otherwise,}\end{array}\right.

and $q_{{}_{0}}$ and $q_{{}_{1}}$ are the lower and upper bounds of the coefficient $q$ respectively. The first term $\int_{\Omega}\frac{\alpha}{2}(|\nabla q|^{2}+q^{2})dx$ is the $H^{1}$ regularization term, the second term $\int_{\Omega}\beta|q|\,dx$ corresponds to the $L^{1}$ regularization, and the third term $\chi(q;q_{{}_{0}},q_{{}_{1}})$ enforces the reconstruction to meet the constrains of the coefficient.

In practice, the $L^{1}$ regularization enhances the reconstruction and helps find geometrically sharp inclusions, but might cause spiky results. The $H^{1}$ 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 $\alpha$ and $\beta$ 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 $\sigma$ and $\mu$ simultaneously. Similarly to the formulation (2.1), we now have $q=(\sigma,\mu)$ , and the operator $L$ is given by

L(u,\mathbf{p},\sigma,\mu)=\begin{pmatrix}-\nabla\cdot\mathbf{p}+\mu\,u\\ \mathbf{p}-\sigma\,\nabla u\end{pmatrix},

(4.3)

where we have introduced an auxiliary vector flux $\mathbf{p}$ and each entry of $L$ is of a first-order form such that two coefficients are separated naturally. Clearly, $L(v,q)$ is still bilinear with respect to the state variables $v=(u,\mathbf{p})$ and coefficients $q=(\sigma,\mu)$ . Using the first-order system, we can then come to the following total least-squares functional:

J(u,\mathbf{p},\sigma,\mu)=\|-\nabla\cdot\mathbf{p}+\mu(x)u-g\|_{L^{2}(\Omega)}^{2}+\|\mathbf{p}-\sigma(x)\,\nabla u\|_{(L^{2}(\Omega))^{d}}^{2}+\|Cu-f\|_{L^{2}(\partial\Omega)}^{2}+\psi_{1}(\sigma)+\psi_{2}(\mu),

(4.4)

where $C$ is the trace operator, and $\psi_{1}(\sigma)=\phi(\sigma;\alpha_{\sigma},\beta_{\sigma},\sigma_{{}_{0}},\sigma_{{}_{1}})$ and $\psi_{2}(\mu)=\phi(\mu;\alpha_{\mu},\beta_{\mu},\mu_{{}_{0}},\mu_{{}_{1}})$ are the corresponding mixed regularization terms of $\sigma$ and $\mu$ defined as in (4.2), $\mu_{{}_{0}}$ , $\mu_{{}_{1}}$ are lower and upper bounds of $\mu$ , and $\sigma_{{}_{0}}$ , $\sigma_{{}_{1}}$ are lower and upper bounds of $\sigma$ . We shall minimize (4.4) over $(u,\mathbf{p},\sigma,\mu)\in L^{2}(\Omega)\times(L^{2}(\Omega))^{d}\times L^{2}(\Omega)\times L^{2}(\Omega)$ , that is, $X=L^{2}(\Omega)\times(L^{2}(\Omega))^{d}$ , $Z=L^{2}(\Omega)\times L^{2}(\Omega)$ .

We will apply the ADI method to solve the least-squares formulation of $v=(u,\mathbf{p})$ and $q=(\sigma,\mu)$ :

\min_{(v,q)\in X\times Z}J(u,\mathbf{p},\sigma,\mu)=\|-\nabla\cdot\mathbf{p}+\mu(x)u-g\|_{L^{2}(\Omega)}^{2}+\|\mathbf{p}-\sigma(x)\,\nabla u\|_{(L^{2}(\Omega))^{2}}^{2}+\|Cu-f\|_{L^{2}(\partial\Omega)}^{2}+\psi_{1}(\sigma)+\psi_{2}(\mu).

(4.5)

Given an initial guess $(u_{0},\mathbf{p}_{0},\sigma_{0},\mu_{0})$ , we find a sequence $(u_{k},\mathbf{p}_{k},\sigma_{k},\mu_{k})$ for $k\geq 1$ as below:

•

Given $\sigma_{k}$ , $\mu_{k}$ , find $u=u_{k+1}$ , $\mathbf{p}=\mathbf{p}_{k+1}$ by solving

\min_{(u,\mathbf{p})\in X}J_{1}(u,\mathbf{p})=\|-\nabla\cdot\mathbf{p}+\mu_{k}(x)u-g\|_{L^{2}(\Omega)}^{2}+\|\mathbf{p}-\sigma_{k}(x)\,\nabla u\|_{(L^{2}(\Omega))^{d}}^{2}+\|Cu-f\|_{L^{2}(\partial\Omega)}^{2},

•

Given $u_{k+1}$ , $\mathbf{p}_{k+1}$ , find $\sigma=\sigma_{k+1}$ , $\mu=\mu_{k+1}$ by solving

\min_{(\sigma,\mu)\in Z}J_{2}(\sigma,\mu)=\|-\nabla\cdot\mathbf{p}_{k+1}+\mu(x)u_{k+1}-g\|_{L^{2}(\Omega)}^{2}+\|\mathbf{p}_{k+1}-\sigma(x)\,\nabla u_{k+1}\|_{(L^{2}(\Omega))^{d}}^{2}+\psi_{1}(\sigma)+\psi_{2}(\mu).

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 $\psi_{1}(\sigma)$ and $\psi_{2}(\mu)$ in (4.5) is convex and weakly lower semicontinuous. As the first term of (4.2) is the $H^{1}$ regularization term, $\psi_{1}(\sigma)$ and $\psi_{2}(\mu)$ are strictly convex by definition. One can also observe that there exists a positive constant $c$ such that

	$\displaystyle\psi_{1}(\sigma)=\int_{\Omega}\frac{\alpha}{2}(\|\nabla\sigma\|^{2}+\sigma^{2})dx+\int_{\Omega}\beta\|\sigma\|\,dx+\chi(\sigma;\sigma_{{}_{0}},\sigma_{{}_{1}})\geq c\\|\sigma\\|^{2}_{H^{1}(\Omega)},$		(4.6)
	$\displaystyle\psi_{2}(\mu)=\int_{\Omega}\frac{\alpha}{2}(\|\nabla\mu\|^{2}+\mu^{2})dx+\int_{\Omega}\beta\|\mu\|\,dx+\chi(\mu;\mu_{{}_{0}},\mu_{{}_{1}})\geq c\\|\mu\\|^{2}_{H^{1}(\Omega)},$		(4.6)

which imply that $\psi_{1}$ and $\psi_{2}$ are coercive in $H^{1}$ -norm.

Next we verify Assumption 2 on the closedness of $L_{q}$ and the coercivity of its graph norm. For fixed $q=(\sigma,\mu)$ , we denote the entries of $L_{q}$ by $L_{q,1}$ , $L_{q,2}$ , i.e., for $\delta v=(\delta u,\delta\mathbf{p})\in dom(L)$ ,

L_{q,1}\delta v=-\nabla\cdot\delta\mathbf{p}+\mu\delta u,\quad L_{q,2}\delta v=\delta\mathbf{p}-\sigma\nabla\delta u\,.

Applying Hölder’s inequality leads to

\displaystyle\|L_{q}\delta v\|_{X}\leq

\displaystyle\|\delta\mathbf{p}\|_{H(\text{div},\Omega)}+\|\delta u\|_{H^{1}(\Omega)}\|\mu\|_{L^{\infty}(\Omega)}+\|\delta u\|_{H^{1}(\Omega)}\|\sigma\|_{L^{\infty}(\Omega)}

(4.7)

for $q$ in the level set $\{q\in Z:\psi(q)=\psi_{1}(\sigma)+\psi_{2}(\mu)\leq c_{0}\}$ for some constant $c_{0}$ . Then we can deduce from (4.7) that $L_{q}$ is a closed linear operator. Next, we verify the coercivity of the graph norm. For simplicity, we consider the model problem with homogeneous Neumann boundary condition $h=0$ for the state variable, and set a side constraint for $dom(L)$ as $\mathbf{p}\cdot\nu=0$ on $\partial\Omega$ . Introduce the following notations

\sigma\nabla u-\mathbf{p}=\tilde{m},\quad-\nabla\cdot\mathbf{p}+\mu\,u=\tilde{g}.

(4.8)

From (4.8) and integration by part, we can derive

\int_{\Omega}\sigma\nabla u\cdot\nabla udx+\int_{\Omega}\mu u\cdot udx=\int_{\Omega}\nabla u\cdot\tilde{m}dx+\int_{\Omega}u\cdot\tilde{g}dx\,,

which implies

\|\nabla u\|_{(L^{2}(\Omega))^{d}}^{2}+\|u\|_{L^{2}(\Omega)}^{2}+\|\nabla\cdot\mathbf{p}\|_{L^{2}(\Omega)}^{2}+\|\mathbf{p}\|_{(L^{2}(\Omega))^{d}}^{2}\leq c\,(\|\tilde{m}\|_{(L^{2}(\Omega))^{d}}^{2}+\|\tilde{g}\|_{L^{2}(\Omega)}^{2})

(4.9)

for some constant $c>0$ when $q=(\sigma,\mu)$ satisfies $\psi_{1}(\sigma)\leq c_{0}$ and $\psi_{2}(\mu)\leq c_{0}$ for some constant $c_{0}>0$ . In this way we verify that the graph norm corresponding to $L_{q}$ defined as

|(u,\mathbf{p})|^{2}_{W,q}=\|\sigma\nabla u-\mathbf{p}\|_{(L^{2}(\Omega))^{d}}^{2}+\|-\nabla\cdot\mathbf{p}+\mu\,u\|_{L^{2}(\Omega)}^{2}

is uniformly coercive.

Then we consider Assumption 3 on the weakly sequentially closedness of operators $L$ and $C$ . It is noted that Assumption 3 is applied in Section 2 for general inverse medium problems to prove that the operator $L$ maps a subsequence of a bounded sequence $\{(u_{n},q_{n})\}_{n=1}^{\infty}$ to a converging sequence $\{L(u_{n},q_{n})\}_{n=1}^{\infty}$ with its limit equal to $L(u,q)$ , where $(u,q)$ is the limit of $\{(u_{n},q_{n})\}_{n=1}^{\infty}$ . In the analysis of the concrete DOT problem (4.1), the corresponding sequence $\{(u_{n},\mathbf{p}_{n},\sigma_{n},\mu_{n})\}_{n=1}^{\infty}$ is actually bounded in a stronger norm than $\|\cdot\|_{X}$ and $\|\cdot\|_{Z}$ , as shown in (4.6) and (4.9). As stated in Remark 1, to prove the well-poseness of the least-squares formulation (4.5), it suffices to verify that for a sequence $\{(u_{n},\mathbf{p}_{n},\sigma_{n},\mu_{n})\}_{n=1}^{\infty}$ bounded in $H^{1}(\Omega)\times H(\text{div},\Omega)\times H^{1}(\Omega)\times H^{1}(\Omega)$ , there exists a subsequence weakly converging to $(u,\mathbf{p},\sigma,\mu)$ , and operator $L$ defined in (4.3) satisfies $\{L(u_{n},\mathbf{p}_{n},\sigma_{n},\mu_{n})\}_{n=1}^{\infty}$ weakly converges to $L(u,\mathbf{p},\sigma,\mu)$ in $X$ .

The verification of the desired property of $L$ is as follows. Given the bounded sequence $\{(u_{n},\mathbf{p}_{n},\sigma_{n},\mu_{n})\}_{n=1}^{\infty}$ , there exists a subsequence, still denoted as $\{(u_{n},\mathbf{p}_{n},\sigma_{n},\mu_{n})\}_{n=1}^{\infty}$ , weakly converging to $(v,q):=(u,\mathbf{p},\sigma,\mu)$ in $H^{1}(\Omega)\times H(\text{div},\Omega)\times H^{1}(\Omega)\times H^{1}(\Omega)$ . Denoting $(u_{n},\mathbf{p}_{n})$ by $v_{n}$ and $(\sigma_{n},\mu_{n})$ by $q_{n}$ , there holds for any $V=(V_{1},V_{2})\in X$ that

\begin{split}&\,\quad\lim_{k\rightarrow\infty}(L(v_{n},q_{n}),V)_{X}-(L(v,q),V)_{X}\\ &=\lim_{k\rightarrow\infty}((L_{q_{n}}(v_{n}-v),V)_{X}+(L(v,q_{n}-q),V)_{X})\\ &=\lim_{k\rightarrow\infty}\big{(}(-\nabla\cdot(\mathbf{p}_{n}-\mathbf{p})+\mu_{n}(u_{n}-u),V_{1})+((\mathbf{p}_{n}-\mathbf{p})-\sigma_{n}\nabla(u_{n}-u),V_{2})\\ &\quad+((\sigma_{n}-\sigma)\nabla u,V_{1})+((\mu_{n}-\mu)u,V_{2})\big{)}\\ &=\,0,\end{split}

where we have used the weak convergence of $\{(u_{n},\mathbf{p}_{n},\sigma_{n},\mu_{n})\}_{n=1}^{\infty}$ and Hölder’s inequality. Therefore, $\{L(u_{n},\mathbf{p}_{n},\sigma_{n},\mu_{n})\}_{n=1}^{\infty}$ weakly converges to $L(u,\mathbf{p},\sigma,\mu)$ in $X$ . On the other hand, as $C$ is the trace operator on $\partial\Omega$ in DOT problem, it is linear and thus weakly sequentially closed, which completes our verification.

5 Numerical experiments

In this section, we carry out some numerical experiments on the DOT problem in different scenarios to illustrate the efficiency and robustness of the proposed two-stage algorithm in this work. Throughout these examples, we shall assume that we apply the Neumann boundary data $h$ on $\partial\Omega$ and measure the corresponding Dirichlet data $f$ to reconstruct the diffusion coefficient $\sigma$ and the absorption coefficient $\mu$ simultaneously. The basic algorithm involves two stages: we apply the direct sampling method (DSM) in the first stage to get some initial approximations of the two unknown coefficients, and then adopt the total least-squares method to achieve more accurate reconstructions of the coefficients.

5.1 Direct sampling method for initialization

For all the numerical experiments, we shall use the DSM in the first stage of our algorithm, in an attempt to effectively locate the multiple inclusions inside the computational domain with limited measurement data. Here we give a brief description of DSM and refer the readers to [19] for more technical details about this DSM that can identify multiple coefficients. The DSM develops two separate families of probing functions, i.e., the monopole and dipole probing functions, for constructing separate index functions for multiple physical coefficients. The inhomogeneities of coefficients can be approximated based on index functions due to the following two observations: the difference of scattered fields caused by inclusions can be approximated by the sum of Green’s functions of the homogeneous medium and their gradients; and the two sets of probing functions have the mutually almost orthogonality property, i.e., they only interact closely with the Green’s functions and their gradients respectively. Thus we can decouple the monopole and dipole effects and derive index functions for separate physical properties.

In practice, if the value of an index function $\phi$ for one physical coefficient at a sampling point $x$ is close to $1$ , the sampling point is likely to stay in the support of inhomogeneity; whereas if $\phi(x)$ is close to $0$ , the sampling point $x$ probably stays outside of the support. Hence the index functions give an image of the approximate support of the inhomogeneity, and we can determine a subdomain $D$ to locate the support from index functions. The subdomain $D$ could be chosen as $D=\{x\in\Omega:\phi(x)\geq\theta\}$ with $\theta$ being a suitable cut-off value, then we restrict the index function $\phi$ on $D$ . We adopt this choice in the numerical experiments to remove the spurious oscillations in the homogeneous background. Once we have the restricted index function $\phi|_{D}$ , we set the value of the approximation $\tilde{\phi}$ in $D$ as $\tilde{\phi}|_{D}=c_{\phi}\phi|_{D}$ , where the constant $c_{\phi}$ is some priori estimate of the true coefficient, and set the value of $\tilde{\phi}$ outside of $D$ as the background coefficient. In this way, we obtain $\tilde{\phi}$ as the initial guess for the total least-squares method for our further reconstruction in the second stage.

5.2 Examples

We present numerical results for some two-dimensional examples and showcase the proposed two-stage least-squares method for inverse medium problems using both exact and noisy data. Here the objective functional $J$ in (4.4) is discretized by a staggered finite difference scheme [23, 24]. The computational domain $\Omega=[0,1]\times[0,1]$ is divided into a uniform mesh consisting of small squares of width $h=0.01$ . The noisy measurements $u^{\delta}$ are generated point-wise by the formula

u^{\delta}(x)=u(x)+\epsilon\eta\max_{x\in\partial\Omega}|u(x)|

where $\epsilon$ is the relative noise level and $\eta$ follows the standard Gaussian distribution. The subdomain $D$ is chosen using the formula

D=\{x\in\Omega:\phi(x)\geq\theta\},

(5.1)

where $\phi$ is the index function from DSM and the cutoff value $\theta$ is taken in the range $(0.4,0.7)$ . As there exist limited theoretical results for the choice of regularization parameters for mixed regularization that we use, the regularization parameters $(\alpha_{\sigma},\beta_{\sigma},\alpha_{\mu},\beta_{\mu})$ in the functional for reconstructions are determined in a trial-and-error manner, and we present them for different examples in Table 1. The maximum number $K$ of alternating iterations is set at 50.

All the computations were performed in MATLAB (R2018B) on a desktop computer.

Example 1.

We consider the discontinuous diffusion and absorption coefficients $\sigma$ and $\mu$ with one inclusion each. The inclusion of $\sigma$ is of width $0.05$ and centered at $(0.25,0.65)$ as shown in Fig. 1(a) ; the inclusion of $\mu$ is of width $0.05$ and centered at $(0.35,0.3)$ as shown in Fig. 1(e). The magnitudes of the coefficients inside the inclusions are $20$ .

We shall use only one set of measurements to reconstruct the diffusion and absorption coefficients. It is shown in Fig. 1(b) and 1(f) that the index functions from the DSM separate inclusions of different physical nature well and give approximated locations, while the exact locations of small inclusions are difficult to detect. If we simply take the maximal points of the index functions in Fig. 1(b) and 1(f) as locations of the reconstructed inclusions, we may not be able to identify the true locations of inhomogeneity. Then we set the subdomain $D$ using information from the first stage by (5.1), see Fig. 1(c), 1(g), and set the value of approximation out of the subdomain as equal to the background coefficients. As in Fig. 1(d), 1(h), this example illustrates that the least-squares formulation in the second stage works very well to improve the reconstruction and provides a much more accurate location than that provided by DSM. With $10\%$ noise in the measurements, the reconstructions remain accurate as in Fig. 2, which shows that the two-stage algorithm gives quite stable final reconstructions with respect to the noise even though the DSM reconstructions become more blurry in Figs. 2(b) and 2(f). The regularization parameters $(\alpha_{\sigma},\beta_{\sigma},\alpha_{\mu},\beta_{\mu})$ are presented in Table 1.

Compared with the reconstructions derived from the DSM, the improvement of the approximations for both $\sigma$ and $\mu$ is significant: the recovered background is now mostly homogeneous, and the magnitude and size of the inhomogeneity approximate those of the true coefficients well. These results indicate clearly the significant potential of the proposed least-squares formulation with mixed regularization for inverse medium problems.

From Table 1 we obtain an insightful observation about the mixed regularization: the magnitude of parameter $\beta$ is much larger than that of $\alpha$ . We can conclude that the $L^{1}$ penalty plays a predominant role in improving the performance of reconstruction for such inhomogeneous coefficients, whereas the $H^{1}$ penalty yields a locally smooth structure.

Noise level	$\epsilon=0$		$\epsilon>0$
Example	$(\alpha_{\sigma},\beta_{\sigma})$	$(\alpha_{\mu},\beta_{\mu})$	$(\alpha_{\sigma},\beta_{\sigma})$	$(\alpha_{\mu},\beta_{\mu})$
1	(1.0e-2, 2.0e-2 )	(5.0e-4, 5.0e-4)	(1.0e-2, 2.0e-2)	(5.0e-4, 1.0e-3)
2.1	(1.0e-3, 5.0e-3)	N/A	(1.0e-3, 1.0e-2)	N/A
2.2	(1.0e-6, 1.0e-3)	N/A	(1.0e-6, 2.0e-3)	N/A
3	(1.0e-3, 1.0e-2)	(1.0e-2, 5.0e-3)	(1.0e-3, 2.0e-2)	(1.0e-2, 5.0e-3)
4	(1.0e-5, 5.0e-4)	N/A	(1.0e-5, 1.0e-3)	N/A

Table 1: The regularization parameters

(\alpha,\beta)

in each example for

\sigma

and

\mu

without noise and with noise.

Example 2.

We implement the algorithm to reconstruct the discontinuous diffusion coefficient $\sigma$ with two inclusions. We assume that $\mu$ is known and equal to the background coefficient $\mu_{0}$ in this example. One set of data is measured to locate two inclusions of $\sigma$ , which are centered at $(0.15,0.5)$ and $(0.5,0.85)$ respectively and of width $0.05$ . $\sigma$ is taken to be $20$ inside the regions as shown in Fig. 3(a).

The reconstructions for $\sigma$ using exact measurements and measurements with $20\%$ noise are shown in Fig. 3. In the first stage, the DSM gives the index function that separates these two inclusions well using only one set of data as in Fig. 3(b), while some inhomogeneity is observed in the background. This phenomenon comes from the ill-posedness of inverse medium problems and also the oscillation of fundamental solutions used in the DSM. Even so, the approximations in Fig. 3(b) and 3(f) still provide the basic modes of the inhomogeneity. We can identify that there are two inclusions, and capture the subdomain $D$ for the second stage in Figs. 3(c) and 3(g). The least-squares formulation with mixed regularization significantly improves the reconstruction: the locations of both inclusions are captured better with clear background and accurate size in Figs. 3(d), 3(h). Comparing with Figs. 3(d) and 3(h), one can observe that the reconstruction deteriorates only slightly in that the left inclusion shrinks a little bit when the noise level $\epsilon$ increases from 0 to $20\%$ . This example verifies that the two-stage algorithm is very robust with respect to the data noise.

Example 3.

In this example we reconstruct two inclusions of diffusion coefficient $\sigma$ that stay very close to each other. The two inclusions of diffusion coefficient $\sigma$ are centered at $(0.45,0.425)$ and $(0.55,0.575)$ respectively and of width $0.1$ as shown in Fig. 4(a). The coefficient $\sigma$ is $20$ in both regions.

The two inclusions in Example 2 are relatively far from each other, while in Example 3 we consider the case when two inclusions are quite close to each other, which is more challenging as it would be difficult to distinguish these two separated inclusions and reconstruct their locations and magnitudes precisely. As shown in Fig. 4(b), the index function from DSM presents limited information on the diffusion coefficient with only one set of data, and shows one connected inclusion. With $2\%$ noise, the index function is blurred a lot as shown in Fig. 4(f). In both noisy and noiseless cases, only one subdomain can be detected from this index function from the first stage. The second stage still presents well separated reconstructions with the size and magnitude that match the exact diffusion coefficient well as shown in Fig. 4(d). For the case with $2\%$ noise, this two-stage algorithm also gives satisfying reconstruction in Fig. 4(h). Compared with Fig. 4(d), it is observable that the left inclusion moves towards x-axis and elongates a little bit, while we can still tell the sizes and locations of two inclusions from the reconstruction. This shows that the least-squares method provides much more details than DSM and is relatively robust with respect the the noise in the measurement.

Note that for the noise level higher than $2\%$ , the DSM can not provide a feasible initial guess for the least-squares method in the second stage with only one set of data, but with an initial guess that reflects the mode of the true coefficient, the least-squares method in the second stage has great tolerance of noise as shown in other examples.

Example 4.

With this example, we reconstruct $\sigma$ and $\mu$ simultaneously, with two inclusions for each coefficient. The inclusions are in the following scenario: The inclusions of diffusion coefficient $\sigma$ are of width $0.1$ and centered at $(0.5,0.25)$ , $(0.5,0.75)$ respectively, and the magnitude inside the region is $20$ ; The inclusions of absorption coefficient $\mu$ are of width $0.1$ and centered at $(0.25,0.5)$ , $(0.75,0.5)$ respectively, and the magnitude inside the region is $20$ , as shown in Figs. 5(a) and 5(e).

In Example 1, we have one inclusion for coefficients $\mu$ and $\sigma$ each, and it is shown that this algorithm can reconstruct both medium coefficients well. Example 4 is more challenging than Example 1, as the existence of two inclusions for each coefficient will influence the reconstruction of the other coefficient. As shown in Figs. 5(b) and 5(f), the index functions separate these inclusions well and give a rough approximation for their locations with only one set of data. However, it can be observed that the maximal points of index functions actually differ from the exact coefficients, and one can not capture the information of inclusions by simply taking square of the index functions. In the second stage of the proposed algorithm, the locations are modified as shown in Figs. 5(d) and 5(h). When we consider the case with $20\%$ noise, DSM provides blurry approximations as shown in Figs. 6(b) and 6(f). In the second stage of reconstruction, Figs. 6(d) and 6(h) present results almost the same as Figs. 5(d) and 5(h), which once again prove the robust performance of the least-squares method.

Example 5.

In this example we reconstruct diffusion coefficient $\sigma$ with ring-shaped square inclusion as shown in Fig. 7(a) with two sets of measurement. The outer and inner side length of the ring are $0.2$ and $0.15$ , and the rectangle ring is centered at $(0.5,0.6)$ . The coefficient $\sigma$ is taken to be $20$ inside the region.

We use two sets of measurement from different directions for the two-stage least-squares method. For this challenging case, the index function from the DSM only reflects an approximated location of inclusion as in Fig. 7(b), and we have no clue about the shape of inclusion from only DSM reconstruction. With two sets of measurement, the least-squares formulation can reconstruct the edges of the ring-shape inclusion as shown in Fig. 7(d). After adding $20\%$ noise, one has the reconstruction (see Fig. 7(h)) that is very similar to the result without noise, which shows the approximation is quite stable with respect to the noise.

References

[1] S. R. Arridge, Optical tomography in medical imaging, Inverse Problems, 15(2):R41, 1999.
[2] A. P. Gibson, J. C. Hebden, and S. R. Arridge, Recent advances in diffuse optical imaging, Phys. Med. Biol., 50(4):R1, 2005.
[3] O. Dorn, H. Bertete-Aguirre, J. G. Berryman, and G. C. Papanicolaou, A nonlinear inversion method for 3d electromagnetic imaging using adjoint fields, Inverse Problems, 15(6):1523, 1999.
[4] G. Bao and P. Li, Inverse medium scattering problems for electromagnetic waves, SIAM J. Appl. Math., 65(6):2049–2066, 2005.
[5] M. Cheney, D. Isaacson, and J. C. Newell, Electrical impedance tomography, SIAM Rev., 41(1):85–101, 1999.
[6] K. Ito and B. Jin, Inverse problems: Tikhonov theory and algorithms, World Scientific, Hackensack, NJ, 2014.
[7] R. V. Kohn and M. Vogelius, Relaxation of a variational method for impedance computed tomography, Commun. Pure Appl. Math., 40(6):745–777, 1987.
[8] A. Wexler, B. Fry, and M. R. Neuman, Impedance-computed tomography algorithm and system, Appl. Opt., 24(23):3985–3992, 1985.
[9] E. Chung, K. Ito and M. Yamamoto, Least square formulation for ill-posed inverse problems and applications, Appl. Anal., 1–15, 2021.
[10] R. Lattes, J.L. Lions, The method of quasi-reversibility. applications to partial differential equations, Elsevier, New York, 1969.
[11] G. Tusnady, I. Csiszar, Information geometry and alternating minimization procedures, Statistics and Decisions, Supp.1:205–237, 1984.
[12] C. L. Byrne, Iterative optimization in inverse problems, CRC Press, 2014.
[13] D. C. Jespersen, A least squares decomposition method for solving elliptic equations, Math. Comput., 31(140):873–880, 1977.
[14] P. B. Bochev and M. D. Gunzburger, Finite element methods of least-squares type, SIAM Rev., 40(4):789–837, 1998.
[15] A. I. Pehlivanov, G. F. Carey, and R. D. Lazarov, Least-squares mixed finite elements for second-order elliptic problems, SIAM J. Numer. Anal., 31(5):1368–1377, 1994.
[16] P. B. Bochev and M. D. Gunzburger, Least-squares finite element methods, volume 166, Springer Science & Business Media, 2009.
[17] K. Ito, B. Jin, and J. Zou, A two-stage method for inverse acoustic medium scattering, J. Comput. Phys., 2013, 237: 211–223, 2012.
[18] Y. T. Chow, K. Ito, and J. Zou, A direct sampling method for electrical impedance tomography, Inverse Problems, 30(9):095003, 2014.
[19] Y. T. Chow, F. Han, and J. Zou, A direct sampling method for simultaneously recovering inhomogeneous inclusions of different nature, arXiv preprint arXiv:2005.05499, 2020.
[20] B. Harrach, On uniqueness in diffuse optical tomography, Inverse Problems, 25(5):055010, 2009.
[21] L. M. Bregman, The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming, USSR Comput. Math. and Math. Phys., 7(3):200–217, 1967.
[22] V. Girault and P. A. Raviart, Finite element methods for Navier-Stokes equations: theory and algorithms, volume 5, Springer Science & Business Media, 2012.
[23] J. Virieux, Sh-wave propagation in heterogeneous media: Velocity-stress finite-difference method, Geophysics, 49(11):1933–1942, 1984.
[24] J. O. A. Robertsson, J. O. Blanch, and W. W. Symes, Viscoelastic finite-difference modeling, Geophysics, 59(9):1444–1456, 1994.
[25] M. Burger and S. Osher, Convergence rates of convex variational regularization, Inverse Problems, 20(5):1411, 2004.
[26] Y. T. Chow, K. Ito, K. Liu, and J. Zou, Direct sampling method for diffusive optical tomography, SIAM J. Sci. Comput., 37(4):1658-A1684, 2015.
[27] Y. T. Chow, K. Ito, and J. Zou, A time-dependent direct sampling method for recovering moving potentials in a heat equation, SIAM J. Sci. Comput., 40(4):2720-A2748, 2018.

$\displaystyle 0$	$\displaystyle=\big{(}L_{q_{k}}^{}(L_{q_{k}}u_{k+1}-g),u_{k+1}-u_{k}\big{)}_{X}+\big{(}C^{}(Cu_{k+1}-f),u_{k+1}-u_{k}\big{)}_{Y}$
	$\displaystyle=\frac{1}{2}\Big{(}\\|L_{q_{k}}u_{k+1}-g\\|_{X}^{2}+\\|Cu_{k+1}-f\\|_{Y}^{2}-\left(\\|L_{q_{k}}u_{k}-g\\|_{X}^{2}+\\|Cu_{k}-f\\|_{Y}^{2}\right)$
	$\displaystyle\quad+\\|L_{q_{k}}(u_{k+1}-u_{k})\\|_{X}^{2}+\\|C(u_{k+1}-u_{k})\\|_{Y}^{2}\Big{)}\,.$	(3.5)