Nonlinear Iterative Hard Thresholding for Inverse Scattering

TL;DR

This paper introduces a nonlinear iterative hard thresholding algorithm for reconstructing sparse scatterers in inverse scattering problems, with analysis of convergence and error, supported by numerical simulations.

Contribution

It presents a novel nonlinear generalization of iterative hard thresholding tailored for inverse scattering, with theoretical convergence analysis.

Findings

Convergence of the proposed method is established through coherence estimates.

Numerical simulations validate the effectiveness of the algorithm.

The method outperforms traditional approaches in sparse inverse scattering reconstruction.

Abstract

We consider the inverse scattering problem for sparse scatterers. An image reconstruction algorithm is proposed that is based on a nonlinear generalization of iterative hard thresholding. The convergence and error of the method was analyzed by means of coherence estimates and compared to numerical simulations.

Figures

Tables

Table 1. Table 1 . The central slice of the reconstructions of model 1 for varying susceptibilities. A comparison of the linear IHT and the fully nonlinear IHT algorithms is shown. The last column plots the error Y err subscript 𝑌 err Y_{\rm err} . The field of view in each image is 5 ​ λ × 5 ​ λ 5 𝜆 5 𝜆 5\lambda\times 5\lambda .

potential	Linear IHT	Nonlinear IHT	Error
${\eta }_0=0.01$
${\eta }_0=0.06$
${\eta }_0=0.1$
${\eta }_0=0.4$

Table 2. Table 2 . The central slice of the reconstructions for model 1 with η 0 subscript 𝜂 0 \eta_{0} =0.06 and 0.1. The reconstructions were conducted with the linear, second and third Born approximations. The Born series converges for the top row and diverges for the bottom row. Note that in both cases, the fully nonlinear IHT algorithm converges. The field of view in each image is 5 ​ λ × 5 ​ λ 5 𝜆 5 𝜆 5\lambda\times 5\lambda .

	$M=1$	$M=2$	$M=3$
${\eta }_0=0.06$
${\eta }_0=0.1$

Equations

$$\mathbf{y}=\Phi \mathbf{x}+\mathbf{ϵ},$$

$${\mathbf{x}}_{n+1}=H_s\left({\mathbf{x}}_n+{\Phi }^{\ast }(\mathbf{y}-\Phi {\mathbf{x}}_n)\right),$$

$$\mathbf{y}=\Phi (\mathbf{x})+\mathbf{ϵ},$$

$${\mathbf{x}}_{n+1}=H_s\left({\mathbf{x}}_n+{\Phi }_{{\mathbf{x}}_n}^{\ast }(\mathbf{y}-{\Phi }_{{\mathbf{x}}_n}{\mathbf{x}}_n)\right),$$

$$(1-{\delta }_s)\Vert \mathbf{x}{\Vert }_2^2\le \Vert \Phi (\mathbf{x}){\Vert }_2^2\le (1+{\delta }_s)\Vert \mathbf{x}{\Vert }_2^2.$$

$$\mu (\Phi )=\underset{j\ne k}{\max }\frac{|⟨{\Phi }_j,{\Phi }_k⟩|}{\Vert {\Phi }_j{\Vert }_2\Vert {\Phi }_k{\Vert }_2}.$$

$${\mathbf{x}}_{n+1}=H_s\left({\mathbf{x}}_n+{\Phi }_{{\mathbf{x}}_n}^{\ast }(\mathbf{y}-{\Phi }_{{\mathbf{x}}_n}{\mathbf{x}}_n)\right),$$

$$\Vert {\mathbf{x}}_n-\mathbf{x}{\Vert }_1\le \rho \Vert {\mathbf{x}}_{n-1}-\mathbf{x}{\Vert }_1+(3s+1)\Vert {\Phi }_{{\mathbf{x}}_n}^{\ast }{\mathbf{e}}_n{\Vert }_{\infty },$$

$$\Vert {\mathbf{x}}_n-\mathbf{x}{\Vert }_1\le {\rho }^n\Vert {\mathbf{x}}_0-\mathbf{x}{\Vert }_1+\sum _{j=0}^n{\rho }^j(3s+1)\Vert {\Phi }_{{\mathbf{x}}_j}^{\ast }{\mathbf{e}}_j{\Vert }_{\infty }.$$

$$\Vert {\mathbf{x}}_n-\mathbf{x}{\Vert }_1\le {\rho }^n\Vert {\mathbf{x}}_0-\mathbf{x}{\Vert }_1+\frac{M(3s+1)}{1-\rho }$$

$${\nabla }^{2}u(\mathbf{x})+k^2(1+\eta (\mathbf{x}))u(\mathbf{x})=-S(\mathbf{x}),$$

$${\nabla }^2{u}_i(\mathbf{x})+k^2{u}_i(\mathbf{x})=-S(\mathbf{x}).$$

$${\nabla }^2{u}_s(\mathbf{x})+k^2{u}_s(\mathbf{x})=-\eta (\mathbf{x})u(\mathbf{x}),$$

$$\underset{r\to \infty }{\lim }r\left(\frac{\partial u_s}{\partial r}-ik{u}_s\right)=0,$$

$${\nabla }^{2}G(\mathbf{x},\mathbf{y})+k^{2}G(\mathbf{x},\mathbf{y})=-\delta (\mathbf{x}-\mathbf{y}),$$

$$G(\mathbf{x},\mathbf{y})=\frac{e^{ik|\mathbf{x}-\mathbf{y}|}}{4\pi |\mathbf{x}-\mathbf{y}|}.$$

$$u_s(\mathbf{x})=k^2{\int }_{\Omega }G(\mathbf{x},\mathbf{y})\eta (\mathbf{y})u(\mathbf{y})d\mathbf{y}.$$

$$u_s(\mathbf{x})=$$

$$+k^4{\int }_{\Omega }{\int }_{\Omega }G(\mathbf{x},{\mathbf{y}}_1)\eta ({\mathbf{y}}_1)G({\mathbf{y}}_1,{\mathbf{y}}_2)\eta ({\mathbf{y}}_2)u_i({\mathbf{y}}_2)d{\mathbf{y}}_{1}d{\mathbf{y}}_2+\cdots ,$$

$$u_s(\mathbf{x})\approx k^2{\int }_{\Omega }G(\mathbf{x},\mathbf{y})\eta (\mathbf{y})u_i(\mathbf{y})d\mathbf{y}.$$

$$u_s=$$

$$u_s=\tilde{\mathcal{G}}(1-\mathcal{V}\mathcal{G}{)}^{-1}\mathcal{V}{u}_i.$$

$$\mathcal{T}=\mathcal{V}(1-\mathcal{G}\mathcal{V}{)}^{-1}=(1-\mathcal{V}\mathcal{G}{)}^{-1}\mathcal{V},$$

$$u_s=\tilde{\mathcal{G}}\mathcal{T}{u}_i.$$

$$u_i(\mathbf{x})=e^{ik\hat{\mathbf{d}}\cdot \mathbf{x}},$$

$$\frac{e^{ik|\mathbf{x}-\mathbf{y}|}}{4\pi |\mathbf{x}-\mathbf{y}|}\sim \frac{e^{ikR}}{4\pi R}{e}^{-ik\hat{\mathbf{x}}\cdot \mathbf{y}}.$$

$$u_s(\mathbf{x})\sim A(\hat{\mathbf{x}})\frac{e^{ikR}}{4\pi R},$$

$$A(\hat{\mathbf{x}})=k^2{\int }_{\Omega }{e}^{-ik\hat{\mathbf{x}}\cdot \mathbf{y}}\eta (\mathbf{y})u(\mathbf{y})d\mathbf{y}.$$

$$Y=A(I-V\Gamma {)}^{-1}VB.$$

$$A_{mn}=e^{-ik{\hat{\mathbf{x}}}_m\cdot {\mathbf{r}}_m},$$

$$B_{mn}=e^{ik{\hat{\mathbf{d}}}_n\cdot {\mathbf{r}}_m},$$

Full text

1. Introduction

1.1. Background

Scattering experiments are of fundamental importance in nearly every branch of physics. They also arise in numerous applied fields including optics, seismology and biomedical imaging. In this paper, we consider the following experimental setup. A wave field is incident on a medium and the resulting scattered field is measured. The spatially-dependent properties of the medium are encoded in coefficients that arise in the wave equation. Depending on the physical setting, the associated inverse scattering problem is to reconstruct one or more coefficients from boundary measurements. Much is known about theoretical aspects of the problem, especially on matters of uniqueness, stability and reconstruction [8, 9, 11, 12, 13, 17, 22]. By uniqueness we mean the injectivity of the forward map from the coefficients, sometimes referred to as the scattering potential, to the scattered field. Stability refers to continuity of the inverse map from the scattered field to the potential. We note that inverse scattering problems are typically ill-posed, which means that the inverse map must be suitably regularized to achieve stable inversion.

The forward map depends nonlinearly on the scattering potential. Linearization of the forward map is referred to as the first Born approximation (FBA). The ISP within the FBA is relatively well understood. The development of reconstruction methods for the nonlinear ISP is an area of active research. These include optimization, qualitative and direct methods. There is also considerable interest in developing reconstruction algorithms for scatterers that are known a priori to be sparse [6, 16, 18, 30]. For instance, similar algorithms to iterative hard thresholding (IHT) have been applied to the ISP in the FBA [14, 15]. We note that IHT has also proven to be effective in many applications other than inverse scattering [3].

In this paper we apply a generalization of IHT to the nonlinear ISP. Our results provide nonlinear corrections to inversion within the FBA. In particular, we characterize the convergence and error of the method. The analysis is primarily carried out within the framework of coherence estimates. Extensions to restricted isometry property (RIP) estimates are also provided. Numerical simulations are used to illustrate the results.

1.2. Related work

A nonlinear version of IHT was proposed in [2]. An analysis of this algorithm was performed by making use of the restricted isometry property (RIP). However, many physical measurement systems of interest do not obey even moderately strong RIP bounds. Moreover, computationally efficient methods to compute RIP constants are not available, thereby preventing the practical verification of RIP bounds. Thus, sparse imaging in deterministic settings is usually characterized by the mutual coherence. We note that coherence calculations for various scattering experiments have been reported [15, 23].

1.3. Organization

The paper is organized as follows. In section 2.1, we state the nonlinear IHT algorithm and provide a convergence analysis based on coherence estimates. Relevant background on inverse scattering is then introduced in  section 2.2. In  section 3 we apply the nonlinear IHT algorithm to inverse scattering. An analysis of the resulting algorithm is presented in section 4. Our results are illustrated with numerical simulations in section 5. A related analysis of nonlinear IHT for inverse scattering using RIP is given in Appendix B.

2. Preliminaries and Background

Throughout the paper, matrices will be denoted by uppercase letters and vectors will be denoted by bold lowercase letters. For a matrix $A$ with entries $A_{ij}$, we denote its $j$th column by $A_{j}$. We denote the adjoint of $A$ by $A^{*}$. For an $N\times 1$ vector $\mathbf{x}$, $D_{\mathbf{x}}$ is the $N\times N$ diagonal matrix with diagonal entries given by $\mathbf{x}$ and $\mathbf{\hat{x}}$ is the vector of unit length in the same direction as $\mathbf{x}$. We define the max norm of a matrix $A$ as $\|A\|_{\max}=\max_{ij}|A_{ij}|$. We also let $\Omega\subset\mathbb{R}^{3}$ denote a bounded domain with smooth boundary $\partial\Omega$.

2.1. Iterative Hard Thresholding

Consider the linear system

[TABLE]

where $\Phi$ is an $M\times N$ matrix, $\mathbf{x}$ is an $s$-sparse vector of length $N$ with $s$ non-zero entries, and $\bm{\epsilon}$ is a vector of length $M$. We will refer to $\Phi$ as the sensing matrix and $\bm{\epsilon}$ accounts for the effects of noise. When $\Phi$ has fewer rows than columns, the linear system is under-determined and the problem of finding or recovering $\mathbf{x}$ is frequently referred to as compressive sensing. Here the matrix $\Phi$ senses $\mathbf{x}$, albeit in a noisy fashion, and the vector $\mathbf{y}$ is a set of noisy observations of $\mathbf{x}$. Despite this problem being under-determined, if there are enough rows in $\Phi$ or observations of $\mathbf{x}$ (roughly $O(s\log N)$), then $\mathbf{x}$ can be recovered from its observations $\mathbf{y}$, provided $\mathbf{x}$ is $s$-sparse, up to the level of the noise $\bm{\epsilon}$.

IHT is an algorithm for sparse signal recovery, first proposed by Blumensath and Davies in [3], and subsequently analyzed and developed in [1, 4, 5]. The IHT algorithm recovers $\mathbf{x}$ by the iterative process

[TABLE]

where $H_{s}(x)$ is the nonlinear thresholding operator that sets all but the largest (in absolute value) $s$ elements of $\mathbf{x}$ to zero. This algorithm is the classical Richardson first-order iteration with the additional application of a thresholding operator to promote sparsity [19].

Since not all compressive sensing problems are linear, it is of interest to consider the nonlinear problem

[TABLE]

where $\Phi$ is a nonlinear function rather than a linear one. Blumensath [2] proposed solving the above problem with an analogous nonlinear IHT algorithm. This algorithm replaces the iterative step in Eq. (2) by

[TABLE]

where $\Phi_{\mathbf{x}_{n}}$ is the linearization of $\Phi$ at the point $\mathbf{x}_{n}$. Blumensath analyzed this algorithm by making use of the restricted isometry property (RIP), which gives a measure of how “close” to orthogonal the sensing matrix is when applied to $s$-sparse vectors. In particular, a matrix $\Phi$ satisfies the RIP of order $s$ if there exists a constant $\delta_{s}$ with $0<\delta_{s}<1$ such that, for all $s$-sparse vectors $\mathbf{x}$,

[TABLE]

See [10] for further details. Rudelson and Vershynin [28] have shown that there are several randomized algorithms for generating matrices that satisfy the RIP, namely a random matrix with iid sub-Gaussian entries or a random submatrix of a matrix with bounded entries with a sufficient number of rows (drawn uniformly at random). There are also deterministic constructions of matrices that satisfy the RIP that have nearly but not quite optimal dimensions.

While the analysis of both linear and nonlinear compressive sensing algorithms using RIP yields strong theoretical results, most scattering experiments with deterministic source and detector geometries do not yield sensing matrices for which RIP estimates are easy to verify either analytically or computationally. Therefore, we utilize a weaker notion than RIP, namely coherence, to analyze the nonlinear IHT algorithm. The coherence of the sensing matrix $\Phi$ is defined by

[TABLE]

Equivalently, if $\Phi$ is column-normalized ($\|\Phi_{j}\|_{2}=1$ for all $j$), then the mutual coherence is the largest (in absolute value) off-diagonal entry in the matrix $\Phi^{*}\Phi$. This quantity is considerably easier to compute, does not depend on the putative sparsity of the underlying data, and does not require sophisticated algebraic constructions (for either random or deterministic matrices), in contrast to the RIP constructions. There are, however, deterministic constructions of RIP matrices based upon simple coherence estimates but these have far from optimal dimensions [26].

There are convergence results for linear IHT expressed in terms of the coherence bound of the sensing matrix and the sparsity level of the unknown vector $\mathbf{x}$ in Eq. (1) (see [29]). Our analysis, using coherence rather than RIP of the sensing matrix, of the nonlinear IHT algorithm shows that it behaves similarly to the linear IHT algorithm under two conditions: (1) the linearization $\Phi_{\mathbf{x}_{n}}$ obeys a coherence bound for all iterations and (2) the error in the linearization is sufficiently small.

We now state our main result on the convergence of nonlinear IHT using coherence rather than RIP. The proof is given in Appendix A.

Theorem 2.1**.**

Let $\mathbf{y}=\Phi(\mathbf{x})+\bm{\epsilon}$ and let $\{\mathbf{x}_{n}\}$ be the sequence generated by the iteration

[TABLE]

where $\Phi_{\mathbf{x}_{n}}$ is the linearization of $\Phi$ at $\mathbf{x}_{n}$. Suppose that $\Phi_{\mathbf{x}_{n}}$ obeys the coherence estimate $\mu(\Phi_{\mathbf{x}_{n}})\leq\mu_{0}$ for all $n\geq 1$. If $\mathbf{x}$ is $s$-sparse and $\mu_{0}<\frac{1}{3s+1}$, then

[TABLE]

where $\mathbf{e}_{n}=\Phi(\mathbf{x})-\Phi_{\mathbf{x}_{n}}\mathbf{x}+\bm{\epsilon}$ and $\rho=\mu_{0}(3s+1)$.

We note that the above result reduces to the known linear IHT convergence result when $\Phi$ is linear [29]. Iterating (8), we obtain for all $n\geq 1$,

[TABLE]

Eq. (9) illustrates the importance of the error term in Theorem 2.1. Since by hypothesis $\rho<1$, the condition that $\|\Phi^{*}_{x_{j}}e_{j}\|_{\infty}$ is bounded for all $j$ leads to a convergent geometric series in (9), which is summarized in the following corollary.

Corollary 2.2**.**

Let $\mathbf{y}=\Phi(\mathbf{x})+\bm{\epsilon}$ and $\{\mathbf{x}_{n}\}$ be the sequence generated by (7). Suppose that $\mathbf{x}$ is $s$-sparse and $\rho=\mu_{0}(3s+1)<1$. If there exists a constant $M$ such that for all $j$, $\|\Phi^{*}_{\mathbf{x}_{j}}\mathbf{e}_{j}\|_{\infty}<M$, then the error estimate

[TABLE]

holds for all $n\geq 1$.

The above results are applied in Section 4.2 to the inverse scattering problem (ISP), which we describe next.

2.2. Inverse scattering problem

Here we describe the formulation of the ISP that we require in this paper. For simplicity, we restrict our attention to scalar waves. Our results extend naturally to other types of scattering problems, though some modifications are required [7].

We begin with the forward scattering problem. We consider the scattering of a time-harmonic scalar wave $u$ from a medium with scattering potential $\eta$, where we assume that $\eta$ is compactly supported in a bounded domain $\Omega$. The field satisfies the Helmholtz equation

[TABLE]

where $k$ is the wavenumber and $S(\mathbf{x})$ is the source term. The total field $u$ can be decomposed into the incident field $u_{i}$ and scattered field $u_{s}$ by $u=u_{i}+u_{s}$. The incident field obeys Eq. (11) in the absence of scattering:

[TABLE]

It follows that the scattered field obeys

[TABLE]

which we supplement with the Sommerfeld radiation condition

[TABLE]

where $r=|\mathbf{x}|$. Eq. (11) can be converted to an equivalent integral equation by introducing a Green’s function. The Green’s function $G$ satisfies

[TABLE]

together with the Sommerfeld radiation condition. In three dimensions, $G$ is given by

[TABLE]

It follows that $u_{s}$ obeys the integral equation

[TABLE]

By iterating the above, we obtain the Born series

[TABLE]

which yields an explicit formula for the scattered field. Convergence of the Born series for scalar waves is considered in [21]. In the case of weak scattering, the scattered field can be well approximated by the first term of the Born series, which yields the first Born approximation

[TABLE]

Multiple scattering is accounted for by including the higher order order terms in the Born series.

It will prove useful to express the Born series in operator notation. We define $\mathcal{G}$ to be the integral operator with kernel $G(\mathbf{x},\mathbf{y})$, where both $\mathbf{x}$ and $\mathbf{y}$ belong to $\Omega$, and define $\tilde{\mathcal{G}}$ to be the integral operator with the same kernel, but with $\mathbf{x}$ outside of $\Omega$. We also let $\mathcal{V}$ be the operator that corresponds to multipication by $k^{2}\eta$. The Born series (2.2) thus becomes

[TABLE]

By summing the geometric series we obtain

[TABLE]

Introducing the $T$-matrix operator which is defined by

[TABLE]

we find that

[TABLE]

The above formulation is quite general and holds for many experimental configurations. For the remainder of the paper, we consider a specific scattering experiment in which the incident field is a plane wave, and measurements of the scattered field are taken in the far field on a ball of radius $R$ with $kR\gg 1$, as shown in Fig. 1. That is, the incident field is given by

[TABLE]

where the unit vector $\mathbf{\hat{d}}$ is the direction of propagation. In the far field, the Green’s function is asymptotically of the form

[TABLE]

Eq. (17) then becomes

[TABLE]

where the scattering amplitude $A$ is defined by

[TABLE]

Thus the scattered field behaves as an outgoing spherical wave with amplitude $A$.

The goal of the ISP is to reconstruct $\eta$ from measurements of $A$. To proceed, we suppose that $N_{d}$ measurements of the scattering amplitude in the directions $\mathbf{\hat{x}}_{1},\dots,\mathbf{\hat{x}}_{N_{d}}$ are acquired for each of $N_{s}$ incident plane waves in the directions $\mathbf{\hat{x}}_{1},\dots,\mathbf{\hat{x}}_{N_{s}}$. We assemble these measurements in a $N_{d}\times N_{s}$ matrix that we denote by $Y$. We also discretize $\Omega$ into $N$ voxels of volume $h^{3}$ with centers $\mathbf{r}_{1},\dots,\mathbf{r}_{N}$. The discretized version of (21) is then of the form

[TABLE]

Here $A$ is an $N_{d}\times N$ matrix, $B$ is an $N\times N_{s}$ matrix, $\Gamma$ is an $N\times N$ matrix, and $V$ is a $N\times N$ diagonal matrix. These matrices are defined as

[TABLE]

It is important to note that the dependence of $Y$ on $V$ is nonlinear. This can be seen by expanding the matrix inverse in (28) to obtain

[TABLE]

which converges when $\|V\Gamma\|<1$. The first Born approximation corresponds to keeping only the first term in the series, so that (28) becomes

[TABLE]

Likewise, the $M$th Born approximation is given by

[TABLE]

Note that we can also introduce the $N\times N$ discretized $T$-matrix, which is defined by $T_{mn}=\mathcal{T}(\mathbf{r}_{m},\mathbf{r}_{n})$. Eq. (28) can then be written as

[TABLE]

In the above discrete setting, the ISP consists of recovering the matrix $V$ from measurements of $Y$. Our goal is to use suitably modified variants of IHT to solve the ISP when we assume that the desired matrix $V$ is sparse. That is, we cast the ISP as a sparse signal recovery problem and use IHT to solve it. Depending on the form of the approximation for the ISP, we use either linear or nonlinear IHT. In the first Born approximation, this consists of solving the linear system (34). We solve this problem by making use of the linear IHT algorithm (2). In the multiple scattering regime, we solve either (35) or (36) using the nonlinear IHT algorithm (4).

3. Inverse Scattering and Iterative Hard Thresholding

3.1. Linear inverse problem

We now apply the linear IHT algorithm to the ISP in the first Born approximation. Hence, we consider using the linearized forward equation (34), which is valid in weakly scattering regimes. This equation, $AVB=Y$ where $V$ is diagonal, can be converted to a standard linear equation in vector form as

[TABLE]

where $\mathbf{y}$ is the $(N_{s}N_{d})$-dimensional vector formed by stacking the columns of $Y$, $\mathbf{v}$ is the $N$-dimensional vector containing the diagonal entries of $V$, and $\Phi$ is the $(N_{s}N_{d})\times N$ matrix that is defined by

[TABLE]

Here $(mn)$ is a composite index where $1\leq m\leq N_{d}$ and $1\leq n\leq N_{s}$.

One only considers using the IHT algorithm if $\mathbf{v}$ is sparse. Note that this corresponds to the diagonal matrix $V=D_{\mathbf{v}}$ being sparse along its diagonal. Eq. (37) can be solved by a straightforward application of the IHT algorithm (2), but for our purposes, we are interested in applying the algorithm directly to the matrix equation (34). While the two linearized forward equations are identical (as are the resulting formulations of IHT), this is not an empty exercise. Using the full matrix form in (34) yields a more tractable form of IHT, especially as we move to the nonlinear formulation of IHT and its analysis in the next section. Additionally, and perhaps more importantly, these insights lead to computational advantages. Fundamentally, these benefits come from treating the sensing matrix $\Phi$ for inverse scattering problems as two smaller separate sensing matrices. Similar observations are made in [24, 25].

To derive the matrix form of IHT for the linear ISP, we first note that

[TABLE]

Thus $\Phi^{*}\Phi=(A^{*}A)\circ(BB^{*})^{T}$, where  $\circ$  denotes the Hadamard (entrywise) product of two matrices.

To formulate IHT in terms of the diagonal matrix $V_{n}$ instead of the vector $\mathbf{v}_{n}$, we take advantage of the natural interplay between Hadamard products and diagonal matrices, namely

[TABLE]

Combining (3.1) and (40), the matrix product $\Phi^{*}\Phi\mathbf{v}$ in (2) can be rewritten as

[TABLE]

Similarly, for the matrix vector product $\Phi^{*}y$, we have

[TABLE]

We introduce the linear “diagonalizing” operator $\mathcal{D}$ that acts on square matrices according to

[TABLE]

and can write the IHT algorithm of matrix equation (34) by

[TABLE]

where the nonlinear thresholding operator $H_{s}^{\mathcal{D}}$ only acts on the diagonal entries of the resulting matrix.

It is instructive to compare the matrix formulation of IHT in (44) with the standard form of IHT in (2). Recall that the dimensions of the matrices $\Phi,A$, and $B$ are $(N_{s}N_{d})\times N,\ N_{d}\times N$, and $N\times N_{s}$ respectively. Thus the computational complexity of IHT in the form (2) is $O(N_{s}N_{d}N)$. This result is obtained by carrying out two sequential matrix-vector products. In contrast, naively forming the matrix $\Phi^{*}\Phi$ requires $O(N_{s}N_{d}N^{2})$ operations. However, in matrix form (44), the matrix multiplications have complexity $O(\max\{N_{d}N^{2},N_{s}N^{2}\})$. Depending on the relative values of $N_{s},N_{d},$ and $N$, this approach can be computationally advantageous.

We note that the general framework of treating the fundamental unknown as a matrix instead of a vector provides flexibility for algorithmic developments. The unknown $V$, while known to be diagonal, has nonzero off-diagonal terms before the diagonalizing operator $\mathcal{D}$ acts during each iteration. [24] explores a similar algorithm with such flexibility. The definition used in (43) is a particular choice, and there may indeed be better choices for defining the diagonalizing operator.

3.2. Nonlinear inverse problem

We now consider nonlinear IHT, which applies the ISP in the multiple-scattering regime. Returning to the forward equation $Y=A(I-V\Gamma)^{-1}VB$, we add additional corrections beyond the first Born approximation via the Born series. Our goal is to develop an IHT algorithm for Eq. (35).

The derivation of nonlinear IHT follows by making the transformation

[TABLE]

in the linear IHT algorithm. Note that we could also rewrite the forward equation as

[TABLE]

and then make the substitution

[TABLE]

This formulation is justified by the reciprocity of sources and detectors, which allows for interchanging $A$ and $B$, while leaving the forward equation unchanged (with the data matrix $Y$ becoming $Y^{T}$). Note that we can still think of having one larger sensing matrix $\Phi_{\mathbf{v}}$ (which analogous to Eq. (4) is the linearization of $\Phi$ at $V=D_{\mathbf{v}}$), which obeys the equation

[TABLE]

Making either of the above substitutions, we obtain the nonlinear IHT algorithm in the form

[TABLE]

or alternatively

[TABLE]

This algorithm is summarized below as Algorithm 1. As $M\to\infty$, we can replace $\tilde{A}$ and $\tilde{B}$ in Eqs. (49) and (49) with

[TABLE]

or this fully nonlinear IHT algorithm can be written in T-matrix form

[TABLE]

Note that $V_{n}^{-1}$ is not well-defined since it is sparse along its diagonal. We thus compute $V_{n}^{-1}$ only for the nonzero entries of $V_{n}$, setting all other matrix elements of $V_{n}^{-1}$ to zero.

We now remark on the computational complexity of the nonlinear IHT algorithm. Since $V$ is $s$-sparse along its diagonal, the product $\Gamma V$ will be block sparse (up to permutation of rows), with only $s$ nonzero rows. Thus, the multiplication $(\Gamma V)^{m}$ requires $O(Nms^{2})$ operations, and the subsequent formation of $\tilde{A}$ requires $O(sN_{d}N)$. Done in correct order, the remaining matrix products again require $O(\max\{N_{d}N^{2},N_{s}N^{2}\})$ operations. Note that in the case when $M=\infty$, due to the sparseness of $V$, computation of the $T$-matrix has only complexity $O(s^{3})$.

4. Coherence Estimates for Inverse Scattering

In this section we apply Theorem 2.1 on convergence of the nonlinear IHT algorithm to the ISP. We begin by restating Theorem 2.1 using the notation from Section 3.

Theorem 4.1**.**

Let $Y=A(I-V\Gamma)^{-1}VB+E$ and $\{V_{n}\}$ be the sequence of diagonal matrices generated by the nonlinear IHT iteration

[TABLE]

where $\tilde{A}_{n}=A\left(\sum_{m=0}^{M-1}(V_{n}\Gamma)^{m}\right)$. Let $\Phi_{\mathbf{v}_{n}}$ be the linearization about $V_{n}$ corresponding to $\tilde{A}_{n}$. Assume that for all $n\geq 1$, the coherence of $\Phi_{\mathbf{v}_{n}}$ is bounded by $\mu_{0}$. If $V$ is $s$-sparse along its diagonal, then

[TABLE]

where $E_{n}=A(I-V\Gamma)^{-1}VB-\tilde{A}_{n}VB+E$ and $\rho=\mu_{0}(3s+1)$.

To proceed, we must obtain bounds on the coherence suitable for use in the nonlinear IHT algorithm. We investigate two specific cases: double scattering ($M=2$) and multiple scattering ($M=\infty$).

4.1. Mutual coherence of sensing matrices

For the linearized problem, it follows from Eq. (3.1) that normalizing the columns of $\Phi$ is equivalent to normalizing the columns of both $A$ and $B^{*}$. We note that normalizing the column $\Phi_{j}$ is performed by dividing all of its entries by the square root of

[TABLE]

Since we have normalized the columns of $A$ and $B^{*}$ (as well as $\Phi$) and because the coherence of a column-normalized matrix is the largest off-diagonal entry in its Gram matrix, we conclude that

[TABLE]

where we have used the fact that $\Phi^{*}\Phi=(A^{*}A)\circ(BB^{*})^{T}$. To extend this result to the nonlinear case, we use that fact that Eq. (48) implies

[TABLE]

where $\tilde{A}$ and $\tilde{B}$ are defined as in Eqs. (49) and (50). Since the mutual coherence is bounded above by 1, we have the immediate (pessimistic) bound that for any number of nonlinear terms $M$,

[TABLE]

The above bound is most likely not optimal. However, it can be used to easily compute (or accurately estimate) the coherence of higher order linearizations based solely on the geometry of the experiment. Moreover, Eq. (58) is applicable to the nonlinear problem.

For the remainder of the paper, we assume that $A$ and $B^{*}$ have been column-normalized. We require the following result on the coherence of products of matrices.

Proposition 4.2**.**

Let $A$ be an $(L\times N)$ matrix with normalized columns and $H$ be an $(N\times P)$ matrix. Then the mutual coherence of the product $AH$ obeys the following estimate:

[TABLE]

Proof.

We first compute an upper bound on $|(H^{*}A^{*}AH)_{ij}|$.

[TABLE]

where (62) is obtained from (61) by replacing the inner products between columns of $H$ and $A$ respectively with the upper bound on the coherence. We now compute a lower bound on $|(H^{*}A^{*}AH)_{ii}|$.

[TABLE]

Now, by definition of coherence,

[TABLE]

where to obtain the final result we have used inequality

[TABLE]

∎

The following corollary will prove to be of use.

Corollary 4.3**.**

Let $A$ be an $(L\times N)$ matrix with normalized columns, and $H$ be an $(N\times P)$ matrix. Furthermore, assume that $H$ has only $s$ nonzero entries in any column. Then the mutual coherence of the product $AH$ obeys the following bound:

[TABLE]

The simplest case of interest arises when $H=I+VG$ and $V$ is $s$-sparse along its diagonal. Suppose that the first $s$ entries of $V$ are nonzero. Then

[TABLE]

It is clear that the maximum value required for the coherence calculation is given by choosing two columns $H_{i},H_{j}$ where $1\leq i,j\leq s$. If we let $\delta=\|VG\|_{\max}$, an upper bound on the coherence can be explicitly computed to be

[TABLE]

Combining this bound with Corollary 4.3 yields the following result.

Corollary 4.4**.**

Let $A$ be an $(L\times N)$ matrix with normalized columns. Let $V$ be a diagonal $(N\times N)$ matrix, with only $s$ nonzero entries along its diagonal. Let $G$ be an $(N\times N)$ matrix with zeros along its diagonal, and let $\delta=\|VG\|_{\max}<1$. Then the mutual coherence of the product $A(I+VG)$ obeys the following bound:

[TABLE]

We note that analogous results can be derived when $M>1$, but will not be used in this paper.

4.2. Application to second Born approximation ($M=2$)

We now apply Theorem 2.1 to the case $M=2$, which corresponds to double scattering, also known as the second Born approximation. Here we directly apply Corollary 4.4. We then obtain the following estimate on the coherence:

[TABLE]

where we require the third term to be less than $1/(3s+1)$ for all iterations $n$ to guarantee convergence.

The above analysis ignores the behavior of the error term $\|\Phi^{*}_{\mathbf{v}_{n}}\mathbf{e}_{n}\|_{\infty}$ in Theorem 2.1, where $\mathbf{e}_{n}=\Phi(\mathbf{v})-\Phi_{\mathbf{v}_{n}}\mathbf{v}+\bm{\epsilon}$. Next, we investigate the behavior of error for $\Phi^{*}_{\mathbf{v}_{n}}$ as defined by section 2.2. To proceed, we note that the error term can be split into three parts as

[TABLE]

The first term measures the error of the linearization of the forward problem, the second term measures the difference between the linearizations at the true solution and the current iteration, and the last term is the error in the measurement system. Note that this last term can account for the error if the underlying model is not $s$-sparse. Let us examine the first term and make use of the definitions of $\Phi$ and $\Phi_{\mathbf{v}}$. We thus obtain

[TABLE]

Each of the terms in the series is of the form $\|(Y^{(j)}\circ Z^{T})\mathbf{v}\|_{\infty}$ for $j\geq 2$, where

[TABLE]

Letting $y^{(j)}=(Y^{(j)}\circ Z^{T})\mathbf{v}$, we obtain the bound

[TABLE]

The largest entry of $Y^{(j)}$ can be further bounded by

[TABLE]

where $\delta=\|V\Gamma\|_{\rm max}$ and $\gamma=\|V\Gamma\|_{1}$, with analogous definitions for $\delta_{n}$ and $\gamma_{n}$. Because $B$ has ones along the diagonal and all off diagonal terms smaller than $\mu(B)$, we can bound the error term by

[TABLE]

To bound the total error, we sum the above geometric series, assuming $s\delta<1$, with the result

[TABLE]

Now turning to the second term in the error and using (40), we find that

[TABLE]

Once again we use $Y=A^{*}AV\Gamma$ and $Z=BB^{*}(I+\Gamma V_{n})^{*}$, and thus obtain

[TABLE]

We compute the max norm of $Y\circ Z^{T}$ by comparing the maximum diagonal and maximum off-diagonal entries. On the diagonal, $Y_{jj}\leq s\delta\mu(A)$, since the diagonal of $V\Gamma$ vanishes. Likewise, $Z_{jj}\leq 1+s\delta_{n}\mu(B)$. The off-diagonal calculations are identical to those needed for bounding the first term of the error. Therefore, for $j\neq k$, we have

[TABLE]

Since $\mu(A),\mu(B),\delta,\delta_{n}\leq 1$, we have

[TABLE]

We can now restate Theorem 2.1 for the case of the second Born approximation.

Theorem 4.5**.**

Let $Y=A(I-V\Gamma)^{-1}VB+E$ and let $\{V_{n}\}$ be the sequence of diagonal matrices generated by the nonlinear IHT iteration

[TABLE]

where $\tilde{A}_{n}=A\left(I+V_{n}\Gamma\right)$. Let $\gamma=\|\Gamma V\|_{1}<1$, $\delta=\|\Gamma V\|_{\max}$, and $\delta_{n}=\|\Gamma V_{n}\|_{\max}$. If $V$ is $s$-sparse along its diagonal, then for all $n\geq 1$,

[TABLE]

where

[TABLE]

Note that when the coefficient $\rho^{(1)}_{n}$ is less than one for all $n$, the above algorithm exhibits linear convergence. We also note that similar analysis can be performed for linear IHT under the first Born approximation. Here, the second error term in Eq. (73) disappears, and the first error term acquires an additional term in the resulting sum.

Corollary 4.6**.**

Thus, we obtain the error estimate for linear IHT

[TABLE]

4.3. Application to multiple scattering ($M=\infty$)

We now consider the nonlinear version of the IHT algorithm. As before, we assume without loss of generality that $\mu(B^{*})\leq\mu(A)$. We note that the linear error term in Eq. (73) is not present. Therefore, we only need to estimate the remaining terms, which are of the form

[TABLE]

Next, we note that

[TABLE]

Then, we use the Hadamard property (40) of the diagonal matrix $(V_{n}-V)$ in (94) to obtain

[TABLE]

We use the same analysis to bound (4.3) by $\|Y\|_{\max}\|Z\|_{\max}\|\mathbf{v}_{n}-\mathbf{v}\|_{1}$, where $Y=(I-\Gamma^{*}V_{n}^{*})^{-1}A^{*}A(I-V_{n}\Gamma)^{-1}$ and $Z=\Gamma V(I-\Gamma V)^{-1}BB^{*}$. Since $A^{*}A$ and $BB^{*}$ are normalized, we have the bounds

[TABLE]

which results in the final error estimate

[TABLE]

Putting this result into Theorem 2.1, we obtain the following result for the full $T$-matrix IHT algorithm.

Theorem 4.7**.**

Let $Y=A(I-V\Gamma)^{-1}VB+E$ and let $\{V_{n}\}$ be the sequence of diagonal matrices generated by the nonlinear IHT iteration

[TABLE]

where $\tilde{A}_{n}=A\left(I-V_{n}\Gamma\right)^{-1}$. Let $\gamma=\|\Gamma V\|_{1}<1$, $\gamma_{n}=\|\Gamma V_{n}\|_{1}$, $\delta=\|\Gamma V\|_{\max}$, and $\delta_{n}=\|\Gamma V_{n}\|_{\max}$. Assuming that $\mu(B^{*})\leq\mu(A)$, for all $n\geq 1$, then

[TABLE]

where

[TABLE]

5. Numerical simulations

In this section we investigate the linear and nonlinear IHT algorithms in numerical simulations. We also analyze the convergence of the algorithms as characterized by Theorems 4.5 and 4.7. The analysis is carried out in three parts. First, we compare the coherence estimates from section 4.1 with the true coherence for nonlinear IHT. Next, we compare the convergence guarantees in Theorems 4.5 and 4.7 with numerical simulations. Finally, we report simulations that explore the use of the nonlinear IHT algorithm in settings beyond which the convergence estimates in Theorems 4.5 and 4.7 hold.

5.1. Far Field Coherence Estimates

The coherence estimates derived in Section 4.1 are quite general and can be specialized for specific experimental geometries. We recall the following coherence in bounds under the assumption that $\mu(B^{*})\leq\mu(A)$:

• •

Linear: $\mu(\Phi)\leq\mu(A)\mu(B^{*})$

• •

Second Born: $\mu(\Phi_{\mathbf{v}})\leq\mu(A(I+V\Gamma))\mu(B)^{*}\leq\mu^{(1)}(A)\mu(B^{*})$

• •

Fully nonlinear: $\mu(\Phi_{\mathbf{v}})\leq\mu(A(I-V\Gamma)^{-1})\mu(B)^{*}\leq\mu(B^{*})$

where $\mu^{(1)}(A)$ is given by (4.5). In this section, we consider the case of far-field scattering as described in Section 2.2. We begin by computing the coherence of the matrix $A$, as defined by (29). Using (6), we find that

[TABLE]

We assume that the measurement directions $\{\mathbf{\hat{x}}_{m}\}$ are uniformly distributed on the unit sphere. For $N_{s}$ sufficiently large, we can approximate the above sum by an integral:

[TABLE]

where $h$ is the spacing between the voxels. This result is compared to the coherence computed for a finite number of directions in Fig. 2. As seen from Eq. (101), the coherence is a decreasing function in the number of directions of illumination. If the incident directions $\mathbf{\hat{x}}$ are chosen in the same manner as the measurement directions $\mathbf{\hat{x}}$, then $\mu(B^{*})=\mu(A)$.

We now turn to the effects of scattering beyond the Born approximation on the coherence. We proceed by studying the experimental scenario of Fig. 2. In this setting, within the second Born approximation, (4.5) yields the bound $\mu(A(I+V\Gamma))\leq 1$ for all $s$ and $\delta_{n}$. This bound is not useful. Likewise, in the previous analysis of the fully nonlinear algorithm, we use the same trivial bound $\mu(A(I-V\Gamma)^{-1})\leq 1$. We look to improve these upper bounds for simple cases.

5.1.1. The case $s=1$

We consider the simplest case of a single scatterer ($s=1$) at the position $\mathbf{r}_{*}$ with potential $\eta_{0}$. Here $V\Gamma$ has only one nonzero row and thus $(I-V\Gamma)^{-1}=(I+V\Gamma)$. By (29) the entries of this matrix are given by

[TABLE]

The coherence is now approximated for large values of $N_{s}$ by

[TABLE]

where, in accordance with Eq. (103), the function $f$ is given by

[TABLE]

Computing the above integrals, we obtain

[TABLE]

where the functions $D_{1}$ and $D_{2}$ are defined as

[TABLE]

This formula makes evident that the coherence of $A(I+V\Gamma)$ is a scaled version of $\mu(A)$ with an additive term. The sinc function (which gives $\mu(A)$) is the only term that depends on the distance between $\mathbf{r}_{j}$ and $\mathbf{r}_{\ell}$. All other terms depend on the details of the scatterer.

The dependence of the coherence on $|\mathbf{r}_{j}-\mathbf{r}_{*}|$ is shown in Fig. 3 for several values of $\eta_{0}$. We have verified numerically that the above maximum is obtained by setting $|\mathbf{r}_{j}-\mathbf{r}_{\ell}|=h$ and $|\mathbf{r}_{j}-\mathbf{r}_{*}|=|\mathbf{r}_{\ell}-\mathbf{r}_{*}|$. We see that the coherence is given by the maximum value of the curve, as indicated by the dots. Interestingly, for stronger scattering ($\eta_{0}=0.2$) the maximum is achieved by taking points close to the scatterer (with the closest distance being $h$), while for weaker scattering, the maximum is obtained at approximately $2h$. Thus stronger scattering increases the coherence and there is a trade off between coherence and nonlinearity, as further illustrated below.

5.1.2. The case $s>$1

For $s>1$, let $\mathcal{I}$ index the locations of the scatterers, with $|\mathcal{I}|=s$. If each scatterer has the same potential $\eta_{0}$, the coherence is given by Eq. (104) with

[TABLE]

which is derived by computing the entries of the matrix $A(I+V\Gamma)$. The required integrals can still be computed analytically, but become quite cumbersome. Thus, we will resort to numerical calculations. Fig. 4 displays the numerically computed coherences for varying numbers of scatterers. In each experiment, $s$ scatterers are uniformly randomly placed inside the unit cube, which has been discretized into a $10\times 10\times 10$ voxel grid, with $N_{s}=500$ and $kh=1.885$. We compare the case when $\eta_{0}$ is fixed for all $s$, and the case when $\eta_{0}$ changes but $\|V\Gamma\|_{1}$ is fixed. In the latter case, the coherence decreases as $s$ increases. Moreover, there is a minimal difference between the coherence for the first and second Born approximations. We also compare the results of Figs. 3 and  4, where for $s=1$ and $\eta_{0}=0.1$, we obtain the analytically computed nonlinear coherence value of 0.540 and the numerically computed value of 0.539. Overall this observation confirms the conservative nature of the mutual coherence estimates.

5.2. Comparison with convergence theorems

We consider the same experimental setup as above, but at higher frequencies with $kh=1.885$. Evidently, convergence of the nonlinear IHT algorithm is not guaranteed. We also note that in the linear case with $\mu(A)=\mu(B^{*})=0.505$, convergence is not guaranteed, since no choice of $s$ satisfies the condition $\rho=\mu(A)\mu(B^{*})(3s+1)<1$. In this section, we investigate some examples that shed some light on the theoretical convergence results. We will see that the theory for the nonlinear IHT algorithm is quite conservative, and that even when the theory does apply, faster convergence than predicted is obtained in practice.

Since the previous choice of parameters fell outside of the theoretical convergence regime, we consider two experiments with extreme choices of parameters for which the theoretical convergence guarantees hold. Here the domain $\Omega$ is discretized into a 10$\times$10$\times$10 grid. For each experiment, three voxels are chosen uniformly at random to have nonzero potential $\eta_{0}$. The sources consist of 400 incident plane waves at uniformly spaced angles over the entire sphere, with the scattering amplitude being measured in the far field for the same angles. Forward data is generated using the coupled-dipole method [27].

For the first experiment, $\Omega$ is a cube with sidelengths $4.75\lambda$, which makes the dimensionless parameter $kh=2.98$. This high frequency setup results in coherence values of $\mu(A)=\mu(B)=0.0525$. The potential $\eta_{0}$ is set to $10^{-5}$. After choosing the three scatterer locations at random (uniformly), we find that $\delta=\|\Gamma V\|_{\max}=0.0395$ and $\gamma=\|\Gamma V\|_{1}=0.0460$. For the second experiment, the cubic region is taken to be smaller with sides of length 4.55$\lambda$ with $kh=2.858$. We then obtain $\mu(A)=\mu(B)=0.098$. The potential has also been increased to $\eta_{0}=10^{-4}$, which gives $\delta=\|\Gamma V\|_{\max}=0.987$ and $\gamma=\|\Gamma V\|_{1}=0.1299$ after choosing the three scatterer locations uniformly at random. For both experiments, the theoretical error bounds $\|\mathbf{v}-\mathbf{v}_{n}\|_{1}$ for linear, second Born approximation and nonlinear IHT are shown in Fig. 5. Also shown, for comparison, are numerical results obtained from the reconstruction algorithms under the same conditions with the thresholding parameter set to $s=3$. In Fig. 5 we plot the dependence of the $L^{1}$ error on the number of iterations of the reconstruction algorithms for both experiments. The results from numerical reconstructions are compared to Theorems 4.5 and 4.6. As may be expected, we find that the theoretical error bounds are conservative.

5.3. Numerical Reconstructions

In this section we illustrate the use of the nonlinear IHT algorithm in settings where convergence cannot be expected. Two models are used to test the algorithm: two spherical scatterers of constant potential (model 1) and one larger spherical scatterer with radially varying potential (model 2). The scatterers are contained in the cubic domain $\Omega$ with sides of length $5\lambda$. In model 1 the scatterers have potential $\eta_{0}$, radii $\lambda/2$, and are separated by a distance of $3\lambda/2$ between their centers. In model 2, the scatterer has radius $5\lambda/4$ and potential that decreases linearly from the value $\eta_{0}$ to 0 at its center. Tomographic slices for each model are shown in Fig. 6.

Forward data was generated using the coupled-dipole method [27], using a $25\times 25\times 25$ discretization. The incident and measurement directions were chosen to be spaced uniformly on half of the unit sphere, with $N_{s}=N_{d}=225$ and $kh=1.2566$. Gaussian white noise was added to the data at a level of 1%. Here the coherence is given by$\mu(A)=\mu(B)=0.764$.

Reconstructions were conducted using a discretization of $\Omega$ into a $21\times 21\times 21$ cubic mesh, thereby avoiding so-called inverse crime. Four versions of the IHT algorithm were investigated: linear ($M=1$), second Born ($M=2$), third Born ($M=3$), and fully nonlinear ($M=\infty$). Each algorithm was run for 100 iterations. In all cases, the iterations are initialized by $V_{0}=0$. For model 1, the threshold level $s$ for the thresholding operator $H_{s}$ was set to 230/9,261 = 2.48%. Note that the threshold limit in $H_{s}$ is not equal to the actual sparsity $s$, which is 0.84% sparse. While not explained by the theoretical analysis, this level of sparsity leads to better performance, especially in the presence of noise. Reconstructions of the central slices (corresponding to a depth of $2.5\lambda$) of model 1 using the linear and fully nonlinear algorithms are displayed in Table 1. At the weakest level of scattering, both algorithms perform nearly identically. However, as the strength of scattering is increased, the linear reconstruction degrades, while the nonlinear algorithm reconstructs the scatterers quite well. At the highest level of scattering, both algorithms fail. Convergence plots for all algorithms are also displayed. For each iteration, the error $Y_{err}$ measures the relative error between the true data $Y$ and data that would be generated by the current reconstruction $V^{(\rm rec)}$. Consistent with with Eq. (28), the error is defined as

[TABLE]

Next we consider the reconstruction of model 2 using the same parameters as for model 1 with $\|V\Gamma\|_{2}=0.9813$. This model is substantially less sparse, with 6.5% sparsity. A thresholding limit of 1,800 was chosen for all reconstructions. Reconstructions of the central slices of model 2 using the linear and fully nonlinear algorithms are displayed in Table 2. Tomographic reconstructions are also shown in Fig. 7. Central slices of the reconstructions for varying number of terms in the Born series are shown in Fig. 8. The convergence of the algorithms is illustrated in Fig. 9. We note that as in model 1, the linear reconstruction algorithm finds the support well and higher orders of scattering recover the interior susceptibilities.

Using the same experimental setup ($N_{s}=N_{d}=225$), $\Omega$ was discretized into 1,000 voxels, where $s$ voxels were chosen uniformly at random to have potential $\eta_{0}$. 50 iterations of the IHT algorithms were performed for $M=1,2,3,4,\infty$. Each experiment was performed with 500 realizations with the hard thresholding limit set to the sparsity level $s$. We define a successful reconstruction as exactly recovering the support of all $s$ voxels. Fig. 10 shows the success rate for two choices of $\eta_{0}$ with varying levels of sparsity.

We note that since the targets are quite sparse, the second Born approximation accounts for the majority of the improvement of the nonlinear reconstruction over linear IHT.

6. Discussion

We have considered the application of nonlinear IHT to inverse scattering. The convergence and error of the method was analyzed by means of coherence estimates and compared to numerical simulations. It was found that in the multiple-scattering regime, the performance of nonlinear IHT is superior to linear IHT. Moreover, numerical evidence suggests that the coherence estimates we have obtained are relatively conservative. Future work will be concerned with average case coherence estimates, which may be anticipated to lead to a more accurate characterization of the convergence of the nonlinear IHT algorithm. It would also be of interest to explore the role of near-field measurements in improving the resolution of reconstructed images. Other areas of future interest include the extension of our results to the setting of the Bremmer series, which is well adapted to the study of layered media and geophysical applications.

Appendix A Proof of Theorem 2.1

The proof follows the general ideas in [29]. First, we assume without loss of generality that the columns of the sensing matrix are normalized so that their $\ell_{2}$ norms are unity. Then, we use the following lemmas:

Lemma A.1**.**

Let $A\in\mathbb{C}^{m\times n}$ be normalized such that $\|A_{i}\|_{2}=1$ for all $i=1,\dots,n$. Then,

[TABLE]

Lemma A.2**.**

Let $\mathbf{z}_{n+1}=\mathbf{x}_{n}+\Phi^{*}_{\mathbf{x}_{n}}(\mathbf{y}-\Phi_{\mathbf{x}_{n}}\mathbf{x}_{n})$. Let $I^{*}$ be the index set of nonzero entries of $\mathbf{x}$, $I^{n}$ be the index set of the $s$ largest entries of $\mathbf{z}_{n}$, and $I^{n}_{+}$ be the index set of the $s+1$ largest entries of $\mathbf{z}_{n}$. Let $S_{n}=I^{n}_{+}\bigcup I^{n-1}\bigcup I^{*}$. Let $\mathbf{w}_{n}$ be the subvector of $\mathbf{z}_{n}$ with components restricted to $S_{n}$. Then, the following inequality holds

[TABLE]

The proofs of these lemmas can be found in [29].

Using $\mathbf{z}_{n+1}=\mathbf{x}_{n}+\Phi^{*}_{\mathbf{x}_{n}}(\mathbf{y}-\Phi_{\mathbf{x}_{n}}\mathbf{x}_{n})$, we find that

[TABLE]

Let $A=I-\Phi^{*}_{\mathbf{x}_{n}}\Phi_{\mathbf{x}_{n}}$, where by hypothesis the linearized matrix $\Phi_{\mathbf{x}_{n}}$ has normalized columns. Then,

[TABLE]

where the last inequality follows from Lemma A.1. The above inequality holds for any subvector of $\mathbf{z}_{n+1}-\mathbf{x}$, so we have

[TABLE]

By combining the inequalities in Eqs. (111) and (116), we obtain

[TABLE]

The theorem follows by iterating the above estimate. ∎

Appendix B RIP Analysis of Nonlinear IHT for Scattering

Here we prove a similar result to Theorem 4.7 using RIP rather than coherence estimates. We first state a related result, which has been slightly simplified for our purposes.

Theorem B.1**.**

[2]** Let $\mathbf{y}=\Phi(\mathbf{x}_{*})$ with $\mathbf{x}_{*}$ $s$-sparse and let $\{\mathbf{x}_{n}\}$ be the sequence generated by the iteration

[TABLE]

where $\Phi_{\mathbf{x}_{n}}$ is the linearization of $\Phi$ at $\mathbf{x}_{n}$. Assume that for all $n\geq 1$, the linearizations $\Phi_{\mathbf{x}_{n}}$ satisfy the RIP for all $s$-sparse $\mathbf{z}_{1}$ and $\mathbf{z}_{2}$

[TABLE]

Further assume that $\Phi$ and $\Phi_{\mathbf{x}_{n}}$ satisfy

[TABLE]

Then, the following inequality holds for all $n\geq 1$

[TABLE]

•

We note that this theorem is conceptually similar to Theorem 2.1, since we require the linearizations to obey a common RIP bound, as well as a bound on the error of the linearizations.

We will assume that we are given RIP bounds for the sensing matrices $A$ and $B$. That is, we let $A$ and $B$ satisfy RIP bounds of order $s$ with constants $\delta_{s}^{A}$ and $\delta_{s}^{B}$, respectively. By invoking the results of [20], the RIP constant for the linear matrix $\Phi$ defined by (38) is given by $\delta_{s}^{2}$, where $\delta_{s}=\max\{\delta_{s}^{A},\delta_{s}^{B}\}$. This result holds for linearized matrices of the form (49). Now, given RIP bounds for the linear sensing matrices, we obtain a bound on the RIP constant for the linearizations. Again, without loss of generality, we will use (49) for the nonlinear IHT, and attach the extra factors in the linearization to the matrix $A$. Let us first consider the linearization using the second Born approximation. Results for the higher order terms follow immediately. Next, we derive an RIP type bound for the matrix $A(I+VG)$, where $V$ is a diagonal matrix that is $s$-sparse. Using the triangle inequality, we have

[TABLE]

Making use of this result, we obtain the inequalities

[TABLE]

where we have introduced the RIP constant for $A$ of order $2s$. Since $\mathbf{x}$ is an $s$-sparse vector and $V$ is $s$-sparse along its diagonal, the vector $(I+VG)\mathbf{x}$ is at most $2s$-sparse. For higher order terms, $\|VG\|_{2}$ is replaced by the $\sum_{j=1}^{M}\|VG\|_{2}^{j}$. Thus, for any order $M$, we have

[TABLE]

This result gives us the constants $\alpha$ and $\beta$ in Theorem B.1. In particular, by squaring the RIP constant, we have

[TABLE]

Now, we bound the error term

[TABLE]

Introducing matrix notation, this term becomes

[TABLE]

which can be rewritten as

[TABLE]

Treating the “vector” as the diagonal matrix $V_{n}-V$, to compute an RIP bound, we are interested in the two matrices $A(I-V_{n}G)^{-1}$ and $\left((I-VG)^{-1}VB\right)^{*}$. The RIP bound for the first matrix follows from Eq. (124):

[TABLE]

Similarly, for the second matrix, we have

[TABLE]

where $m$ is the smallest (in absolute value) nonzero entry of $V$. Once again using the results of [20], the RIP constant is given by

[TABLE]

which provides the constant $C$ in the theorem. We can now summarize our results in the following theorem.

Theorem B.2**.**

Let $\mathbf{y}=\Phi(\mathbf{v})$ be the scattering model, where $\Phi(\mathbf{v})={\rm vec}(AV(I-\Gamma V)^{-1}B)$. Let $\gamma=\|\Gamma V\|_{2}<1/2$. Suppose $A$ and $B$ satisfy the RIP on $2s$-sparse vectors with constants $\delta_{2s}^{A}$ and $\delta_{2s}^{B}$, respectively. Define $\delta_{2s}=\max\{\delta_{2s}^{A},\delta_{2s}^{B}\}$ and

[TABLE]

Then the nonlinear IHT algorithm converges if

[TABLE]