# Efficient combination of a 3D Quasi-Newton inversion algorithm and a vector dual-primal finite element tearing and interconnecting method

Ivan Voznyuk, Amelie Litman, Hervé Tortel

► To cite this version:

Ivan Voznyuk, Amelie Litman, Hervé Tortel. Efficient combination of a 3D Quasi-Newton inversion algorithm and a vector dual-primal finite element tearing and interconnecting method. Inverse Problems, 2015, 31 (8), pp.085005. 10.1088/0266-5611/31/8/085005. hal-01279541

HAL Id: hal-01279541

https://amu.hal.science/hal-01279541

Submitted on 5 Nov 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Efficient combination of a 3D Quasi-Newton
inversion algorithm and a vector Dual-Primal Finite
Element Tearing and Interconnecting method

I. Voznyuk¹, A. Litman¹, H. Tortel¹

1 Aix Marseille Université, CNRS, Centrale Marseille, Institut Fresnel, UMR 7249, 13397 Marseille, France

E-mail: amelie.litman@fresnel.fr

**Abstract.** A Quasi-Newton method for reconstructing the constitutive parameters of three-dimensional penetrable scatterers from scattered field measurements is presented. This method is adapted for handling large-scale electromagnetic problems while keeping the memory requirement and the time flexibility as low as possible. The forward scattering problem is solved by applying the Finite-Element Tearing and Interconnecting Full-Dual-Primal (FETI-FDP2) method which shares the same spirit as the Domain Decomposition methods for Finite Element methods (FEM). The idea is to split the computational domain into smaller non-overlapping sub-domains in order to simultaneously solve local sub-problems. Various strategies are proposed in order to efficiently couple the inversion algorithm with the FETI-FDP2 method: a separation into *permanent* and *non-permanent* subdomains is performed, iterative solvers are favorized for resolving the interface problem and a *marching-on-in-anything* initial guess selection further accelerates the process. The computational burden is also reduced by applying the adjoint state vector methodology. Finally, the inversion algorithm is confronted to measurements extracted from the 3D Fresnel database.

Submitted to: *Inverse Problems*
# 1. Introduction

Quantitative inverse scattering algorithms attempt to estimate from scattering experiments the physical parameters and features (position, form, size and complex permittivity) of an unknown target. Quantitative microwave imaging is faced with two main challenges: its ill-posedness and its non-linearity. The former jeopardizes the robustness of the reconstruction algorithms and the quality of the results and the latter results in a high computational cost. Indeed, unlike qualitative methods, quantitative ones solve the exact non-linear electromagnetic inverse problem, which requires the solution of a system of coupled equations. For this aim, either a global optimization procedure is applied [1, 2] or a local one [3, 4, 5], into which the forward solver plays a key role.

Various works take profit of the Finite Element Method (FEM) in order to solve inverse problems in different scientific domains, such as in optical imaging [6], electroencephalography imaging [7], electrocardiographic imaging [8], elasticity imaging [9], electrical impedance tomography [10], eddy-current imaging [11] and, of course, in microwave imaging [12, 13, 14, 15, 16, 17, 18, 19]. Thanks to its flexibility and its capabilities of managing complex geometrical configurations, this type of modelling scheme is more than appropriate to tackle biomedical and non-destructive testing applications. Nevertheless, the numerical resolution of the Helmholtz equation in heterogeneous media at high wave number is a challenging problem for the classical finite element method. Indeed, the discretization mesh in FEM is directly proportional to the wavelength $λ_0$. When the frequency increases, the number of field unknowns increases very rapidly, especially in three-dimensional configurations. For example, a $30λ_0^3$ computational domain corresponds to approximately $0.510^6$ field unknowns and $0.510^6$ permittivity unknowns, while a domain of $90λ_0^3$ meshed this time for a frequency which is simply doubled leads to almost $1.510^6$ field unknowns and $1.410^6$ permittivity unknowns. Even if the associated matrix system is sparse, its inversion as well as its storage rapidly overwhelms even the largest resources that are currently available.

A way to overcome these difficulties is to apply a domain decomposition (DDM) technique [20]. The principal idea of DDM is to split the entire computational domain into smaller non-overlapping sub-domains and to solve a sequence of similar sub-problems on these sub-domains. Among a variety of DDMs, the Finite-Element Tearing and Interconnecting (FETI) method [21, 22] and its electromagnetic counter part, the Finite-Element Tearing and Interconnecting-Dual Primal Electromagnetic method (FETI-DPEM) [23, 24, 25, 26, 27, 28, 29, 30], are shown to be powerful techniques with an excellent scalability.

Very few microwave inversion schemes take profit from the FETI-DPEM method to handle large-scale problems, apart eventually [31] where a two-dimensional (2D) configuration is investigated. In this work, we focus on the three-dimensional (3D) case, because the computational burden here is more important than in 2D and leads us to the necessity of well-defined numerical strategies. To that end, an updated version of
the FETI-DPEM method, named FETI-FDP2, has been developed and thoroughly
tested [32, 33]. Our aim here is to efficiently coupled it with an inversion scheme.

We have set our sights on quantitative imaging algorithms, without any use of *a-priori* information, to recover the value of the complex permittivity in every point of a given investigation domain. As the number of the permittivity unknowns is increasing drastically at high wave number, we must tune the FETI-FDP2 method in order to make this process more flexible in terms of time and memory consumption. At first, we replace most the frontal solver by an iterative solver for inverting the large matrix system which comes into play within the FETI-FDP2 method. Indeed, the sparse matrices when factorized are not sparse any more and can not be easily stored through compressed storage algorithms. We also exploit the *marching-on-in-anything* technique [34, 35, 36] to select an initial field estimate which is the closest from the true field value. Then, we separate the computational sub-domains into *permanent* and *non-permanent* ones in order to reduce the number of permittivity unknowns and to pre-compute as many elements as possible.

A cost function depending on the dielectric characteristics of the unknown target
is then defined and iteratively minimized. Two approaches currently prevail in the
literature, when a local minimum is searched for. In the first approach, called the
Contrast-Source inversion (CSI) method, the unknowns correspond to the current
distribution and to the permittivity. This method has been extended to various problems
[3, 37, 38, 39]. Here, we focus on the second, so-called "conventional" approach, where
the electrical field comes into play in order to satisfy the Helmholtz equation at each
iteration of the inversion process. In the inversion part, the only unknowns are the
permittivity values in each tetrahedron of the *non-permanent* subdomains.

As the dimensions of the targets are small compared to the wavelengths, their
scattered field amplitudes are very low and close to the signal-to-noise ratio [8]. It
is thus of even more interest to provide robust inversion algorithms. In particular, we
can point out a Broyden-Fletcher-Goldfarb-Shanno (BFGS) Quasi-Newton optimization
algorithm with line search [4] [40], which distinguishes itself through an approximation
of the Hessian matrix in the Newton correction step with a matrix that does not involve
the explicit computation of second order derivatives. Even if it requires additional
storage, this BFGS scheme is known for presenting a faster convergence rate than the
conjugate-gradient methods which can be seen as memory-less BFGS methods [41, 42].
To save up from additional inversion iterations and thus additional forward problems
computations, we select the BFGS method instead of the conjugate-gradient method.
But to compensate for the additional memory requirement, we use the limited-memory-
BFGS version [43] where instead of storing the full-dense approximation of the inverse
of the Hessian, we only save few vectors that represent the approximation implicitly.

Finally, taking advantage of the Lagrangian formalism and the definition of an ad-
hoc adjoint field, we calculate the gradients with respect to the permittivity. Thus,
we obtain an iterative algorithm which involves two full forward scattering problems at
each iteration of the inversion procedure.
The article is organized as follows. In Section 2, the 3D free-space inverse scattering configuration is detailed. The classical finite element method is recalled as well as the main steps of the FETI-FDP2 method. Some figures are also provided to illustrate the computational burden associated to such forward problem schemes. In Section 3, the inverse problem is stated and the expression of the associated adjoint field is recalled. As the measurements are performed in far-field, a near-field-to-far-field (NF-FF) transform [44] [45] is applied and its implication on the inversion scheme is explained. In Section 4, the details and the various steps of our proposed FETI-based inversion algorithm are described. The correct implementation of this inversion scheme is first assessed by means of synthetic noiseless dataset and by comparing the reconstructions with the ones obtained with a classical FEM algorithm. Finally, in Section 5, the inversion algorithm is confronted to measurements extracted from the 3D Fresnel database [46] [5]. Results obtained with either a classical FEM-based inversion scheme and the FETI-FDP2-based inversion scheme are compared and discussed. Conclusions follow.

## 2. Forward scattering problem

### 2.1. Electromagnetic configuration

Let us consider an isotropic three-dimensional inhomogeneous dielectric object embedded in an isotropic background with known permittivity $\varepsilon_b(\mathbf{r}) = \varepsilon_0\varepsilon_{br}(\mathbf{r})$, where $\varepsilon_0$ is the permittivity of vacuum (Figure 1). The embedding medium as well as the object under consideration is assumed to be infinite and non-magnetic ($\mu(\mathbf{r}) = \mu_0\mu_r = \mu_0$). In the following, even if the modelling capabilities enable it, we will assume that the embedding background is homogeneous, i.e., $\varepsilon_b(\mathbf{r}) = \varepsilon_b$. The time convention is in $\exp(j\omega t)$.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Geometry of the problem. (Left) A three-dimensional set of objects inside an investigation test domain $\mathcal{D}$ and a computational domain $\Omega$. (Right) Computational domain $\Omega$ when a partitioning has been performed with $N_\Omega = 10$.</figcaption></figure>

The vectorial electromagnetic field $\mathbf{E}_{sc}(\mathbf{r})$ scattered by the object is entirely determined by the characteristics of the impinging wave $\mathbf{E}_{inc}(\mathbf{r})$ and the complex
permittivity distribution $\varepsilon(\mathbf{r}) = \varepsilon_0\varepsilon_r(\mathbf{r})$ which varies from $\varepsilon_b$ only within the object. In the computational domain $\Omega$ (Figure 1), the scattered field $\mathbf{E}_{sc}$ satisfies the following equation, associated with Silver-Muller radiation conditions,

$$ \mathcal{H}(\mathbf{E}_{sc}, \varepsilon_r) = \nabla \times \left( \frac{1}{\mu_r} \nabla \times \mathbf{E}_{sc} \right) - k_0^2 \varepsilon_r(\mathbf{r}) \mathbf{E}_{sc}(\mathbf{r}) = \mathbf{J}_{sc}(\mathbf{r}) \quad \forall (\mathbf{r}) \in \Omega(1) $$

where $k_0$ corresponds to the vacuum wave-number. The secondary induced currents are defined as $\mathbf{J}_{sc}(\mathbf{r}) = k_0^2 [\varepsilon_r(\mathbf{r}) - \varepsilon_{br}] \mathbf{E}_{inc}(\mathbf{r})$.

As the sources are located in far-field, the incident field is modelled as a plane wave $\mathbf{E}_{inc}(\mathbf{r}; \mathbf{k}_s) = \mathbf{p}_s \exp(-j\mathbf{k}_s \cdot \mathbf{r})$ where $\mathbf{p}_s$ corresponds to the polarization state of the incident field and $\mathbf{k}_s$ to its wave-vector. The scattered field is collected on a probing line $\Gamma^{far}$ also located in far-field. In order to limit the computational domain $\Omega$, we take advantage of the near-field-to-far-field (NF-FF) transform [44] [45] to derive the scattered field $\mathbf{E}_{sc}^{far}(\mathbf{r}; \mathbf{k}_s)$ in far-field,

$$ \mathcal{F}(\mathbf{E}_{sc}) = \int_{\Sigma} \mathcal{Q}(\mathbf{E}_{sc}(.; \mathbf{k}_s))(\mathbf{r}; \mathbf{r}') d\mathbf{r}' = \mathbf{E}_{sc}^{far}(\mathbf{r}; \mathbf{k}_s) \quad \forall \mathbf{r} \in \Gamma^{far} \quad (2) $$

with

$$ \begin{aligned} \mathcal{Q}(\mathbf{E}_{sc}(.; \mathbf{k}_s))(\mathbf{r}; \mathbf{r}') = & j\omega\mu_0[\mathbf{n}' \times \nabla \times \mathbf{E}_{sc}(\mathbf{r}'; \mathbf{k}_s)]G_0(\mathbf{r}; \mathbf{r}') \\ & + [\mathbf{n}' \cdot \mathbf{E}_{sc}(\mathbf{r}'; \mathbf{k}_s)]\nabla' G_0(\mathbf{r}; \mathbf{r}') + [\mathbf{n}' \times \mathbf{E}_{sc}(\mathbf{r}'; \mathbf{k}_s)] \times \nabla' G_0(\mathbf{r}; \mathbf{r}') \end{aligned} \quad (3) $$

The closed surface $\Sigma$ is chosen inside $\Omega$ such that it completely encloses the scatterers and thus the investigation test domain $\mathcal{D}$. $\mathbf{n}'$ denotes the unit normal vector to $\Sigma$ at $\mathbf{r}'$ and points toward the exterior region. $G_0(\mathbf{r}; \mathbf{r}')$ is the classical 3D free-space scalar Green function.

## 2.2. Classical Finite Element method

In order to simulate electromagnetic wave scattering in various types of complex environments, a Finite Element Method (FEM) has been implemented [44] [47].

The computational domain $\Omega$ is first bounded and Radiation Boundary Conditions (RBC) are applied in order to mimic infinite free-space [44]. The computational domain is then discretized into small tetrahedrons. A first-order vector basis element scheme is employed, where the permittivity is assumed to be constant in each tetrahedron cell. The unknown components $E_{sc}$ of the vector $\mathbf{E}_{sc}$ are directly proportional to the number of edges $N_{edges}$ in $\Omega$. These values are computed by inverting the discretized version of the weak form associated to Equation (1)

$$ \mathbb{K}E_{sc} = f_{sc} \quad (4) $$

where $\mathbb{K}$ is the sum of the stiffness matrix, the mass matrix and the external boundary conditions matrices on $\partial\Omega$ and $f_{sc}$ represents the discretized version of the weak form associated to the induced current terms $\mathbf{J}_{sc}$. Contrary to the method of moments where the linear system is dense, the matrix $\mathbb{K}$ is sparse and can be easily stored thanks to compressed storage algorithms. The linear system is inverted with a frontal
method [48] specially adapted for sparse matrices. Once $E_{sc}$ is computed inside $\Omega$, it can be transformed into $E_{sc}^{far}$ thanks to the near-field to far-field transform, where $\Sigma$ corresponds to the surface of a cubic box of slightly smaller dimensions than $\Omega$.

## 2.3. FETI-FDP2 Finite Element method

The finite element method previously described has been further improved to address large scale models while preserving the versatility of the FEM method. To that end, an updated version of the Finite Element Tearing and Interconnecting-Dual Primal Electromagnetic method (FETI-DPEM) [25], named FETI-FDP2, has been developed and thoroughly tested [32] [33].

The general principle of the FETI methods is to divide the entire computational domain $\Omega = \bigcup_{i=1}^{N_{\Omega}} \Omega^i$ into smaller non-overlapping subdomains with $\Omega^i \cap \Omega^j = \emptyset, \forall i \neq j$. An arbitrary mesh partitioner [49] yielding strongly irregular surfaces with varying curvatures is employed, leading to subdomains $\Omega^i$ with no specific planar interfaces (Figure 1). In each of these subdomains $\Omega^i$, local linear systems can now more easily be inverted. The discretized version of the weak form associated to Equation (1) then reads

$$ K^i E_{sc}^i = f_{sc}^i - D^i \Lambda_{sc}^i \quad i = 1, \dots, N_{\Omega} \qquad (5) $$

where $K^i$ is the sum of the stiffness matrix, the mass matrix and the external and internal boundary conditions matrices on $\partial\Omega$ and $\partial\Omega^i$. $E_{sc}^i$ comprises the unknown components of the vector $E_{sc}$ in $\Omega^i$ and is directly proportional to the number of edges $N_{edges}^i$ in $\Omega^i$. $f_{sc}^i$ represents the discretized version of the weak form associated to the induced current terms $J_{sc}$ for elements in $\Omega^i$. $D^i$ is a Boolean matrix which extracts only the edges of $\Omega^i$ which are on $\partial\Omega^i$. As for the dual Lagrange multipliers $\Lambda_{sc}^i$, they are involved in the unknown boundary conditions which are imposed at the internal interfaces between adjacent subdomains. These adjacent subdomains are glued at their common interfaces thanks to appropriate boundary conditions and constraints imposed to the field and its derivatives, leading to the following matrix system

$$ \Lambda_{sc}^{i \to j} + \Lambda_{sc}^{j \to i} = -W^{i \leftrightarrow j} E_{sc}^{j \to i} \quad i, j = 1, \dots, N_{\Omega} \text{ with } i \neq j \qquad (6) $$

where the matrix $W^{i \leftrightarrow j}$ directly translates the boundary conditions on $\partial\Omega^i \cap \partial\Omega^j$. $E_{sc}^{j \to i}$ (resp. $E_{sc}^{i \to j}$) is the subset of $E_{sc}^i$ associated to the edges of $\Omega^i$ (resp. $\Omega^j$) which belong to the interface $\partial\Omega^i \cap \partial\Omega^j$.

By combining the previous linear systems, we end up with the so-called Interface Problem,

$$ F_{sc} \Lambda_{sc} = d_{sc} \qquad (7) $$

where the vector $\Lambda_{sc} = [\Lambda_{sc}^1, \dots, \Lambda_{sc}^{N_{\Omega}}]^T$ contains all the Lagrange multipliers. The explicit definition of the matrix $F_{sc}$ and the right-hand-side term $d_{sc}$ can be found in [33].

For a given illumination $(p_s, k_s)$ and a given permittivity distribution $\epsilon_r$, the forward scattering problem requires first to solve the interface problem provided
in Equation (7). Then, the solution inside each subdomain $\Omega^i$ can be evaluated independently by solving Equation (5). When the size of the computational domain $\Omega$ increases with respect to the wavelength, the resolution of the various linear systems requires more and more memory storage and computational resources. Nevertheless, it is possible by increasing $N_\Omega$ to limit the number of discretization cells in each subdomain $\Omega^i$ (Table 1).

<table><caption>**Table 1.** Number of elements when $\Omega$ is a parallepipedic domain discretized with a fixed mesh step and separated into $N_\Omega = 25$ non-overlapping subdomains.</caption><thead><tr><th>$f$ (GHz)</th><th>Size of $\Omega$ in $\lambda_b^3$</th><th>$N_{edges}$</th><th>Average $N_{edges}^i$</th><th>$N_\Lambda$</th></tr></thead><tbody><tr><td>4</td><td>3.3 × 3.3 × 2.7</td><td>168 743</td><td>6 749</td><td>41 350</td></tr><tr><td>6</td><td>4.0 × 4.0 × 3.0</td><td>291 182</td><td>11 647</td><td>59 066</td></tr><tr><td>8</td><td>4.7 × 4.7 × 3.3</td><td>461 155</td><td>18 446</td><td>83 836</td></tr></tbody></table>

With such a strategy, Equation (5) can still be inverted with a direct frontal method specially adapted for sparse matrices [48]. On the contrary, the size of the Interface Problem increases with the frequency and the number of sub-domains, as it is directly linked to the number of boundary conditions and constraints applied on the internal interfaces $\partial\Omega^i \cap \partial\Omega^j$. Figure 2 shows an example of the requirement in terms of memory consumption and time when the FETI-FDP2 method is employed. The FEM case is more or less corresponding to the case when $N_\Omega = 1$. It is obvious that the memory requirement for the FETI-FDP2 method is much less than the one associated with the FEM method, even if there is a trade-off to find between a large $N_\Omega$ (a large collection of little sub-domains easily invertible (Equation 5) combined with a large Interface Problem (Equation 7)) or a small $N_\Omega$ (a small collection of large sub-domains combined with a small Interface Problem).

Instead of a frontal solver for solving Equation (7), we prefer an iterative solver, such as the Generalized Minimal Residual Method (GMRES) algorithm [50], in order to limit the memory storage even if we lose on the computational time. Moreover, for the time being, all the implementation is performed in a sequential way. We thus do not expect a real gain on the computation time. This gain will certainly be achieved once the entire algorithm will be parallelized as the resolution of the local sub-problems can be easily distributed in an independent way on several processors. As any iterative algorithm, a sequence of iterates $\{\Lambda_{sc,s}^{(0)}(\epsilon_r), \dots, \Lambda_{sc,s}^{(g)}(\epsilon_r), \dots, \Lambda_{sc,s}^{(G_s)}(\epsilon_r)\}$ is generated with the GMRES algorithm, for a given initial guess $\Lambda_{sc,s}^{(0)}(\epsilon_r)$ and a stopping criterion $\eta$. $G_s$ corresponds to the final number of GMRES iterations. Once the field is computed everywhere in $\Omega$, the NF-FF transform provides $E_{sc}^{far}$.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Evolution with respect to the number of sub-domains $N_\Omega$ of the memory consumption of the FETI-FDP2 method and the time requested to perform a $F_{sc}$ matrix-vector product. The computational domain $\Omega$ is an almost $90\lambda^3$ parallelepiped, discretized with 15 points per wavelength leading to $1.5 \times 10^6$ edges and $1.4 \times 10^6$ tetrahedrons.</figcaption></figure>

### 3. Inverse scattering problem

The scattered field measurements $E^{mes}$ are acquired on a probing line $\Gamma^{far}$ far from the target position. The receiving antenna, located at $\mathbf{r}_r$, is assumed to be linearly polarized along $\mathbf{q}_r$, $\mathbf{q}_r$ being a unit vector. The field actually measured for each receiving point is thus $E_{sc}^{far}(\mathbf{r}_r; \mathbf{k}_s) \cdot \mathbf{q}_r$. Two possibilities occur to compare the simulated field and the measured field. Either the far-field measurement is back-propagated from $\Gamma^{far}$ into the computational domain $\Omega$ with a far-field-to-near-field (FF-NF) transform, or the simulated field is extracted from $\Omega$ and projected onto $\Gamma^{far}$ thanks to the NF-FF transform. We have chosen the second option for several reasons: (i) the integral operator defined in Equation (7) is known to be ill-posed and suffers some real issues when one tries to invert it, taking into account that the number of receivers is limited and the measurements are unfortunately noisy, (ii) the scattered field is a vectorial field and to be able to fully capture it would require to measure all the polarization cases $\mathbf{q}_r$ which might not be available in practical situations. The misfit criterion is thus performed on the probing line and corresponds, for a given permittivity distribution $\varepsilon_r$, to

$$ \mathcal{J}(\varepsilon_r) = \sum_{s=1}^{N_{src}} \sum_{r=1}^{N_{rec}} w_{s,r} |E_{s,r}^{mes} - E_{sc}^{far}(\mathbf{r}_r; \mathbf{k}_s) \cdot \mathbf{q}_r|^2 \quad (8) $$

where $w_{s,r}$ corresponds to weighting coefficients enabling to balance each measurement point with respect to its own measurement uncertainty [51]. $N_{rec}$ (resp. $N_{src}$) corresponds to the maximum number of combinations of $(\mathbf{q}_r, \mathbf{r}_r)$ (resp. $(\mathbf{p}_s, \mathbf{k}_s)$) used in the experiment.

In our inverse problem scheme, the aim is to recover the permittivity distribution which minimizes the above misfit criterion while fulfilling the various equations (1) (2)
(3). In order to take into account all constraints, the following Lagrangian functional is introduced

$$ \mathcal{L}(\mathbf{E}_{\text{sc}}^{\text{far}}, \mathbf{E}_{\text{sc}}, \varepsilon_r, \mathbf{U}, \mathbf{P}) = \mathcal{J}(\varepsilon_r) + \sum_{s=1}^{N_{\text{src}}} \text{Re} \langle \mathbf{U} | \mathcal{F}(\mathbf{E}_{\text{sc}}(\cdot; \mathbf{k}_s)) - \mathbf{E}_{\text{sc}}^{\text{far}}(\cdot; \mathbf{k}_s) \rangle_{\Gamma^{\text{far}}} + \sum_{s=1}^{N_{\text{src}}} \text{Re} \langle \mathbf{P} | \mathcal{H}(\mathbf{E}_{\text{sc}}(\cdot; \mathbf{k}_s), \varepsilon_r) - \mathbf{J}_{\text{sc}}(\cdot; \mathbf{k}_s) \rangle_{\Omega} \quad (9) $$

where $\mathbf{U}$ and $\mathbf{P}$ are Lagrange multipliers. $\langle \cdot | \cdot \rangle_{\Gamma^{\text{far}}}$ corresponds to the $L_2$ inner product on $\Gamma^{\text{far}}$ and $\langle \cdot | \cdot \rangle_{\Omega}$ to the $L_2$ inner product on $\Omega$.

As we select an iterative optimization scheme based on a gradient-descent direction algorithm, we need to compute the gradient $\nabla_{\varepsilon_r} \mathcal{J}(\varepsilon_r)$ of the misfit criterion. This quantity can be directly evaluated when the Karush-Kuhn-Tucker (KKT) conditions [52] are met, which enables to find the saddle-point of $\mathcal{L}$. At this saddle-point, Equations (1) and (2) are satisfied. At this saddle-point, new equations arise, linking the scattered field and the Lagrange multipliers to provide the definition of an adjoint field $\mathbf{P}_{\text{adj}}(\cdot; \mathbf{k}_s)$ for each illumination pair $(\mathbf{p}_s, \mathbf{k}_s)$ such that

$$ \mathcal{H}(\mathbf{P}_{\text{adj}}(\cdot; \mathbf{k}_s), \varepsilon_r) = \mathbf{J}_{\text{adj}}(\mathbf{r}; \mathbf{k}_s) \quad (10) $$

with

$$ \mathbf{J}_{\text{adj}}(\mathbf{r}; \mathbf{k}_s) = -\mathcal{F}^{\dagger,*} \left( 2 \sum_{r=1}^{N_{\text{rec}}} w_{s,r} [E_{s,r}^{\text{mes}} - \mathbf{E}_{\text{sc}}^{\text{far}}(\mathbf{r}_r; \mathbf{k}_s) \cdot \mathbf{q}_r]^* \delta(\mathbf{r} - \mathbf{r}_r) \mathbf{q}_r \right) \quad (11) $$

The operator $\mathcal{F}^\dagger$ corresponds to the adjoint operator of $\mathcal{F}$, such that $\langle \mathbf{u} | \mathcal{F}(\mathbf{v}) \rangle_{\Gamma^{\text{far}}} = \langle \mathcal{F}^\dagger(\mathbf{u}) | \mathbf{v} \rangle_{\Sigma}$. The notation $\mathbf{u}^*$ corresponds to the complex conjugate of $\mathbf{u}$. The adjoint field thus satisfies a similar equation as $\mathbf{E}_{\text{sc}}$ (see Equation (1)) apart from the right-hand-side term, where receivers now act as sources and emit the discrepancy existing between the simulated and measured fields.

Unfortunately, the adjoint operator $\mathcal{F}^\dagger$ is not straightforward to retrieve and we are more or less facing the same issues as when one tries to search for an effective FF-NF transform. Instead of computing $\mathcal{F}^\dagger$, we follow an other strategy. As the receivers are in far-field, we construct incident plane waves arising from the direction $\mathbf{k}_r$. The adjoint source field term is then computed by summing up these plane waves with appropriate weighting factors

$$ \mathbf{J}_{\text{adj}}(\mathbf{r}; \mathbf{k}_s) = -2 \sum_{r=1}^{N_{\text{rec}}} w_{s,r} [E_{s,r}^{\text{mes}} - \mathbf{E}_{\text{sc}}^{\text{far}}(\mathbf{r}_r; \mathbf{k}_s) \cdot \mathbf{q}_r] \mathbf{q}_r \exp(-j\mathbf{k}_r \cdot \mathbf{r}) \quad (12) $$

The adjoint field is then computed as explained in Section 2.2 or in Section 2.3. When the FETI approach is followed, the associated adjoint interface problem (Equation 7) is solved with an iterative GMRES algorithm leading to an other sequence of Lagrange multipliers $\{\Lambda_{\text{adj},s}^{(0)}(\varepsilon_r), \dots, \Lambda_{\text{adj},s}^{(h)}(\varepsilon_r), \dots, \Lambda_{\text{adj},s}^{(H_s)}(\varepsilon_r)\}$, where $H_s$ denotes the final number of GMRES iteration. Once the scattered field and the adjoint field are computed, the gradient of $\mathcal{J}$ is easily derived as

$$ \nabla_{\varepsilon_r} \mathcal{J}(\varepsilon_r)(\mathbf{r}) = -k_0^2 \sum_{s=1}^{N_{\text{src}}} \mathbf{P}_{\text{adj}}^*(\mathbf{r}; \mathbf{k}_s) \cdot [\mathbf{E}_{\text{inc}}(\mathbf{r}; \mathbf{k}_s) + \mathbf{E}_{\text{sc}}(\mathbf{r}; \mathbf{k}_s)]^* \quad (13) $$
Since a solution of the optimization problem has to be sought numerically, a parameter representation of the complex permittivity function $\epsilon_r(\mathbf{r})$ is needed to obtain a finite number of optimization variables. We thus search for the permittivity values in each of the mesh tetrahedrons. We could have used a coarser grid for representing the permittivity unknown or a non-local basis [15] but we select for the time being the simplest approach. Nevertheless, as the number of measurements is limited and much smaller than the number of tetrahedrons, we limit the search to a subdomain $\mathcal{D}$ of the computational domain $\Omega$. As in any iterative scheme, a sequence of iterates $\{\epsilon_r^{(0)}, \dots, \epsilon_r^{(v)}, \dots, \epsilon_r^{(V)}\}$ is generated for a given initial guess $\epsilon_r^{(0)}$. These permittivity maps may only differ from $\epsilon_b$ within the investigation domain $\mathcal{D}$.

The sequence of permittivity maps is constructed thanks to a Quasi-Newton algorithm

$$ \epsilon_r^{(v+1)} = \epsilon_r^{(v)} + \alpha^{(v)} \mathbb{B} \nabla_{\epsilon_r} \mathcal{J}(\epsilon_r^{(v)}) \quad (14) $$

where $\mathbb{B}$ is an approximation of the inverse of the Hessian $[\nabla_{\epsilon_r}^2 \mathcal{J}(\epsilon_r)]$ and $\alpha^{(v)}$ is the step length. Among all the different approximations for computing the matrix $\mathbb{B}$, we follow the classical Broyden-Fletcher-Goldfarb-Shanno (BFGS) method [40] as it is known to converge faster than conjugate gradient schemes even if it requires additional storage. As we are dealing with 3D large-scale problems, the memory consumption aspect can not be neglected, especially as the Hessian matrix has dimension $N_{tets} \in \mathcal{D} \times N_{tets} \in \mathcal{D}$, where $N_{tets} \in \mathcal{D}$ is the number of tetrahedrons inside $\mathcal{D}$. We have thus adopted the limited-memory version of the BGFS algorithm [43].

The inversion iterations are stopped if the cost function $\mathcal{J}(\epsilon_r)$ is becoming lower than a stopping criterion $\xi$, or if the cost function stagnates or if a maximum number of iterations (here 100) is reached.

## 4. Efficient FETI-based inversion algorithm

The aforementioned inversion algorithm requires to solve two forward problems at every iteration step of the optimization process. These are for the direct and the adjoint problems. Several inversion schemes based on classical Finite Element methods have already been implemented. But high-frequency electromagnetic scattering problems call for fine meshes, and therefore lead to large-scale systems of equations. For such problems, solving the forward problem with a direct method entails memory and CPU requirements that rapidly overwhelms even the largest resources that are currently available. Thus, in order to make this process more flexible, we have implemented the FETI-FDP2 method previously proposed and discussed. Our aim is not just to plug the FETI-FDP2 method to replace the FEM method, but it is also to show how a Domain Decomposition Finite Element method can be efficiently combined with an iterative optimization algorithm and to point out the different strategies that we have followed in order to limit the computational burden and reduce the memory storage.
### 4.1. Comparison criteria

In order to quantitatively appraise the effectiveness of the implementation strategy, we compare the results obtained with the FETI-FDP2 approach with the ones that are retrieved when combining the inversion algorithm with a classical FEM method. Two quantitative criteria are thus introduced for that purpose. The first one relies on the convergence of the cost function $J$, which shows the correct behaviour of the two inversion algorithms and the number of iterations required by each scheme to converge. The second criterion corresponds to the discrepancy between the permittivity map computed with the FEM forward problem and the FETI-FDP2 forward problem,

$$ \mathcal{N}(\varepsilon_r^{(v)}) = \|\varepsilon_{r,FEM}^{(v)} - \varepsilon_{r,FETI}^{(v)}\|_{\Omega}^2 \quad (15) $$

where $\|\cdot\|_{\Omega}^2$ is the $L_2$ norm on $\Omega$. This criterion is either computed all along the iterations or for the last iteration of each method.

### 4.2. Synthetic configuration

A specific synthetic configuration is constructed in order to compare the various schemes. Additional inversion results are also presented based on experimental datasets in Section 5. This synthetic configuration corresponds to the TWO CUBES object from the Fresnel database [46]. The choice of this target is dictated by its size and its simplicity. The TWO CUBES object consists of two identical cubes with permittivity $\varepsilon_r = 2.25$ and a side length of 0.025 m. They are organized such that they share a common corner point. The positions, the number and the polarizations of the emitters and receivers are exactly the ones described in the Fresnel database and are situated at approximately 2 m from the target. The background medium is air. We choose the operating frequency equal to 4 GHz ($\lambda_b = 0.075$ m). The computational domain $\Omega$ is a cube with a side length of 0.23 m ($\sim 3.1\lambda_b$). The investigation domain $\mathcal{D}$ is a sphere with a radius of 0.05 m centred in $\Omega$. The finite element discretization is set to 10 points per wavelength, leading to $N_{tets} \in \mathcal{D} = 15197$ in $\mathcal{D}$. The initial guess $\varepsilon_r^{(0)}$ corresponds to $\varepsilon_{rb}$ apart from a little cube with a 0.01 m side length, with one corner at the center of the set-up and a permittivity equal to 1.1. All the tests in this section are performed on synthetic noiseless scattered fields. We are thus typically in a full-inverse crime situation, even if the FEM forward problem is employed to compute the simulated fields.

### 4.3. Permanent and non-permanent information

The inversion scheme and the FETI-FDP2 forward problem are fully iterative techniques. If we want to efficiently implement the FETI-FDP2 method, we have to conditionally divide it into two parts: the *permanent* and *non-permanent* parts. The subdomains $\Omega^i$ are also classified as *permanent* and *non-permanent* subdomains. Indeed, the *permanent* subdomains keep the same permittivity values and the same geometrical features for all the iterative scheme as they are not intersecting with $\mathcal{D}$.
The *permanent* step thus handles all the computational steps related to information that do not change during the calculations. This corresponds at first to all the geometrical features, such as the mesh generation (lists of nodes, edges and elements), the mesh geometrical information extraction (matrices such as $D^i$), the electromagnetic configuration (positions and orientations of the emitters and receivers, plane waves illuminations, etc). For the *permanent* subdomains, the matrices $K^i$ are computed once for all as well as some parts of the $F_{sc}$ matrix. This *permanent* step is just performed a single time at the beginning of the inversion iterative scheme.

The *non-permanent* step handles the parts that depend on the updated values of $\epsilon_r^{(v)}$. The corresponding matrices $K^i$, the *non-permanent* part of the $F_{sc}$ matrix as well as the right-hand-side terms for both the direct and the adjoint problems are thus updated at each inversion iterative step.

### 4.4. Forward problem initialization

For each iteration $v$ of the inversion process and for each illumination ($\mathbf{p}_s, \mathbf{k}_s$), the initial values $\Lambda_{sc,s}^{(0)}(\epsilon_r^{(v)})$ and $\Lambda_{adj,s}^{(0)}(\epsilon_r^{(v)})$ must now be properly defined. We investigate three different approaches based on the so-called *marching-on-in-anything* technique [34, 35, 36] which has already been shown to accelerate the convergence procedures. For $\bullet \in \{\text{sc, adj}\}$, these approaches are based on

(i) a *blank* initialization: $\Lambda_{\bullet,s}^{(0)}(\epsilon_r^{(v)}) = 0$,

(ii) a *marching-on-in-anything* permittivity initialization: $\Lambda_{\bullet,s}^{(0)}(\epsilon_r^{(v)}) = \Lambda_{\bullet,s}^{(K_s)}(\epsilon_r^{(v-1)})$,

(iii) a *marching-on-in-anything* source initialization: $\Lambda_{\bullet,s}^{(0)}(\epsilon_r^{(v)}) = \Lambda_{\bullet,s-1}^{(K_s)}(\epsilon_r^{(v)})$.

where $K_s$ is either $G_s$ if the forward problem is solved for computing the scattered field, or $H_s$ if it is for the adjoint field. In the original *marching-on-in-anything* technique, the initialization is performed by combining several previous solutions. In our case, as the sources are far apart from each other, the previous computed fields will be far from the one which is currently sought-after. There is little chance that the process will be accelerated by combining more than one previous solution. In order to provide a fair comparison with the other initialization possibilities, we have decided to only use the last field and not to combine several of them.

In all cases, the GMRES stopping criterion $\eta$ is set to $10^{-2}$. The results associated to the various strategies are presented in Figures 3 and 4. The FETI-based inversion algorithm often stops earlier than the FEM-one. It can be due to the value selected for the GMRES stopping criterion $\eta$, which is quite large with respect to the machine precision. This choice has been done in order to speed up the GMRES iterative algorithm, but it implies that the computed field will be different from the one computed with the FEM based on a direct solver. As the scattered field and the adjoint field are approximated, the gradient is as well approximated. When the gradient is small, that is, at the end of the inversion iterative process, these approximations come into play. Tiny variations of the permittivity profile will not be fully reflected in the computed
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Evolution of $J(\epsilon_r)$ (in dB) for the TWOCUBES target with noiseless fields simulated at 4 GHz. The FETI-FDP2 method is computed with various GMRES initialization strategies [blank (□), marching-on-in-anything permittivity (○) and marching-on-in-anything source (Δ)] or with FEM (—).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** (Left) Geometrical configuration of the TWOCUBES. A line crossing the diagonals of the two cubes as well as the shared corner is represented. (Right) Profile cut along this crossing line of the final reconstructed real parts of $\epsilon_{r,FEM}$ and $\epsilon_{r,FETI}$ with different GMRES initialization strategies, obtained from simulated noiseless scattered fields at 4 GHz.</figcaption></figure>

field, due to the choice of $\eta$. Thus, the inversion process will stop earlier than the FEM-inversion scheme as the cost function will reach a plateau. In spite of this, the solutions obtained with the various initialisation strategies are very similar to the one obtained with a FEM-based algorithm. It is thus difficult to derive the best initialisation so far if we only compare the reconstructed permittivity maps.

Let us now consider the number of iterations of the FETI-FDP2 method at each iteration step of the inversion algorithm, that is the number $G_s$ (resp. $H_s$) associated to the direct (resp. adjoint) problem. To simplify, only their average values $G$ and $H$ with respect to the number of sources are plotted in Figure 5, knowing that $G_s$ and $H_s$ vary only from up to $\pm 1$ iteration from these average values.

It appears that the most efficient initialization strategy differs with the problem to solve: the *marching-on-in-anything* permittivity is more efficient for the forward problem while the *marching-on-in-anything* source is more efficient for the adjoint problem. This might be due to the fact that, for the direct problem, the sources are located far apart and the fields change from one source to another more than from one permittivity map
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Average number versus the number of iterations of the inversion algorithm of the GMRES iterations obtained either with the *blank* initialization (—), the *marching-on-in-anything* permittivity initialization (□), or the *marching-on-in-anything* source initialization (○). (Left) *G* values for the direct field computation and (Right) *H* values for the adjoint field computation.</figcaption></figure>

to another. Of course, such results should be mitigated if more sources were present. On the contrary, for the adjoint problem, the incident adjoint field is generated by the receivers emitting all together the discrepancy between the measured and simulated scattered field. Thus, from one source to the other one, there are no changes in the emitters location, only the amplitudes slightly differ. This seems the reason why the *marching-on-in-anything* source initialization works better in this case.

Given the above convergence results, we have selected the following strategy: a *marching-on-in-anything* source initialization for the adjoint problem and a *marching-on-in-anything* permittivity initialization for the forward problem, with the possibility in this latter case to switch to a *blank* initialization if the relative error of the first GMRES iteration is larger than 1.0.

### 4.5. Forward problem stopping criterion

The stopping criterion $\eta$ of the FETI-FDP2 method must also be properly chosen. In order to obtain the same solution as the FEM one, which uses a direct solver, one should set $\eta$ to the computer precision. Unfortunately, as practice shows, the multi-sources calculation is going to be tremendously slow. Moreover, taking a very small $\eta$ does not guarantee the inversion scheme to converge towards the true permittivity distribution. Finally, the forward and adjoint problem computations should be performed taking into account that the measured fields might be noisy and an "infinite" precision is not systematically required. Increasing $\eta$ should mainly influence the gradient computation.

Different stopping criteria are used. The convergence behaviour of $\mathcal{J}$ shown in Figure 6 indicates that for all these stopping criteria, the cost function follows more or less the same evolution, even if it tends to stop earlier when $\eta$ is above $10^{-2}$. In all cases, the inversion algorithm converges to a correct solution (see Figure 7) regardless of the stopping criterion. Of course, even in an inverse crime situation, the final solution is not identical to the true solution as the number of measurement points is much smaller than the number of unknowns and the problem is ill-posed.

As the cost function behaviour only differs after a certain number of iterations, we
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Evolution of $\mathcal{J}(\varepsilon_r)$ (in dB) for the TWO Cubes target with noiseless fields simulated at 4 GHz. The FETI-FDP2 method is computed with various GMRES stopping criteria $\eta$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Profile cut along the crossing line shown in Figure 4(a) of the final reconstructed real parts of $\varepsilon_{r,\text{FEM}}$ and $\varepsilon_{r,\text{FETI}}$ with different GMRES stopping criteria $\eta$, obtained from simulated noiseless scattered fields at 4 GHz.</figcaption></figure>

have thus adopted a strategy for dynamically selecting $\eta$ according to the value of $\mathcal{J}(\varepsilon_r)$. When $\mathcal{J}(\varepsilon_r^{(v)})$ is above $10^{-2}$, we set $\eta = \mathcal{J}(\varepsilon_r^{(v)})$. Once $\mathcal{J}(\varepsilon_r^{(v)})$ gets below $10^{-2}$, we fix $\eta$ to $10^{-2}$. With such a strategy, applied both to the computation of the forward problem and the adjoint problem, we reduce the number of GMRES iterations at the beginning of the inversion process while still ensuring a correct computation of the gradients at the end of the inversion algorithm.

**5. Inversion of experimental data**

The aforementioned inversion algorithm is now confronted to measurements extracted from the 3D Fresnel database [46]. This database involves a set of homogeneous targets, measured in free-space and in far-field, for various positions of sources and receivers, and for various frequencies. For the two inversion algorithms (FEM and FETI-based), only the results obtained at 4 GHz are presented herein. Indeed, we want to compare the behaviour of the two inversion algorithms and thus to limit the size of the computational
domain in order to make it still manageable with a classical FEM forward solver. The investigation test-domain $\mathcal{D}$ is a sphere centred at (0, 0, 0.025) m with a diameter of 0.16 m (~ 2.1 $\lambda_b$). This test-domain has been selected as it enables to comprise all the targets from the database. The computational domain $\Omega$ is a cubic domain, centred at (0, 0, 0.0) m with a side-length of 0.26 m (~ 3.5 $\lambda_b$). The initial guess $\varepsilon_r^{(0)}$ differs from $\varepsilon_{rb}$ only within a little cube of 0.01 m side length where the permittivity is set to 1.1. The discretization mesh parameters are described in Table 2. All computations have

<table><caption>**Table 2.** Discretization information when $\Omega$ is a cubic domain of 0.26 m side-length, discretized with 10 points per wavelength, at 4 GHz.</caption><thead><tr><th>$N_{\text{nodes}}$</th><th>$N_{\text{edges}}$</th><th>$N_{\text{tets}} \in \Omega$</th><th>$N_{\text{tets}} \in \mathcal{D}$</th><th>$N_{\Lambda}$</th><th>$N_{\Omega}$</th><th>$N_{\Omega}^{\text{permanent}}$</th></tr></thead><tbody><tr><td>51 050</td><td>362 280</td><td>302 244</td><td>77 539</td><td>89 308</td><td>38</td><td>12</td></tr></tbody></table>

been performed on an Intel(R) Xeon(R) CPU X5570 @ 2.93GHz, with 48 GB of RAM, with no parallel programming specificities. The FEM-based inversion algorithm requires 15.3 GB of operative memory to reconstruct the objects with the mesh described above. The FETI-FDP2 scheme takes only 5.7 GB. As expected (see Section 2.3), the FETI-FDP2-inversion scheme requires much less memory than the FEM-inversion scheme as there is no need to invert and store a large $N_{\text{edges}} \times N_{\text{edges}}$ matrix. Instead, the largest matrix system is associated to the Interface problem and consists of a $N_{\Lambda} \times N_{\Lambda}$ matrix. As the time requirement highly depends on the number of GMRES iterations, we will discuss it for each object under reconstruction.

## 5.1. Target with cubes

The permittivity map reconstructed for the TWOCUBES target from the measured scattered fields is presented in Figure 8. The corresponding cost function is shown in Figure 9. A specific cut of the permittivity profile is visible in Figure 10.

The behaviour of the cost function is similar for the two algorithms, as well as the associated permittivity reconstructed maps. Even if all the objects from the Fresnel database have no imaginary part, our inversion algorithm is based on the independent search of both the real part and the imaginary part. As no a-priori constraints are set on the unknowns, it may happen that the real part and the imaginary part of the permittivity take non-physical values. Nevertheless, the reconstructed values are relevant and in particular, the imaginary part is close to 0. The average time for one inversion iteration is equal to 1030 sec for the FETI-FDP2-based scheme against 450 sec for the FEM solution. As no parallelization procedure is followed, the FETI-FDP2-based scheme is slower than the FEM-based one. Indeed, each forward/adjoint problem requires the inversion of $N_{\Omega}$ small linear systems sequentially.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** (Left) Geometry of the TWOcUBES target compared to the investigation domain $\mathcal{D}$ and the computational domain $\Omega$. A crossing black line is also presented. (Right) Vertical and horizontal slices, performed at the cube centres, of the real part of the permittivity map reconstructed with the FETI-FDP2 inversion scheme from measurements acquired at 4 GHz.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Evolution of the least squares data fit cost function $J(\epsilon_r)$ (in dB) over the iterations for the TWOcUBES target, using the FEM- and FETI-based inversion algorithms with measured fields at 4 GHz.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** (Left) Real and (Right) Imaginary profile cut of the permittivity reconstructed with the FEM (—) and FETI-based (---) algorithm, from scattered fields measured at 4 GHz. The actual profile (......) along the line crossing the diagonals of the two cubes as shown in Figure (8) is also plotted.</figcaption></figure>

## 5.2. Targets with spheres

Other objects from the Fresnel database are mainly constructed with spheres. This is the case of the TWOSPHERES, the CUBESPHERES and the MYSTER targets. Their
associated reconstructions are shown in Figures 11, 12 and 13.

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Reconstruction of the TWO SPHERES target obtained with the FETI-based inversion algorithm from fields measured at 4 GHz. (Left) 3D iso-surface of the reconstructed permittivity map. (Right) Vertical and horizontal slices of the real part of the permittivity map. The slices are passing through the centres of the spheres. The semi-transparent object inside represents the actual boundaries of the target.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Reconstruction of the CUBES PHERES target obtained with the FETI-based inversion algorithm from fields measured at 4 GHz. (Left) 3D iso-surface of the reconstructed permittivity map. (Right) Vertical and horizontal cuts of the real part of the permittivity map. The slices are passing through the centre of the cube of spheres. The semi-transparent object inside represents the actual boundaries of the target.</figcaption></figure>

The localisation and the range of permittivity of the various targets are correctly recognized. The profiles do not faithfully follow the actual profiles but this is expected due to the presence of the noise, to the fact that the scattering operator is a low-pass filter and that there are much more unknown parameters than degrees of freedom available within the measured fields [53]. When the size of the individual spheres is
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Reconstruction of the MYSTER target obtained with the FETI-based inversion algorithm from fields measured at 4 GHz. (Left) 3D iso-surface of the reconstructed permittivity map. The semi-transparent object inside represents the actual boundaries of the target. (Right) Slices of the real part of the relative permittivity at five different altitudes (top: $z = 0$ mm and $z = 13.74$ mm, center: $z = 17.985$ mm, down: $z = 22.23$ mm and $z = 35.97$ mm). The actual boundaries are also superimposed.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14.** (Top) Evolution of the cost function $J(\epsilon_r)$ (in dB) along the iterations for the FEM- and FETI-based inversion algorithm with measured fields at 4 GHz. (Bottom) Average number of GMRES iterations for the direct field and adjoint field computation. (Left) CUBESPHERES target and (Right) MYSTER target.</figcaption></figure>

getting too small with respect to the wavelength (for example for the CUBESPHERE, each sphere has a diameter approximatively equal to $0.2 \lambda_b$), the inversion algorithm fails to properly separate them. One should also note that we set $w_{s,r} = 1$ and thus the measurement uncertainty associated to each measurement point is not taken into account, which degrades the reconstructed maps [54]. As the spatial coverage of each individual sphere is larger than the effective one, the reconstructed permittivity value is lower than expected in order to provide in the end the same scattered power.
In terms of convergence, the FEM and FETI-based inversion algorithms are behaving almost identically (Figure 14). The FEM-based inversion converges slightly better than the FETI-FDP2-based inversion as the linear system is inverted with a direct solver instead of an iterative solver, and provides results computed up to the machine precision instead of the $\eta$ precision selected for the GMRES algorithm. The average numbers G and H of GMRES iterations are always quite low (below 5) leading to an average time computation per inversion iteration of 2650 sec for the CUBESPHERES and to 1980 sec for the MYSTER target when the FETI-FDP2 algorithm is used (Figure 14), and to 450 sec when the FEM-based inversion is used. The average time per iteration for the FEM-based inversion does not change as a direct solver is used to invert (Equation 4). Instead, the time associated to the FETI-DP2-based inversion varies due to the GMRES iterative solver. This computational time is slightly larger than for the TWOCUBES target as the targets with spheres have overall dimensions larger than the TWOCUBES ones. This computational time can be drastically reduced thanks to a parallel implementation of the entire inversion algorithm.

## 6. Conclusion

In the current article, we have shown how we have efficiently coupled a finite-element forward problem with a Quasi-Newton optimization algorithm to solve a three-dimensional inverse scattering problem. Taking advantage of the Lagrangian formalism and the definition of an *ad-hoc* adjoint field, we computed the gradient of the data misfit criterion and incorporated it into the inversion scheme. In order to render this process less depending on the memory, we applied a domain-decomposition technique to solve the various forward problems. The incorporation of the FETI-FDP2 method has been performed not in a brute-force way but in a specific manner such as to reduce as much as possible the memory and the computational time associated to each part of the entire inversion algorithm. The FETI-based inversion algorithm has been then successfully tested on various objects from the Fresnel database. In the near future, this inversion scheme will be extensively exploited to invert measurements acquired in aspect-limited configurations [55]. Additional *a-priori* information will be introduced based on the target properties (homogeneity [56], limited spatial support [18], etc) in order to improve the spatial resolution of the reconstructed permittivity map. Finally, large-scale configurations will also be studied. To that end, a parallelization of the entire scheme will be performed.

## References

[1] M. Benedetti, M. Donelli, and A. Massa. Multicrack detection in two-dimensional structures by means of GA-based strategies. *IEEE Transactions on Antennas and Propagation*, 55(1):205–215, January 2007.

[2] P. Paul, F.G. Guimaraes, and J.P. Webb. Reducing the computational cost of inverse scattering
problems with evolutionary algorithms. *IEEE Transactions on Magnetics*, 45(3):1514–1517, March 2009.

[3] P.M. van den Berg and R.E. Kleinman. A contrast source inversion method. *Inverse Problems*, 13(6):1607, December 1997.

[4] A. Franchois and A. G. Tijhuis. A quasi-Newton reconstruction algorithm for a complex microwave imaging scanner environment. *Radio Science*, 38(2):8011, 2003.

[5] A. Litman and L. Crocco. Testing inversion algorithms against experimental data: 3D targets. *Inverse Problems*, 25(2):020201, February 2009.

[6] S. R. Arridge. Optical tomography in medical imaging. *Inverse Problems*, 15(2):R41, April 1999.

[7] G. Marin, C. Guerin, S. Baillet, L. Garnero, and G. Meunier. Influence of skull anisotropy for the forward and inverse problem in EEG: simulation studies using FEM on realistic head models. *Human Brain Mapping*, 6(4):250–269, 1998.

[8] M. Seger, G. Fischer, R. Modre, B. Messnarz, F. Hanser, and B. Tilg. Lead field computation for the electrocardiographic inverse problem–finite elements versus boundary elements. *Computer Methods and Programs in Biomedicine*, 77(3):241–252, March 2005.

[9] A.A. Oberai, N.H. Gokhale, and G.R. Feijo. Solution of inverse problems in elasticity imaging using the adjoint method. *Inverse Problems*, 19(2):297, April 2003.

[10] G. Xu, H. Wu, S. Yang, S. Liu, Y. Li, Q. Yang, W. Yan, and M. Wang. 3-D electrical impedance tomography forward problem with finite element method. *IEEE Transactions on Magnetics*, 41(5):1832–1835, May 2005.

[11] M. Soleimani, W.R.B. Lionheart, A.J. Peyton, Xiandong Ma, and S.R. Higson. A three-dimensional inverse finite-element method applied to experimental eddy-current imaging data. *IEEE Transactions on Magnetics*, 42(5):1560–1567, May 2006.

[12] I.T. Rekanos, T.V. Yioultsis, and T.D. Tsiboukis. Inverse scattering using the finite-element method and a nonlinear optimization technique. *IEEE Transactions on Microwave Theory and Techniques*, 47(3):336–344, March 1999.

[13] O. Cmielewski, H. Tortel, A. Litman, and M. Saillard. A two-step procedure for characterizing obstacles under a rough surface from bistatic measurements. *IEEE Transactions on Geoscience and Remote Sensing*, 45(9):2850–2858, September 2007.

[14] R. Lencrerot, A. Litman, H. Tortel, and J.-M. Geffrin. Measurement strategies for a confined microwave circular scanner. *Inverse Problems in Science and Engineering*, 17(6):787–802, September 2009.

[15] R. Lencrerot, A. Litman, H. Tortel, and J.-M. Geffrin. Imposing Zernike representation for imaging two-dimensional targets. *Inverse Problems*, 25(3):035012, March 2009.

[16] W. Rachowicz and A. Zdunek. Application of the FEM with adaptivity for electromagnetic inverse medium scattering problems. *Computer Methods in Applied Mechanics and Engineering*, 200(2932):2337–2347, July 2011.

[17] C. Schwarzbach and E. Haber. Finite element based inversion for time-harmonic electromagnetic problems. *Geophysical Journal International*, 193(2):615–634, May 2013.

[18] A. Litman, R. Lencrerot, and J.-M. Geffrin. Combining spatial support information and shape-based method for tomographic imaging inside a microwave cylindrical scanner. *Inverse Problems in Science and Engineering*, 18(1):19–34, January 2010.

[19] A. Zakaria, I. Jeffrey, and J. LoVetri. Full-vectorial parallel finite-element contrast source inversion method. *Progress In Electromagnetics Research*, 142:463–483, 2013.

[20] A. Quarteroni. *Domain Decomposition Methods in Science and Engineering: The Sixth International Conference on Domain Decomposition, June 15-19, 1992, Como, Italy*. American Mathematical Soc., 1994.

[21] C. Farhat and J. Mandel. The two-level FETI method for static and dynamic plate problems part I: An optimal iterative solver for biharmonic systems. *Computer Methods in Applied Mechanics and Engineering*, 155(1-2):129–151, March 1998. WOS:000073170300008.

[22] C. Farhat, P. S. Chen, J. Mandel, and F. X. Roux. The two-level FETI method part II: Extension
to shell problems, parallel implementation and performance results. *Computer Methods in Applied Mechanics and Engineering*, 155(1-2):153-179, March 1998. WOS:000073170300009.

[23] C. Farhat, A. Macedo, and M. Lesoinne. A two-level domain decomposition method for the iterative solution of high frequency exterior Helmholtz problems. *Numerische Mathematik*, 85(2):283-308, April 2000.

[24] Y. Boubendir, X. Antoine, and C. Geuzaine. A quasi-optimal non-overlapping domain decomposition algorithm for the Helmholtz equation. *Journal of Computational Physics*, 231(2):262-280, January 2012.

[25] Y. Li and J.-M. Jin. A vector dual-primal finite element tearing and interconnecting method for solving 3-D large-scale electromagnetic problems. *IEEE Transactions on Antennas and Propagation*, 54(10):3000-3009, October 2006.

[26] Y.-J. Li and J.-M. Jin. A new dual-primal domain decomposition approach for finite element simulation of 3-D large-scale electromagnetic problems. *IEEE Transactions on Antennas and Propagation*, 55(10):2803-2810, October 2007.

[27] Y.-J. Li and J.-M. Jin. Implementation of the second-order ABC in the FETI-DPEM method for 3d EM problems. *IEEE Transactions on Antennas and Propagation*, 56(8):2765-2769, August 2008.

[28] V. Dolean, S. Lanteri, and R. Perrussel. A domain decomposition method for solving the three-dimensional time-harmonic Maxwell equations discretized by discontinuous Galerkin methods. *J. Comput. Phys.*, 227(3):2044-2072, January 2008.

[29] V. Dolean, M. Gander, and L. Gerardo-Giorda. Optimized Schwarz methods for Maxwell's equations. *SIAM Journal on Scientific Computing*, 31(3):2193-2213, January 2009.

[30] M.-F. Xue and J.-M. Jin. Nonconformal FETI-DP methods for large-scale electromagnetic simulation. *IEEE Transactions on Antennas and Propagation*, 60(9):4291-4305, September 2012.

[31] E.A. Attardo, G. Vecchi, and L. Crocco. Contrast source extended Born inversion in noncanonical scenarios via FEM modeling. *IEEE Transactions on Antennas and Propagation*, 62(9):4674-4685, September 2014.

[32] I. Voznyuk, H. Tortel, and A. Litman. Scattered field computation with an extended FETI-DPEM2 method. *Progress In Electromagnetics Research*, 139:247-263, 2013.

[33] I. Voznyuk, H. Tortel, and A. Litman. 3D electromagnetic scattering computation in free-space with the FETI-FDP2 method. *IEEE Transactions on Antennas and Propagation*, accepted, 2015.

[34] A.G. Tijhuis and A.P.M. Zwamborn. Marching on in anything: Solving electromagnetic field equations with a varying physical parameter. In *Ultra-Wideband, Short-Pulse Electromagnetics 5*, pages 655-662. Springer US, January 2002.

[35] A.G. Tijhuis, M.C. van Beurden, and A.P.M Zwamborn. Iterative solution of field problems with a varying physical parameter. *Elektrik, Turkish Journal of Electrical Engineering and Computer Sciences*, 10(2):163-183, 2002.

[36] A.G. Tijhuis, M.C. van Beurden, and A.P.M. Zwamborn. Iterative solution of field problems with a varying physical parameter. In *Ultra-Wideband, Short-Pulse Electromagnetics 6*, pages 223-230. Springer US, January 2003.

[37] A. Abubakar and T.M. Habashy. Application of the MR-CSI method for three-dimensional imaging of the triaxial induction measurements. *IEEE Transactions on Geoscience and Remote Sensing*, 48(6):2613-2619, June 2010.

[38] M. D'Urso, T. Isernia, and A.F. Morabito. On the solution of 2-D inverse scattering problems via source-type integral equations. *IEEE Transactions on Geoscience and Remote Sensing*, 48(3):1186-1198, March 2010.

[39] P.-A. Barriere, J. Idier, J.-J. Laurin, and Y. Goussard. Contrast source inversion method applied to relatively high contrast objects. *Inverse Problems*, 27(7):075012, July 2011.

[40] R. Byrd, P. Lu, J. Nocedal, and C. Zhu. A limited memory algorithm for bound constrained
optimization. SIAM Journal on Scientific Computing, 16(5):1190–1208, September 1995.

[41] J. Nocedal and S.J. Wright. *Numerical Optimization*. Springer Series in Operations Research. 2006.

[42] G.A. Newman and P.T. Boggs. Solution accelerators for large-scale three-dimensional electromagnetic inverse problems. *Inverse Problems*, 20(6):S151, 2004.

[43] C. Zhu, R.H. Byrd, P. Lu, and J. Nocedal. Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization. *ACM Trans. Math. Softw.*, 23(4):550–560, December 1997.

[44] J.-M. Jin. *The Finite Element Method in Electromagnetics*. Wiley-IEEE Press, New York, 2 edition, May 2002.

[45] W.C. Chew. *Waves and fields in inhomogeneous media*. Van Nostrand Reinhold, 1990.

[46] J.-M. Geffrin and P. Sabouroux. Continuing with the Fresnel database: experimental setup and improvements in 3D scattering measurements. *Inverse Problems*, 25(2):024001, February 2009.

[47] J.L. Volakis, A. Chatterjee, and L.C. Kempel. Finite element method for electromagnetics: antennas, microwave circuits, and scattering applications. IEEE Press, June 1998.

[48] P. R. Amestoy, I. S. Duff, and J. Y. L'Excellent. Multifrontal parallel distributed symmetric and unsymmetric solvers. *Computer Methods in Applied Mechanics and Engineering*, 184(24):501–520, April 2000.

[49] G. Karypis and V. Kumar. A parallel algorithm for multilevel graph partitioning and sparse matrix ordering. *International Conference on Parallel Processing*, pages 113–122, 1995.

[50] Y. Saad and M.H. Schultz. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Stat. Comput., 7(3):856–869, July 1986.

[51] C. Eyraud, J.-M. Geffrin, P. Sabouroux, P.C. Chaumet, H. Tortel, H. Giovannini, and A. Litman. Validation of a 3D bistatic microwave scattering measurement setup. *Radio Science*, 43(4):RS4018, 2008.

[52] D. Bertsekas. *Constrained Optimization and Lagrange Multiplier Methods*. Academic Press, 1982.

[53] O.M. Bucci and G. Franceschetti. On the degrees of freedom of scattered fields. *IEEE Transactions on Antennas and Propagation*, 37(7):918–926, July 1989.

[54] C. Eyraud, J.-M. Geffrin, A. Litman, and J.-P. Spinelli. A large 3D target with small inner details: A difficult cocktail for imaging purposes without a priori knowledge on the scatterers geometry. *Radio Science*, 47(5):RS0E23, 2012.

[55] S. Nounouh, C. Eyraud, A. Litman, and H. Tortel. Quantitative imaging with incident field modelling from multistatic measurements on line segments. *IEEE Antennas and Wireless Propagation Letters*, 14(1):253–256, 2015.

[56] A. Litman, D. Lesselier, and F. Santosa. Reconstruction of a two-dimensional binary obstacle by controlled evolution of a level-set. *Inverse Problems*, 14(3):685, June 1998.