An Accelerated Nonlinear Contrast Source Inversion Scheme for Sparse Electromagnetic Imaging

An efficient nonlinear contrast source inversion scheme for electromagnetic imaging of sparse two-dimensional investigation domains is proposed. To avoid generating a sequence of linear sparse optimization problems, the non-linearity is directly tackled using the nonlinear Landweber (NLW) iterations. A self-adaptive projected accelerated steepest descent (A-PASD) algorithm is incorporated to enhance the efficiency of the NLW iterations. The algorithm enforces the sparsity constraint by projecting the result of each steepest descent iteration into the $L_{1}$ -norm ball and selects the largest-possible iteration step without sacrificing from convergence. Numerical results, which demonstrate the proposed scheme’s accuracy, efficiency, and applicability, are presented.


I. INTRODUCTION
Developing numerical schemes for solving the electromagnetic (EM) inverse scattering problem on (spatially) sparse investigation domains has become an active research topic in the last decade.This is simply because such sparse domains naturally occur in many applications such as non-destructive testing, crack detection, and see-through wall imaging [1], [2].Indeed, recently, sparsity-constrained regularization schemes have been successfully developed for nonlinear EM imaging [3]- [8].These algorithms use L 0 /L 1 -norm of the solution as the penalty term in the nonlinear minimization problem required by the inversion process.Examples of these algorithms include distorted Born iterative methods (DBIMs) that include also L 2 -norm of the solution in the penalty term [3], [4], a preconditioned inexact Newton contrast-source (INCS) scheme [6], and a scheme that makes use of nonlinear Landweber (NLW) iterations [7].All these methods incorporate different thresholding functions to enforce the sparsity constraint [9]- [11].Numerical results presented in the these papers demonstrate that the sparsity (L 0 / L 1 -norm) constrained inversion schemes produce sharper and more accurate reconstructions than "traditional" inversion schemes making use of only L 2 -norm-constrained regularization.
Unlike the INCS scheme, the NLW iterations do not call for generation and solution of a sequence of linear sparse optimization problems (i.e., tackle the non-linearity directly) and require selection/tweaking of a smaller set of parameters [12].
On the other hand, the NLW iterations, both the traditional arXiv:1911.06020v1[eess.SP] 14 Nov 2019 version used for L 2 -norm-constrained regularization and the thresholded version, often suffer from slow convergence which prohibits their direct application to electrically large investigation domains.
To overcome this bottleneck of slow convergence, a projected accelerated steepest descent (PASD) algorithm has been used to solve the three-dimensional (3D) EM inverse scattering problem [8], [13], [14].To enforce the sparsity constraint, this algorithm projects the result of each steepest descent iteration into the L 1 -norm ball.The steepest descent iterations are faster than NLW iterations since they allow for selection of a larger iteration step without sacrificing from accuracy and convergence.Consequently, the resulting PASD algorithm achieves significant speed up over the thresholded NLW iterations.However, even with the increased convergence, the computational cost of every iteration is still a limiting factor in applying the PASD scheme to electrically large investigation domains.This high computational cost is a result of the matrix inversions required to compute the forward operator and its (first-order) Frechet derivative [8].
To this end, in this work, the overall computational cost of the nonlinear EM inversion is further reduced by incorporating a self-adaptive version of the PASD algorithm [15] and the contrast source (CS) formulation of EM scattering [16], [17].
The self-adaptive PASD scheme (A-PASD) chooses the largest possible iteration step in a recursive manner (with smaller number of operations) and increases the convergence of the original PASD algorithm.The CS formulation uses both the contrast and equivalent currents as the unknowns to be solved for in the nonlinear optimization.Even though the total number of unknowns is increased, the forward operators and its Frechet derivative do not call for matrix inversions.The resulting A-PASD-CS scheme benefits from the faster convergence rate of the A-PASD algorithm and the lower computational cost (per iterations step) that comes with the CS formulation.

A. Contrast Source Formulation
Let S represent a 2D investigation domain embedded in an unbounded homogeneous medium.The permittivity and permeability in S and the background medium are {ε(r), µ 0 } and {ε 0 , µ 0 }, respectively.S is surrounded by N T number of line sources (Fig. 1).It is assumed that S is excited by these sources individually, i.e., by one of them at a time.Let E inc i (r), i = 1, .., N T , represent the transverse-magnetic (to z), TM z , polarized electric field, generated by the i th source.Upon excitation, equivalent volume electric current J i (r) is induced in S; J i (r) generates the scattered electric field E sca i (r), which is expressed as Here, G(r, r ) = H 2 0 (k 0 |r − r |)/(4j) is the 2D scalar Green function, k 0 = ω √ ε 0 µ 0 is the wavenumber in the background medium, and ω is the frequency.Additionally, J i (r) is related to the total electric field E i (r) where τ (r) = ε(r)/ε 0 − 1 is the contrast.Substituting (1) for E sca i (r), r ∈ S, and multiplying the resulting equation with τ (r) yield: The coupled set of equations ( 1) and ( 2) are known as the CS formulation [16], [17].Let are the measurement locations.Then, the CS equations ( 1) and (2) can be numerically solved for unknowns J i (r) and τ (r) given the next two sections.

B. Discretization
To facilitate the numerical solution, first, S is discretized into N number of square elements.Let r n , n = 1, . . ., N represent the centers of these elements.Then, the unknowns J i (r) and τ (r) are approximated as where where 2) and evaluating the resulting equations at r m , m = 1, . . ., N yield: where and operator D[.] generates a diagonal matrix with entries equal to the entries of the vector at its argument.Equations ( 4) and ( 5) are cascaded together for all excitations i = 1, .., N T in a more compact form as: where Here, subscript 't' represents the transpose operation.

C. Sparsity Regularized Nonlinear Optimization Problem
Equation ( 6) is a nonlinear equation in z; but also it is illposed and cannot be solved accurately/efficiently without using a regularization method [1], [2].In this work, it is assumed that the investigation domain is sparse, i.e., many entries of τ , and consequently the same entries of Ji , i = 1, .., N T , are zero.Therefore, to alleviate the ill-posedness, ( 6) is cast in the form of a sparsity-constrained nonlinear optimization problem: In (7), the nonlinear least squares minimization, i.e., the data 2 is constrained by the condition z 0 ≤ l 0 , which ensures that the total number of non-zero entries in z (given by z 0 ), is less than the positive integer l 0 .Solution of (7) provides the sparsest result possible, however, since it is a non-convex optimization problem, this solution is NPhard and computationally not feasible [10], [18].Best convex approximation to (7) is obtained by constraining the data misfit with the L 1 -norm of z: In ( 8), l 1 is a positive real number that should be estimated based on the prior knowledge of the sparsity level of z and the values of its entries.A plethora of algorithms has been developed to solve the L 1 -norm constrained nonlinear optimization problem [9]- [11].In this work, the scheme of nonlinear Landweber iterations (NLW) is selected from this group of algorithms to solve (8); at every iteration of this scheme, a thresholding function is applied to enforce the sparsity constraint.The NLW scheme is preferred here because it requires a smaller number of parameters to be "tuned" heuristically in comparison to the inexact Newton and Born iterative methods that have been previously developed for solving the inverse scattering problem on sparse domains [5], [6].However, it has also been shown that the classical NLW scheme converges fast at the beginning of the iterations, but then it overshoots the L 1 -norm penalty and takes a large number of iterations to correct back [13].

Scheme
In this work, the slow convergence of the NLW scheme is alleviated using an self-adaptive version of the PASD algorithm [14], [15].The original PASD scheme achieves acceleration by confining the solution search within a particular L 1norm ball, while maintaining a large iteration step size without sacrificing from accuracy and convergence [14].Its selfadaptive version, which is proposed in [15] and abbreviated as A-PASD in this paper, further increases the convergence by controlling the step size in a recursive/adaptive manner.This approach starts with a larger step size and decreases it only when necessary.Application of the A-PASD to the solution of the nonlinear problem (8) yields in the following algorithm: Step 1 : Initialize l 1 , α, ψ, µ ∈ (0, 1), δ ∈ (0, 1) ρ ∈ (0, 1), λ (k≥1) > 0, γ (0) = β (0) > 0, and z(0) Step 2 : for k = 0, 1, 2, ....
They are estimated heuristically by running some numerical tests on a given problem set-up before the A-PASD iterations start [14], [15].This is not a costly operation since these parameters are estimated only once and are not updated during the iterations.
Step 2 and its sub-steps describe the iterations."do loop" that starts at Step 2.2 adjusts the iteration step size in an adaptive/recursive manner until the condition [15], [19], [20] is satisfied.
Step 2.2.2 is the fixed point iteration of the solution, where It is important to note here that z(k) contain samples of the contrast τ (r) and equivalent currents J i (r), which have values that are orders of magnitude different from each other.The effect of this scaling mismatch is observed in the singular values of ∂ z T † (•), and consequently, the iterations converge very slowly.In this work, a left-right diagonal preconditioning scheme (that is similar to the one used in [6]) is applied at Step 2.2.2 to alleviate the effect of this scaling mismatch and increase the convergence rate.
The purpose of Step 2.3 is to increase the iteration step size if it gets too small.This is checked by the condition Note that there could be several positive sequences satisfying λ (k≥1) > 0; in this work, λ (k) = λ 0 , where λ 0 is a positive constant.
Note that the mathematical derivation and justification of the conditions, steps, and projection operation used by the A-PASD algorithm are not provided in here but can be found in [13]- [15].In this work, this algorithm is applied to accelerate the solution of the EM inverse problem constructed using the CS formulation.The numerical results presented in the next section demonstrate that the resulting inversion scheme is indeed faster and more accurate than existing CSbased EM inversion schemes.

III. NUMERICAL RESULTS
This section demonstrates the accuracy and efficiency of the proposed scheme via numerical experiments.In all examples considered here, Four different EM inversion schemes are compared: (i) A CS inversion scheme that uses a multiplicative total variation term together with the cost function to enhance the reconstruction quality of piece-wise constant profiles [16], [17].
This scheme also uses a preconditioned conjugate gradient method to accelerate the convergence.In the rest of this section, this scheme is referred to as ECSI-TV.(ii) The ECSI-TV scheme that enforces an additional positivity constraint on the solution.This scheme is referred to as ECSI-TV-P.(iii) A preconditioned inexact Newton scheme that incorporates the CS formulation.This scheme is referred to as SP-INCS [6], and (iv) A-PASD-CS algorithm proposed in this work.For all simulations, the quality of reconstruction is measured using where τt stores the samples of the contrast reconstructed at the iteration corresponding to the execution time t.

A. Coaxial
In this example, the investigation domain of dimension   converges faster than the other three schemes and produces a more accurate reconstruction at the point of convergence.

B. Austria
In this example, the investigation domain of dimension

IV. CONCLUSION
The A-PASD algorithm and the CS formulation are incorporated and the resulting scheme is used for solving the nonlinear EM inverse scattering problem on sparse 2D investigation domains.This scheme benefits from the faster convergence rate of the A-PASD and lower computational cost of the CS formulation.Numerical results demonstrate that the proposed algorithm is faster than existing CS-based inversion schemes.
Integration of the TV regularization with the proposed scheme and its extension to 3D investigation domains is currently underway.
and p n (r) is the pulse basis function on element n with support S n .p n (r) is nonzero only for r ∈ S n with unit amplitude.Substituting the first expansion in (3) into (1) and evaluating the resulting equations at r R m , m = 1, . . ., N R yield:

At Step 1 , 2 2 2
several parameters are initialized: α and ψ should be selected to satisfy α ≥ sup z∈B ∂ z T (z) 2 and ψ ≥ sup {z,ū}∈B 2 ∂ 2 z T (z, ū) 2 ū , where ∂ z T (z) and ∂ 2 z T (z, ū) are the first-and second-order Frechet derivatives [6] and B = { z 1 ≤ l 1 } is a ball of radius l 1 [14].Parameters µ, γ (k) , and λ (k) help to reduce the number of "trials" in finding the iteration step size β (k) /r (see Steps 2.2.1, 2.2.2, and 2.3) [15].Parameter r = max 2α, 2ψ Γ{z (0) } , where Γ(z) = 0.5 ȳ − T (z) 2 is the data misfit, and helps in adjusting the iteration step size [see Step 2.2.2, (9), and (10)].Parameters δ and ρ, which are on the right hand sides of ( Step 1, vector x stores the absolute of value of the entries of the input vector z sorted in descending order.Steps 2 and 3 are needed to compute the threshold level used at Step 4. Here, Thr χ (•) is the complex soft-thresholding function defined as[8] generated synthetically: First, (5) with τ ref is solved for Ji , then Ji is inserted into (4), and finally 25dB Gaussian noise is added to the result to yield E mea i (r R m ).Here, {τ ref } n = τ ref (r n ), n = 1, . . ., N , are the samples of the actual contrast τ ref (r) being reconstructed.

7. 5
m×7.5 m includes a "coaxial"-shaped scatter.The relative permittivity of the inner cylinder and the outer shell is 1.8 and 2.5 respectively.The investigation domain is discretized using N = 2500 cells with dimension ∆d = 0.15 m.The parameters of the transmitter-receiver configuration are N T = 8 , N R = 16, and f = 125 MHz.The investigation domain (as represented by τ ref ) and transmitter-receiver configuration are shown in Fig. 2(a).The sparseness level (the ratio of number of non-zero entries in τ ref to N ) is 4%.The parameters of the A-PASD-CS scheme are provided in TableI.

Fig. 2 (Fig. 2 .
Fig.2(e) plots err t versus execution time t for all the schemes (t-axis is in log-scale).These results show that A-PASD-CS
1.25 x 10 −2 0.2 0.5 0.8 0.25 Lossy Austria 0.39 1.0 x 10 −2 0.2 0.5 0.8 0.25 the point of convergence.But ECSI-TV-P can not be used for investigation domains involving dispersive/lossy scatterers (i.e., when τ (r) ref is complex) while A-PASD-CS maintains its convergence rate and accuracy for these cases.This is demonstrated next.The contrast of the Austria scatterer is made complex by introducing a conductivity of 5m/S.All other configuration parameters are kept the same.Fig. 4 (a) and (b) show Re{τ ref } and Im{τ ref }, respectively.The updated A-PASD-CS parameters are provided in Table I .

Fig. 4 (
Fig.4 (c)-(e) and (d)-(f) show Re{τ t=25 s }-Im{τ t=25 s } obtained using ECSI-TV and A-PASD-CS, respectively.As before, the images are obtained at the time of convergence and the corresponding error levels are err t=25 s = 60% and err t=25 s = 39%.Fig. 4(f) plots err t versus execution time t for both the schemes.These results show that A-PASD-CS indeed maintains its convergence rate and accuracy for investigation domains with complex contrast.

Fig. 4 .
Fig. 4. Investigation domain with the lossy Austria scatterer as represented by (a) Re{τ ref } and (b) Im{τ ref }.(c) Real part and (d) imaginary part of the solution obtained by ECSI-TV.(e) Real part and (f) imaginary part of the solution obtained by A-PASD-CS.(g) Reconstruction error errt versus execution time t for the two schemes.

TABLE I PARAMETERS
OF THE A-PASD-CS SCHEME APPLIED TO THE TWO EXAMPLES.