A non-iterative sampling approach using noise subspace projection for EIT

This study concerns the problem of the reconstruction of inclusions embedded in a conductive medium in the context of electrical impedance tomography (EIT), which is investigated within the framework of a non-iterative sampling approach. This type of identification strategy relies on the construction of a special indicator function that takes, roughly speaking, small values outside the inclusion and large values inside. Such a function is constructed in this paper from the projection of a fundamental singular solution onto the space spanned by the singular vectors associated with some of the smallest singular values of the data-to-measurement operator. The behavior of the novel indicator function is analyzed. For a subsequent implementation in a discrete setting, the quality of classical finite-dimensional approximations of the measurement operator is discussed. The robustness of this approach is also analyzed when only noisy spectral information is available. Finally, this identification method is implemented numerically and experimentally, and its efficiency is discussed on a set of, partly experimental, examples.


Introduction
Electrical impedance tomography (EIT) is an imaging technique for the reconstruction of objects embedded in a given conductive background medium . Applications range over a broad spectrum such as non-destructive material testing or tumor detection in medical imaging. This approach aims at determining the internal electrical conductivity map γ of the perturbed domain considered from boundary measurements of a current f and the associated electric potential u. These measurements represent respectively the Neumann and Dirichlet boundary data of the corresponding problem of diffusion. This inverse problem is mathematically ill-posed and its solution generally requires the knowledge of boundary data (provided by the measurements) that are 'overdetermined' relative to what is normally necessary for solving a well-posed forward (i.e. direct) problem. The Neumann-to-Dirichlet operator γ : f → u| ∂ is linear, however the operator γ → γ is nonlinear and it turns out that the inverse conductivity problem is severely ill-posed [1,2], lacking stability with respect to the input data.
In such situations, linearization techniques are often too restrictive, either in the context of physical configurations they can accommodate or the information they can provide. Moreover, the minimization-based approaches that exploit the data through a misfit cost function and have a potential of overcoming the latter restrictions unfortunately bear significant computational cost associated with repeated solutions to the forward problem. Traditional gradient-based optimization is a computationally reasonable alternative for solving the featured class of inverse problems, however, their performance depends on choosing adequately the initial guess (location, geometry, conductivity) of the hidden objects. For an overview on the subject one can refer to the review articles [3][4][5] and the references therein.
Over the past two decades, the above considerations led to the paradigm shift in mathematical theories of inverse problems initiating the development of the so-called qualitative methods for non-iterative object reconstruction from remote measurements (see e.g. [6] in the context of inverse scattering). These techniques, which provide a powerful alternative to the customary minimization approaches and linearization approximations, are commonly centered around the construction of an indicator function, that depends on an interior sampling point. Such an indicator function is normally designed to reach extreme values when the sampling point belongs to the support of the hidden flaw (or the set thereof), thereby providing a computationally-effective platform for geometric defect reconstruction. Among the diverse field of methods using approaches that can be classified as probe or sampling techniques [7] one may mention the so-called factorization method [8], the probe method and the point source method [9,10], the topological sensitivity approach [11,12] and the MUltiple SIgnal Classification (MUSIC) algorithm [13,14] among the most prominent examples.
The approach adopted in the present study is based on the non-iterative method initially proposed in [15] and [16][17][18] for the inverse conductivity problem, and which has been generalized in a variety of models [19,6]. This method can be considered as a generalization of the MUSIC algorithm in [20] and it is strongly connected with the factorization method presented in [21][22][23][24][25] for EIT. Its interest lies in the fact that it avoids the issue of the nonlinearity and it does not require any a priori information on the topology or the conductivity of the hidden object(s). Based on the solution of a linear equation that features a singular solution to the diffusion equation, this approach allows the detection of the geometrical support of the inclusions in a non-iterative framework. The aim of the present study is the proposition of a robust convergence criterion in order to characterize the range of the compact linear operator related to the above mentioned integral equation. By making use of the projection onto the associated noise subspace, this approach constitutes an alternative to the commonly used Picard series convergence criterion to characterize the operator range. This reconstruction scheme has to be related to the recent original proposition made in the article [26] for inverse scattering and which has been justified in [27] as a generalization of the MUSIC algorithm to extended objects. A comparable approach has also been successfully implemented numerically in [28] for inverse scattering problems in electromagnetism. In the MUSIC algorithm, a multi-static response matrix is computed and an indicator function is derived from the projection of a given fundamental singular test function onto its nullspace. In the present study, a similar argument is used by replacing the projection onto the range of the operator synthesizing the measurements, which is commonly associated with the Picard criterion, by a projection onto the space spanned by some of the singular vectors corresponding to small singular values. Moreover, one can mention that the approach proposed is not restricted to any particular background geometry.
The paper is organized as follows. An overview of the forward and inverse conductivity problem is given in section 2 with an emphasis on the MUSIC algorithm for infinitesimal inclusions and its extension to larger inhomogeneities by the factorization method. The convergence of the Picard series is discussed in this latter context. In section 3 the projection of an appropriate test function in the noise subspace of the measurement operator is introduced and the new criterion is constructed and discussed. In connection with practical problems which involve discrete measurements, the quality of common finite-dimensional approximations of the Neumann-to-Dirichlet operator is evaluated in section 4. Finally, section 5 presents a set of numerical results, that feature both synthetic and experimental data, to assess the efficiency and robustness of a reconstruction scheme that employs the aforementioned criterion.

Preliminaries
Let ⊂ R d , d = 2, 3, denote a bounded and connected background domain, with Lipschitz boundary ∂ and unit reference conductivity. Consider a finite union I = K j=1 j of disjoint inclusions j ⊂ such that \I is connected, with corresponding real-valued conductivities γ j ∈ L ∞ ( , R), j = 1, . . . , K. One defines the scalar conductivity map as For simplicity we further assume that 0 < c γ j < 1. This assumption could be changed into 1 < γ j C (compare remark 1).
Applying a current distribution f ∈ L 2 (∂ ) which satisfies ∂ f dS = 0, the potential u that arises in is solution of the following problem with imposed Neumann condition where n is the outward unit normal on the boundary ∂ . On introducing the Sobolev space H 1 ( ) = ϕ ∈ H 1 ( ) : ∂ ϕ dS = 0 , the boundary value problem (3) is now interpreted in a variational sense: we seek u ∈ H 1 ( ) such that Existence and uniqueness of solution follow from a Poincaré inequality and the Lax-Milgram lemma. Then, on denoting L 2 (∂ ) = ϕ ∈ L 2 (∂ ) : ∂ ϕ dS = 0 , one introduces the Neumann-to-Dirichlet (NtD) map as : (4), together with its counterpart 1 : which corresponds to the reference problem where the conductivity γ is set to the unit background value everywhere in , i.e. with γ j = 1 for j = 1, . . . , K. The NtD maps are bounded when acting on H −1/2 (∂ ) into H 1/2 (∂ ), with the subscript indicating a mean-free property over ∂ , and they are compact operators on L 2 (∂ ). Finally, denote the measurement operator, or the relative NtD map, as = − 1 . The electrical impedance tomography problem consists in determining the conductivities γ j from . It has been proved in [29][30][31][32] that it is possible to reconstruct exactly the distribution of conductivity γ within the domain considered from the knowledge of the complete, i.e. infinite dimensional, Neumann-to-Dirichlet operator. Experimentally, on using M ∈ N currents densities f m ∈ L 2 (∂ ) applied on ∂ with the corresponding u m ∈ H 1 ( ) solutions of the problem (4), and u 1m ∈ H 1 ( ) in the homogeneous case (5), the data accessible to measurement are then the traces on ∂ of these potentials. From the corresponding measured voltage densities u m | ∂ and u 1m | ∂ for m = 1, . . . , M, one can finally form the discretized relative NtD operator M which will be analyzed in detail in section 4.

Non-iterative sampling approaches
This paper investigates a non-iterative sampling approach which aims at reconstructing geometrically the unknown set I of inclusions by extracting the information synthesized in the measurement operator. To do so, the idea is to probe the range of the relative NtD operator with a fundamental solution of the diffusion equation in which exhibits a singular behavior at a chosen sampling point z ∈ as z varies over the domain of interest.

MUSIC for infinitesimal inhomogeneities.
For completeness, we first comment on the case where the inclusions are characterized by a common scaling parameter ε > 0, i.e. j ≡ ε j = z j + εˆ j with the centers z j of the inhomogeneities and their normalized shapeŝ j for j = 1, . . . , K. Let ≡ ε denote the corresponding Neumann-to-Dirichlet map. In the limit ε → 0 it has been proven in [33,18] that the relative NtD operator ≡ ε = ε − 1 converges to a finite-rank operatorˆ , ε − ε dˆ L 2 (∂ )→L 2 (∂ ) = O ε d+ 1 2 . Let us denote by N(·, z) the Green's function for the Laplace operator in with Neumann boundary conditions, i.e. N(·, z) satisfies ∂ N(·, z) dS = 0 and solves for a point z ∈ and where δ is the Dirac delta function. Then the range ofˆ is given by where the vectors e k , k = 1, . . . , d, constitute an orthonormal basis of R d . Due to the identity (8), the operatorˆ has maximal rank dK. The main result for imaging revolves around the characterization of the point-like inhomogeneities at z j by establishing a relationship between the range of the operator and the traces on the boundary ∂ of dipole functions. To state such a relation, let G z,d = d · ∇ z N(·, z) denote a test function expressed in terms of a given arbitrary unit vector d ∈ S d−1 = {x ∈ R d , |x| = 1} and the Neumann function N(·, z). From (7) one finds that G z,d is harmonic in \{z} with homogeneous Neumann boundary condition on ∂ and that ∂ G z,d dS = 0. Note that G z,d is singular at the point z. Theorem 1 constitutes the key result for the justification of the MUSIC algorithm for the reconstruction of point-like inhomogeneities: A discrete version of the operatorˆ , the so-called multi-static response matrix, is computed from a finite number of imposed currents and measured voltages. Then, a discretized version of the test function g z,d is projected onto the noise subspace which is the orthogonal complement of the space spanned by the eigenfunctions associated with the significant eigenvalues. The norm of this projection is expected to be small when z ∈ {z 1 , . . . , z K } and large everywhere else. Precisely this feature is then used as an indicator for the points z 1 , . . . , z K .

Sampling for extended inclusions.
The previous approach can be transposed to the case of extended inhomogeneities using the method introduced in [17] and implemented numerically in [16]. In this reconstruction algorithm, a suitable factorization of the operator plays a central role. Key results are stated hereafter and reference to [34] can be made for detailed proofs.
Defining the contrast q < 0 by γ = 1 + q, then for f ∈ L 2 (∂ ), one has that f = ( − 1 ) f = (u − u 1 )| ∂ with u and u 1 satisfying (9) due to equations (4) and (5). Let L 2 (I ) d denote the space of the vector-valued functions in I whose components are in L 2 (I ), and define the operators A : Notably, the adjoint operator of A is given by A * : Remark 1. The operator A is compact and T is self-adjoint and coercive when q < 0. In the case where the conductivity (1) is such that γ j > 1 then q = γ − 1 > 0 and the operator T has to be defined by T h = q(h − ∇w) with w satisfying (10). From the characterization (11) of the adjoint operator A * , the above proposition implies that the range of the operator 1/2 consists of these functions that are harmonic in \I and have homogeneous Neumann boundary conditions on ∂ . The test function G z,d is such a function if the sampling point z belongs to the support I of the inhomogeneities. This is the key theorem on which the Factorization method in the context of EIT is based, and it extends the results previously given in the context of point-like inhomogeneities. The self-adjoint and compact operator possesses an eigensystem {λ j , ψ j } with positive eigenvalues λ j sorted here in decreasing order and eigenfunctions ψ j ∈ L 2 (∂ ) for j ∈ N, A function g ∈ L 2 (∂ ) satisfies g ∈ R 1/2 if and only if the series ∞ j=1 λ j −1/2 (g, ψ j ) L 2 (∂ ) ψ j converges in L 2 (∂ ). This yields the well-known Picard criterion, that has been successfully employed in a number of studies for inclusion detection.

Noise subspace projection approach
In this section, a non-iterative sampling approach based on the projection of the test function g z,d on the noise subspace of the operator is presented. This approach finds its roots in the recent mathematical justification of the MUSIC algorithm for the reconstruction of extended objects in inverse scattering problems [26,27]. In this context, the so-called signal subspace coincides with the space spanned by the eigenfunctions associated with the largest eigenvalues of the data-to-measurement operator, classically the finite-dimensional multi-static response matrix. In particular, given δ > 0 and M δ such that any The aim of this section is to show that it is possible to construct an element which belongs to the noise-subspace of the relative NtD operator and such that its inner product with the test function g z,d is arbitrarily small when the sampling point z lies inside I and large when z ∈ \I. Let M * , M * ∈ N such that M * > M * > 0 and define the density together with the functionĥ M * ,M * z,d such that Remark 2. On introducing P M * ,M * : L 2 (∂ ) → L 2 (∂ ) the orthogonal projection onto the space spanned by the eigenfunctions ψ j for j = M * , . . . , M * , then, from the definition (12), one . Given a threshold value δ > 0 representing the level above which the eigenvalues λ j can be measured with a satisfactory precision, then if M * is sufficiently large one has and thus, by construction, the functionĥ M * ,M * z,d (13) belongs to the noise subspace of the Neumann-to-Dirichlet operator.
with x ∈ , then on denoting B(z, α) the open ball of center z and radius α > 0, the following key result holds. Proof. From the definition (14) one has for x, z ∈ and thus the Cauchy-Schwarz inequality and the definition (12) yield , then the previous inequality together with corollary 1 and when z ∈ \I, from the cororally 1, one has that h M * ,M * z,d In case (a) the convergence h M * ,∞ Hence, the two limits are uniform in x and z on compact sets due to Dini's theorem. In case (b) one can also choose the radius α uniformly in z on compact sets.
From the definition (7) of the Green's function N(·, x), the solution u 1 to the Laplace equation in satisfying (5) but with imposed Neumann boundary conditionh ∈ L 2 (∂ ), is given by Moreover, since the equation (14) can be recast in for any d ∈ S d−1 and z ∈ , one finally obtains that [Gh](x) = d · ∇u 1 (x) where u 1 is the solution (17). Therefore theorem 3 implies that, if the current densityĥ M * ,M * z,d is applied on the boundary ∂ of the domain, then the perturbation due to the inclusions is negligible asĥ M * ,M * z,d belongs to the noise subspace of the operator . Moreover, the corresponding solution potential u 1 in the reference configuration is such that |d ·∇u 1 (x)| < ε for the vector d ∈ S d−1 and at any point x ∈ I, while this current density does not vanish in \I. The physical interpretation is that this potential in is characterized by current streamlines nearly orthogonal to the chosen direction d in the inclusions j . The idea of constructing non-interacting excitations that can avoid some objects while illuminating some other parts of the domain has emerged in [26] where the notion of non-scattering waves was introduced.
From theorem 3 and the identity (16), the norm of h M * ,M * z,d given by immediately satisfies the following characterization.

Error estimates
In numerical examples and when dealing with measurements, the full Neumann-to-Dirichlet operator is never at hand. To this end we investigate in this section the approximation quality of certain finite-dimensional approximations of within the framework of the continuum model (4).
Consider a set { f m } M m=1 ⊂ L 2 (∂ ) of M ∈ N linear independent current densities f m in L 2 (∂ ) that are applied on ∂ to generate the solutions u m ∈ H 1 ( ) and u 1m ∈ H 1 ( ) of the problems (4) and (5) respectively. The trace (u m − u 1m )| ∂ of these potentials is then measured using a projection operator onto a finite-dimensional space. For example, the functions { f m } might be chosen as indicator functions of subsets of the boundary, yielding a crude model for boundary electrodes. The projection operator modeling the measurements might be chosen as a projection onto a similar space spanned by indicator functions, modeling extended electrodes, or as a Lagrange interpolation projection, modeling point electrodes.
In the following, we consider the finite-dimensional spaces Further, we introduce finite-dimensional spaces G N = span{g n } N n=1 ⊂ L 2 (∂ ) spanned by N linearly independent functions g n ∈ L 2 (∂ ) and associated bounded projections for some s > 0.
Both projections P F M and P G N are not assumed to be orthogonal projections, and G N is not required to consist of mean-free functions. Typically, P G N might be an interpolation projection (for s large enough to have the interpolation operation well-defined) or a L 2 (∂ )-orthogonal projection onto G N (for s = 0). In the following, we make use of an explicit representation of P G N via bounded linear forms g n : H s (∂ ) → R, If, for instance, P G N is a Lagrange interpolation projection, then g n (g) is the point evaluation of g in one of the interpolation nodes; if P G N is an orthogonal projection onto span{g n } N n=1 , then g n (g) = (g, g n ) L 2 (∂ ) . For the basis {g n } N n=1 of G N we choose a dual basis g * n N n=1 ⊂ L 2 (∂ ) such that g n , g * n L 2 (∂ ) = δ n,n for n, n = 1, . . . , N. We introduce the finite-dimensional current-to-voltage map by only considering current patterns in F M , and by measuring the resulting potentials u m − u 1m = f m using the projection P G N . The finite-dimensional linear operator corresponding to this model is The operator NM is bounded from L 2 (∂ ) into L 2 (∂ ) if the boundary ∂ is sufficiently regular such that the projection P G N is bounded on the image space of (see the subsequent proposition 2). The finite-dimensional operator NM is characterized by the entries of the matrix NM ∈ R N×M , Indeed, using the dual basis g * For a definition of the fractional Sobolev spaces H s (∂ ) we refer to, e.g., [36] for either smooth domains and to [37, chapter 3] for polygons in case that |s| > 1.

Proposition 2.
Assume that there is s > 0, N 0 ∈ N, and C > 0 such that If ⊂ R 2 is a convex polygon such that all interior angles are less than π/(k − 1) for 2 k ∈ N, then the latter estimate holds for s k − 1/2.

Proof. Obviously,
We start to estimate the operator norm of the second term. To this end, we recall the factorization = A * TA from proposition 1. The product TA is bounded from L 2 (∂ ) into L 2 (I ) d and the bounded linear operator A * : L 2 (I ) d → L 2 (∂ ) maps h ∈ L 2 (I ) d to the solution v| ∂ of (11), that is, If we assume that the boundary is a domain with smooth boundary, basic elliptic regularity theory (see, e.g., [36, chapter 4]) implies that v is smooth in a neighborhood of the boundary, and that the mapping h → v| ∂ is bounded from If is a convex polygon such that all interior angles are less than π/(k−1) for 2 k ∈ N, then theorem 1.7 in [35] implies that v is H k regular in a neighborhood of the boundary, and the above-described definition of the spaces H s (∂ ) implies that v| ∂ ∈ H k−1/2 (∂ ). In this case, the subsequent estimates are correct if s k − 1/2.
We exploit the latter regularity result, estimating that where we used (23). Next we consider the first term of (24) and start with If we choose t = −s, then assumption (23) implies that We use the factorization = A * TA and the boundedness of A * : L 2 (I ) d → H t (∂ ) for 0 t another time, to see by duality that A is bounded from H −t (∂ ) into L 2 (I ) d for all t > 0 in case that is a smooth domain. If is a convex polygon with interior angles less than π/(k − 1), then the last statement holds true for k + 1/2 > t 0. Hence, is bounded from H −s (∂ ) into H s (∂ ) and it remains to show that the operator norms P G N H s (∂ )→L 2 (∂ ) are uniformly bounded. Again, we exploit the second estimate from (23), to obtain that P G N H s (∂ )→L 2 (∂ ) 1 + C(s)N −s . This implies the claimed estimate for NM − .

Illustrating examples
We give three examples how the estimate of proposition 2 can be applied to obtain error estimates between finite and infinite-dimensional current-to-voltage maps.

Example 1.
A simple, yet important, example is the special case when ⊂ R 2 is the unit circle, and when a finite set F M of M trigonometric polynomials is used to discretize the boundary currents. Let us also assume that finitely many Fourier coefficients corresponding to F M model the measurements.
Then P F M is the orthogonal projection onto F M = span{sin(mφ), cos(mφ), m = 1, . . . , M}. It is well known that P F M satisfies (23) for any s 0. In consequence, the finitedimensional approximation M , defined by M = P F M P F M , converges super-exponentially to as M → ∞.
A function f ∈ H 1 (∂ ) is continuous due to Sobolev's embedding lemma. By the mean-value theorem, A duality argument finally yields that The operator M = P * F M P F M models an experimental setting where current is injected via the electrodes E M m and where the mean-value of the resulting potential is measured on E m . Using the above proposition 2 we conclude that M converges linearly to as M → ∞.

Example 3.
A final class of projectors naturally arises from the discretization of the boundary term on the right of the variational problem (4) using finite elements. Here, we assume for simplicity that ⊂ R 2 is a convex polygon.
Denote by V p h 1 a discontinuous finite element space of piecewise polynomials of degree p ∈ N on a shape-regular triangulation T 1 of with mesh size h 1 > 0 (see [39, chapter 4.1.3] for a construction and further details). The projector P p h 1 maps f ∈ H s (∂ ) into V p h 1 using local L 2 -orthogonal projections on each triangle T of the triangulation: If T ∈ T 1 and if L p T is the orthogonal projection in L 2 (T) onto the polynomials of degree p on T, then Note that the latter definition implies that P p It is well known (see [39, theorem 4.3.19]) that for s 0 the estimate holds. A duality argument and the local L 2 -orthogonality of the projection P p h 1 show that Assume that we discretize the input currents f of (4) using piecewise polynomials in V p h 1 , and that we measure point values of the resulting potentials on the nodes of a possibly different shape-regular grid T 2 with mesh size h 2 > 0. We interpret these point values using Lagrangian piecewise linear and globally continuous finite elements. The corresponding interpolation operator is denoted as Q 1 h 2 , the measurements are hence Q 1 h 2 (u| ∂ ) where u is the potential corresponding to an input current in V p h 1 . It is well-known that g−Q 1 h 2 g L 2 (∂ ) Ch s 2 g H s (∂ ) for 1/2 < s 2, see, e.g., [39, theorem 4.3.20]. Note here that the interpolation Q 1 h 2 g does not necessarily integrate to zero.
The operator Q 1 h 2 P p h 1 corresponds to a finite-dimensional current-to-voltage matrix, and the proposition 2 implies that if all interior angles of are less than π/(k − 1) for 2 k ∈ N, Hence, the approximation error even tends to zero quadratically if p and k can be chosen large enough.

Approximation of the spectrum
For implementations of noise subspace projections, it is important to know how the spectrum of the matrix NM is related to the spectrum of the operator . Recall that the eigenvalues of , sorted in decreasing order according to their multiplicity, are denoted as λ j , j ∈ N. Let us, for simplicity, assume here that NM maps L 2 (∂ ) into L 2 (∂ ) (if NM maps L 2 (∂ ) into L 2 (∂ ), then we would have to extend and NM by zero to operators on L 2 (∂ ) to be able to compare their spectra). We write σ ( NM ) = λ NM j , j ∈ N where the eigenvalues λ NM j are sorted in decreasing order according to their multiplicity. The spectrum of the finite dimensional operator NM , defined in (19), contains zero as an eigenvalue with infinite multiplicity, and a finite number of further eigenvalues.

Remark 3.
The matrix M does in general not possess an eigenvalue decomposition, especially when it is perturbed by noise. Numerically, one has to resort to the singular value decomposition and use corresponding perturbation results (see, e.g., [41] or [27, lemma 5.1.]).

Indicator function
Consider the discretized version M ∈ R M×M of the operator using an orthonormal basis of M ∈ N current densities { f m } M m=1 applied on ∂ , M = N, and F M = G N (compare (20)). Since M might no longer be normal, its singular value decomposition is introduced as where the superscript T denotes the transpose operation. In (25) The robustness of the reconstruction scheme is moreover based on the ability to represent accurately the singularity of the featured fundamental solution G z,d for each sampling point of a testing grid commonly designed to probe the entire domain . This issue can be critical for a numerical implementation of the method since singular solutions are commonly poorly represented in the standard computational platforms. However, this drawback can be circumvented by directly taking advantages of the closed-form solution of the dipole potential in R d . For a given point z ∈ and a unit vector d ∈ R d consider the dipole potential z,d (ξ) = −d · ∇ ξ N(ξ, z) which is the harmonic function in R d \{z} given by where 2 = 2π and 3 = 4π . Owing to the property of the Neumann function (7), one can conclude that the difference G z,d − z,d is harmonic in \{z} and solves a boundary value problem with the imposed Neumann condition −∇ z,d · n on ∂ which, on using the definition of the Neumann-to-Dirichlet map, finally means that the test function is given by where the constant c is to be adjusted to ensure the condition ∂ g z,d dS = 0. Replacing the test function g z,d given by (27) by a suitable discretized version g z,d ∈ R M , the Picard criterion from corollary 1 is replaced by the truncated finite sum Linear regression. Some of the previous studies that have been concerned by a numerical implementation of the method (see e.g. [7]) have been based on the expectation that both the numerator {|g T z,d u m | 2 } and the denominator {σ m } exhibit an exponential decay as sequences indexed by m. Using linear regressions of these terms one is able to obtain the decay rate of the summand {σ −1 m |g T z,d u m | 2 } as a function of the sampling point z. The inclusion set I is then determined by the points where this decay rate takes the smallest values, meaning that the featured series does not converge at these points.
Partial sum. Another approach that has been used (see e.g. [22,25]) is directly based on the equation (28), as explained in the comments below, and it consists in directly computing the truncated Picard series as a function of z. The expected behavior is that, for a sufficiently large M providing a satisfactorily approximation of the Neumann-to-Dirichlet operator , this partial sum is large outside I and small inside.
However, our own experience has led to ambivalent conclusions. On the one hand, the expected decay of the numerators of the Picard coefficients has generally not been observed. In most cases, the linear regression performed poorly even with the recourse to algorithm such as the RANdom SAmple Consensus (RANSAC [43]). On the other hand, in our simulation the number of significantly large singular values is generally low, which entails that the computation of either the decay rate from linear regressions or of the partial sum is generally not very stable and may vary significantly with the number of selected singular values. Moreover, on using the sum of the few M t th first values of the Picard coefficients, the indicator function obtained is relatively sensitive to measurement noise and one has to introduce a threshold level to determine whether a given point is inside or outside of the inclusion in order to retrieve the binary character of the initial criterion. The typical situation is illustrated by the example of a single homogeneous L-shaped inclusion of conductivity γ I = 0.01. Figure 1 shows the values of the numerators, denominators and of the coefficients themselves of the Picard series computed with d = (1; 0), when the chosen sampling point lies either outside of the inclusion ( figure 1(a)) or inside ( figure 1(b)). The reader may refer to section 5.3 for details about the numerical settings.
The approach proposed in this study has been partly motivated by these observations. Interestingly, similar conclusions were given in [28] in the context of inverse scattering in electromagnetism where the authors have exposed the need for a new criterion to characterize the range of the discrete far-field operator. The algorithm they have introduced was designated as the SVD-tail method, and makes use of the singular vectors associated with the smallest singular values rather than those associated with the signal subspace. This technique is very close to the method proposed in [26,27] where an indicator function is constructed from the projection of an appropriate singular test function onto the noise subspace of the far-field operator. The approach presented in this section is not based on an accurate reconstitution of the entire noise subspace, which would potentially require the computation of a large number of very small singular values, but rather on an approximation on a few right singular vectors v m associated with the space N * . In the ensuing implementation of the method, the levels M * and M * are chosen manually, by typically detecting the abrupt changes in the behavior of the series m → σ m .
As a consequence of the previously discussed difficulties plaguing the characterization of the signal subspace projection, one introduces the following indicator function  The proposed method aims at improving the stability of the reconstruction scheme by adopting a regularizing approach to the inverse problem considered through the projection onto a subspace of expected larger dimension compare to the dimension of the signal subspace on which previous studies have focused.

Effect of noisy data
Since the indicator function is constructed from the eigenvectors spanning the noise subspace, the questions that naturally arise concern the stability of the reconstruction and whether it is preferable to work with (29) rather than with the truncated Picard series (28). In this section, estimates are provided in order to quantify how the error on the data propagates to the reconstructions provided by the two indicator functions. Comparatively, the study [27] has addressed the stability issues arising in the context of inverse scattering problems when incorrect spectral information is used in the Picard series.
Taking into account an additive measurement noise δ > 0, the original measurement matrix M is replaced by a noisy counterpart δ M such that 1 for any interior and exterior points z i and z e . In other words, if using only the correct spectral information leads to an indicator function exhibiting a contrast between the inclusions and the background which is larger than the ratio δ 2 /σ M σ δ M , then the reconstruction from data corrupted by the additive noise δ cannot be misinterpreted.

Numerical examples
In this section, a set of numerical examples which employ synthetic data are presented to assess the efficiency and robustness of the method. The following results come within the scope of the framework described by the example 3 of section 4. Both direct problems (4) and (5) are implemented in a conventional finite elements-based computational platform (Cast3M) to simulate data in various configurations. Given the number n el of triangular elements associated with a maximal mesh size h and piecewise-linear finite elements, the discretized background domain considered is the square h = [0; 1] × [0; 1] with unit conductivity and which may contain an homogeneous inclusion (or a set thereof) of conductivity γ I < 1. The reference solution u 1h together with u h in the presence of the inclusion(s) are normalized by the constrain u h (x 0 ) = u 1h (x 0 ) = 0 at an arbitrary point x 0 ∈ ∂ h . A set of M = 144 equidistributed unit nodal current densities {f m } with disjoint supports on the domain boundary, is constructed such that it constitutes an orthonormal basis. Keeping implicit the dependance to the mesh size h, the discrete version M of the relative Neumann-to-Dirichlet operator, is computed according to (20) where the entries g n (u m − u 1m ) are consistently substituted by the orthonormal projection f T n (u m − u 1m ).
The indicator function (29) is implemented from the singular value decomposition of the matrix M and the discrete version g z,d of the test function featured in theorem 2. Even if d is arbitrary in corollary 2, we avoid fixing the dipole direction a priori, therefore we introduce a set of unit vectors d k = (cos θ k ; sin θ k ) where θ k = (k − 1)π /8 and k = 1, . . . , 8. To circumvent the limitations due to the difficulties plaguing in the computation of the Green's function on the geometry considered, the formula (27) is employed. This latter expression requires the evaluation of the scalar products of the analytical solution (26) and its normal derivative on the boundary with the basis functions {f m }, as well as the knowledge of the discrete version of the operator 1 which is given by the entries f T n u 1m . The sampling point z varies on a regular sampling grid prob h of size n prob × n prob points where n prob = 40 provides a sufficient resolution. Although it might not be a computationally optimal choice, but rather to illustrate here the performance of the proposed noise-subspace projection method in a non userchosen parameters approach, the contribution of each underlying dipole direction considered Then the following graphs (figure 2 to 5) show, for clarity, the unitary indicator function Note that the cost of computing an indicator function derived from the proposed method, i.e. typically z → I M * ,M * (z), is optimized by fixing arbitrarily a direction d and choosing manually the parameters M * and M * according to theorem 3. This latter approach might be preferred for 3D problems or when n prob is so large that care must be exerted to reduce the numerical complexity.
Single inclusion. figure 2 shows the results obtained for the identification of a single large inclusion of radius r.

Multiple inclusions.
In these examples a number of K inclusions are placed in the background domain. The inclusions have a common radius value r = 0.05 and conductivity γ I = 0.01.

Non-convex inclusion.
A single L-shaped inclusion of conductivity γ I = 0.01 is identified on figure 4.
Noisy data. The noise-free configuration presented on figure 3(c) is now considered in a situation where the data are corrupted by an additive noise. The measurement matrix is perturbed according to (30) and the standard deviation σ n of the noisy matrix is chosen in order to achieve the desired value of δ.  Discussion.
The numerical examples presented in this section show that the proposed indicator function enables a qualitative identification of the set of inclusions embedded in the reference domain considered. In the absence of measurement noise, figures 2 and 3 shows that a relatively small set of inclusions, characterized by different sizes or conductivities, can be located with a satisfactory precision by the maxima of the function (29). The graphs obtained are relatively smooth, in accordance with the discussion in remark 4. This is relatively typical of the qualitative sampling approaches which are based in practice on low-dimensional approximations of the measurement operators. On figure 3(c), the reconstruction provided is less contrasted between interior and exterior points, however four spikes clearly visible permit one to evaluate the exact number of inclusions and distinguish them geometrically. Note that the quality of this reconstruction can be improved by increasing the number M of injected currents. Figure 4 highlights that, as the conductivity model and the measurements employed are of static type, the method captures the convex envelop of the non-convex object considered; however, pronounced values are obtained inside the inclusion. Finally, the performance of the indicator function in capturing the number of inclusions and their locations in a noisy configuration is satisfactory for small values of the noise level, as shown on figure 5(a), and the reconstruction is still informative for the larger value δ = 0.05 ( figure 5(b)). Thus, the use of polluted spectral information has enabled a qualitatively satisfactory and stable identification. As expected, its quality decreases for larger values of the noise level ( figure 5(c)). For comparison with the conventional Picard criterion, figure 6 shows the results obtained using the standard truncated series (28) for the configuration of figure 5. With reference to (8) as well as, e.g., [22,25], the parameter M is chosen so that σ M +1 is the first singular value below the expected measurement error, therefore M = 8 for the examples considered. As previously, the contributions of chosen dipole directions d k , k = 1, . . . , 8, are summed and a unitary version of the indicator function is plotted. For a small amount of noise, such as δ = 0.01, the image figure 5(a), that employs the noise subspace, provides a better reconstruction than figure 6(a) which uses only the signal subspace. However, for a higher amplitude of the noise, e.g. δ = 0.1, figure 5(c) is only slightly more informative than figure 6(c). Therefore, the discrepancy between figure 5 and 6 highlights the importance of using the noise subspace to obtain satisfactory and stable reconstructions, a proposition which may appear counter-intuitive at first. These results suggest that using the noise subspace is indeed crucial to enhancing the quality of such qualitative identifications of inclusions, even when handling noisy data.
It is worth underlying as a final remark that the approach developed in this paper is based on the relative Neumann-to-Dirichlet operator, which necessitates recourse to reference measurements corresponding to the defect-free configuration. This is not an issue in the synthetic data-based examples presented here, however it might be critical to use alternatively a reference-free sampling approach, for example along the lines of the method proposed in [25].

Experimental result
This paragraph presents two preliminary experimental results provided by Choquet and Alaterre in [47]. The experimental domain was a 20 cm × 28 cm carbon-paper sheet (from www.pasco.com) with currents and potentials measured from standard laboratory sources and meters as proposed in [48]. The inclusion was a circular cut in the carbon-paper, centered at x 0 = (20, 7) with radius r = 3 cm, and the discretized version of the operator M was estimated from M = 15 current densities with disjoint supports on the boundary. The measurements were done first on the reference domain and then on the domain with the inclusion. For the computation of the indicator function (29) the numerical approach described in section 5.3 was used to estimate the discrete version of g z,d . These satisfactory initial results involving a single large inclusion prove the feasibility of the method from an experimental point of view. More measurement points and precision would contribute to increasing the accuracy of the reconstruction.

Conclusion
This study concerns the development of a qualitative approach for the identification of inclusions embedded in a conducting background domain given a set of imposed currents on the domain boundary and the measurement of the corresponding external voltages. This setting provides the access to the relative Neumann-to-Dirichlet operator which, by synthesizing the measurements, encapsulates the available information on the medium internal structure. Rather than exploiting the eigenfunctions associated with the largest eigenvalues of the data-tomeasurement operator and which commonly span its signal subspace, the approach developed in this paper is based on the extraction of information from its noise subspace. An indicator function is constructed on the alternative projection of an appropriate test function onto this latter space. This approach can be interpreted as an extension to the case of extended inclusions of the MUSIC algorithm which is dedicated to the identification of point-like objects. This paper aims at discussing the new indicator function which is based on the construction of a suitable subset of the noise subspace of the measurement operator from the behavior of the computed singular values. In particular, it has been proved that this function provides a binary criterion allowing to discriminate whether a given sampling point lies within the exterior or the interior region relatively to the sought inclusion. The intended contributions of such a criterion are twofold. Firstly, it avoids the question of determining whether the Picard series converges or not, which has been a long-standing issue for a practical use of the factorization method. Secondly, in noisy environments, the projection in the large set of those eigenfunctions below the noise level can indeed improve the quality and stability of the reconstruction, a result which may appear counter-intuitive at first. Moreover, the quality of common discrete approximations of the measurement operator has been provided in order to justify that this method is practically amenable to finite-dimensional settings, relevant to practical implementations. In this context, a stability analysis is conducted to discuss the quality of the reconstruction when it employs noisy spectral information. However, obtaining precise estimates remains an open problem and it would require the development of an appropriate perturbation theory. Finally, a set of numerical and experimental results is presented to assess the performance of the proposed method.