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

A Quasi-Newton method for reconstructing the constitutive parameters of three-dimensional (3D) 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. 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.


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 illposedness 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 FEM. Indeed, the discretization mesh in FEM is directly proportional to the wavelength b λ . When the frequency increases, the number of field unknowns increases very rapidly, especially in three-dimensional (3D) configurations. For example, a 30 b 3 λ computational domain corresponds to approximatively 0. 5 10 6 × field unknowns and 0. 5 10 6 × permittivity unknowns, while a domain of 90 b 3 λ meshed this time for a frequency which is simply doubled leads to almost 1.5 10 6 × field unknowns and 1.4 10 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 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 developped and thoroughly tested [32,33]. Our aim here is to efficiently couple 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 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] 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 method, the unknowns correspond to the current distribution and to the permittivity. This method has been extended to various problems [3,[36][37][38]. Here, we focus on the second, socalled '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,39], 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 [40,41]. 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 [42] 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 FEM 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 [43,44] 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 [5,45]. Results obtained with either a classical FEM-based inversion scheme and the FETI-FDP2-based inversion scheme are compared and discussed. Conclusions follow.

Electromagnetic configuration
Let us consider an isotropic 3D inhomogeneous dielectric object embedded in an isotropic background with known permittivity r r , where 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 r ( ( ) ) In the following, even if the modelling capabilities enable it, we will assume that the embedding background is homogeneous, i.e., where k 0 corresponds to the vacuum wave-number. The secondary induced currents are defined as As the sources are located in far-field, the incident field is modelled as a plane wave E r k p k r ( ; ) exp( j · ) where p s corresponds to the polarization state of the incident field and k s to its wave-vector. The scattered field is collected on a probing line far Γ also located in far-field. In order to limit the computational domain Ω, we take advantage of the NF-FF transform [43,44] The closed surface Σ is chosen inside Ω such that it completely encloses the scatterers and thus the investigation test domain . n′ denotes the unit normal vector to Σ at r′ and points toward the exterior region. G r r ( ; ) 0 ′ is the classical 3D free-space scalar Green function.

Classical FEM
In order to simulate electromagnetic wave scattering in various types of complex environments, a FEM has been implemented [43,46]. The computational domain Ω is first bounded and radiation boundary conditions are applied in order to mimic infinite free-space [43]. 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 E sc are directly proportional to the number of edges N edges in Ω. These values are computed by inverting the discretized version of the weak form associated to equation (1) where  is the sum of the stiffness matrix, the mass matrix and the external boundary conditions matrices on Ω ∂ and f sc represents the discretized version of the weak form associated to the induced current terms J sc . Contrary to the method of moments where the linear system is dense, the matrix  is sparse and can be easily stored thanks to compressed storage algorithms. The linear system is inverted with a frontal method [47] specially adapted for sparse matrices. Once E sc is computed inside Ω, it can be transformed into E sc far thanks to the near-field to far-field transform, where Σ corresponds to the surface of a cubic box of slightly smaller dimensions than Ω.

FETI-FDP2 FEM
The FEM 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 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 An arbitrary mesh partitioner [48] yielding strongly irregular surfaces with varying curvatures is employed, leading to subdomains i Ω with no specific planar interfaces (figure 1). In each of these subdomains i Ω , local linear systems can now more easily be inverted. The discretized version of the weak form associated to equation (1) then reads where i  is the sum of the stiffness matrix, the mass matrix and the external and internal boundary conditions matrices on Ω ∂ and i Ω ∂ . E i sc comprises the unknown components of the vector E sc in i Ω and is directly proportional to the number of edges N i edges in i Ω . f i sc represents the discretized version of the weak form associated to the induced current terms J sc for elements in i Ω . i  is a Boolean matrix which extracts only the edges of i Ω which are on i Ω ∂ .

Inverse Problems 31 (2015) 085005 I Voznyuk et al
As for the dual Lagrange multipliers i sc Λ , 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 ∂ . By combining the previous linear systems, we end up with the so-called interface problem, d , sc sc sc contains all the Lagrange multipliers. The explicit definition of the matrix sc  and the right-hand-side term d sc can be found in [33]. For a given illumination p k ( , ) s s and a given permittivity distribution r ε , the forward scattering problem requires first to solve the interface problem provided in equation (7). Then, the solution inside each subdomain i Ω can be evaluated independently by solving equation (5). When the size of the computational domain Ω 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 Ω to limit the number of discretization cells in each subdomain i Ω (table 1). With such a strategy, equation (5) can still be inverted with a direct frontal method specially adapted for sparse matrices [47]. 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 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 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 Ω (a large collection of little sub-domains easily invertible (equation (5)) combined with a large interface problem (equation (7)) or a small N Ω (a small collection of large sub-domains combined and 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 [49], in order to limit the memory storage even if we loose 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 is generated with the GMRES algorithm, for a given initial guess ( ) s r sc, ε and a stopping criterion η. G s corresponds to the final number of GMRES iterations. Once the field is computed everywhere in Ω, the NF-FF transform provides E sc far .

Inverse scattering problem
The scattered field measurements E mes are acquired on a probing line far Γ far from the target position. The receiving antenna, located at r r , is assumed to be linearly polarized along q r , q r being a unit vector. The field actually measured for each receiving point is thus E r k q ( ; ) · r s r sc far . Two possibilities occur to compare the simulated field and the measured field. Either the far-field measurement is back-propagated from far Γ into the computational domain Ω with a far-field-to-near-field (FF-NF) transform, or the simulated field is extracted from Ω and projected onto 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 illposed 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, it would require to measure all the polarization cases 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 r ε , to where w s r , corresponds to weighting coefficients enabling to balance each measurement point with respect to its own measurement uncertainty [50]. N rec (respectively N src ) corresponds to 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)- (3). In order to take into account all constraints, the following Lagrangian functional is introduced  where U and P are Lagrange multipliers. · · far 〈 | 〉 Γ corresponds to the L 2 inner product on far Γ and · · 〈 | 〉 Ω to the L 2 inner product on Ω. As we select an iterative optimization scheme based on a gradient-descent direction algorithm, we need to compute the gradient ( ) r r   ε ε of the misfit criterion. This quantity can be directly evaluated when the Karush-Kuhn-Tucker conditions [51] are met, which enables to find the saddle-point of . At this saddle-point, equations (1) and (2) The operator †  corresponds to the adjoint operator of  , such that The notation u* corresponds to the complex conjugate of u. The adjoint field thus satisfies a similar equation as E sc (see equation (1)) apart from the righthand-side term, where receivers now act as sources and emit the discrepancy existing between the simulated and measured fields.
Unfortunately, the adjoint operator †  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 †  , we follow an other strategy. As the receivers are in farfield, we construct incident plane waves arising from the direction k r . The adjoint source field term is then computed by summing up these plane waves with appropriate weighting factors 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 Since a solution of the optimization problem has to be sought numerically, a parameter representation of the complex permittivity function r ( ) 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  of the computational domain Ω. As in any iterative scheme, is generated for a given initial guess r ε . These permittivity maps only differ from b ε within the investigation domain . The sequence of permittivity maps is constructed thanks to a Quasi-Newton algorithm α is the step length. Among all the different approximations for computing the matrix , we follow the classical BFGS method [39] 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 N tets tets   × ∈ ∈ , where N tets  ∈ is the number of tetrahedrons inside . We have thus adopted the limited-memory version of the BGFS algorithm [42].
The inversion iterations are stopped if the cost function ( ) r  ε is becoming lower than a stopping criterion ξ, or if the cost function stagnates or if a maximum number of iterations (here 100) is reached.

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 FEMs 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 FEM 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.

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 Inverse Problems 31 (2015) 085005 I Voznyuk et al criteria are thus introduced for that purpose. The first one relies on the convergence of the cost function  , 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, This criterion is either computed all along the iterations or for the last iteration of each method.

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 TWOCUBES object from the Fresnel database [45]. The choice of this target is dictated by its size and its simplicity. The TWOCUBES object consists of two identical cubes with permittivity 2.25 ). The investigation domain  is a sphere with a radius of 0.05 m centred in Ω. The finite element discretization is set to 10 points per wavelength, leading to N 15 197 tets  = ∈ in . The initial guess r (0) ε corresponds to rb ε apart from a little cube with a 0.01 m side length, with one corner at the center of the setup 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.

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 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 .
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 i  ), the electromagnetic configuration (positions and orientations of the emitters and receivers, plane waves illuminations, etc). For the permanent subdomains, the matrices i  are computed once for all as well as some parts of the 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 r v ( ) ε . The corresponding matrices i  , the non-permanent part of the sc  matrix as well as the right-handside terms for both the direct and the adjoint problems are thus updated at each inversion iterative step. Λ ε must now be properly defined. We investigate three different approaches based on the so-called marching-on-in-anything technique [34,35] which has already been shown to accelerate the convergence procedures. For • {sc, adj} ∈ , these approaches are based on (i) a blank initialization: (iii) a marching-on-in-anything source initialization: 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 η is set to 10 2 − . The results associated to the various strategies are presented in figure 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 η, 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 field, due to the choice of η. 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 (respectively H s ) associated to the direct (respectively 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 ±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 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.

Forward problem stopping criterion
The stopping criterion η 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 η to the computer precision. Unfortunately, as practice shows, the multi-sources calculation is going to be tremendously slow. Moreover, taking a very small η 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 η should mainly influence the gradient computation.
Different stopping criteria are used. The convergence behaviour of  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 η 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 have thus adopted a strategy for dynamically selecting η according to the value of ( ) gets below 10 2 − , we fix η 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.

Inversion of experimental data
The aforementioned inversion algorithm is now confronted to measurements extracted from the 3D Fresnel database [45]. 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  is a sphere centred at (0, 0, 0.025) m with a diameter of 0.16 m ( 2.1 b λ ∼ ). This test-domain has been selected as it enables to comprise all the targets from the database. The computational domain Ω is a cubic domain, centred at (0, 0, 0.0) m with a side-length of 0.26 m ( 3.5 b λ ∼ ). The initial guess r (0) ε differs from 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 been performed on an Intel(R) Xeon(R) CPU X5570 at 2.93 GHz, 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 N edges edges × matrix. Instead, the largest matrix system is associated to the interface problem and consists of a N N × Λ Λ matrix. As the time requirement highly depends on the number of GMRES iterations, we will discuss it for each object under reconstruction.

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 happen that the real part and the imaginary part of the permittivity take non-physical values.  38 12 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 s for the FETI-FDP2based scheme against 450 s for the FEM solution. As no parallelization procedure is followed,

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. 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 [52]. When the size of the individual spheres is getting too small with respect to the wavelength (for example for the CUBESPHERE, each sphere has a diameter approximatively equal to 0.2 b λ ), the inversion algorithm fails to properly separate them. One should also note that we set w 1 s r , = and thus the measurement uncertainty associated to each measurement point is not taken into account, which degrades the reconstructed maps [53]. 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 η 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  (4)). Instead, the time associated to the FETI-DP2-based inversion varies due to the GMRES iterative solver. This computational time is sligthly 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.

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 3D 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-FPD2 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 [54]. Additional a priori information will be introduced based on the target properties (homogeneity [55], limited spatial support [18], etc) in order to improve the spatial resolution of the reconstructed permittivity map. Finally, largescale configurations will also be studied. To that end, a parallelization of the entire scheme will be performed.