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

Cédric Bellis, Andrei Constantinescu, Thomas Coquet, Thomas Jaravel, Armin Lechleiter

► To cite this version:

Cédric Bellis, Andrei Constantinescu, Thomas Coquet, Thomas Jaravel, Armin Lechleiter. A non-iterative sampling approach using noise subspace projection for EIT. Inverse Problems, 2012, 10.1088/0266-5611/28/7/075015. [hal-00757299](https://hal-00757299)

HAL Id: hal-00757299

https://polytechnique.hal.science/hal-00757299

Submitted on 26 Nov 2012

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
# A non-iterative sampling approach using noise subspace projection for EIT

Cédric Bellis¹, Andrei Constantinescu², Thomas Coquet³,
Thomas Jaravel³ and Armin Lechleiter⁴

¹ Department of Applied Physics and Applied Mathematics, Columbia University,
New York, USA

² Solid Mechanics Laboratory, Ecole Polytechnique, Palaiseau, France

³ Department of Mechanics, Ecole Polytechnique, Palaiseau, France

⁴ Center for Industrial Mathematics, University of Bremen, Bremen, Germany

E-mail: ceb2186@columbia.edu,
andrei constantinescu@lms.polytechnique.fr,
lechleiter@math.uni-bremen.de

**Abstract.** 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 article 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. An introductory overlook to the forward and inverse conductivity problems is followed by the exposition of the so-called Picard criterion as a characterization of the range of the relative Neumann-to-Dirichlet operator. The construction of the novel indicator function based on the noise subspace projection is then introduced and its behavior analysed. 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 examples.

*Keywords:* Electrical impedance tomography, noise subspace, sampling method.

## 1. 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 appears to be mathematically ill-posed [1], and its resolution 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 $\Lambda_{\gamma}: f \mapsto u|_{\partial\Omega}$ is linear, however the operator that maps $\gamma \mapsto \Lambda_{\gamma}$ is non-linear and it turns out that the inverse conductivity problem is severely ill-posed. Indeed, in a general configuration, given two distributions $\gamma, \gamma' \in H^{2+s}(\Omega)$ with $s > d/2$ in dimension $d=2, 3$, the standard logarithmic stability result [2] holds and it is expressed in terms of the corresponding Dirichlet-to-Neumann operators as

$$ \|\gamma - \gamma'\|_{L^{\infty}(\Omega)} \le \beta |\log \|\Lambda_{\gamma}^{-1} - \Lambda_{\gamma'}^{-1}\|_{H^{1/2}(\partial\Omega) \to H^{-1/2}(\partial\Omega)}|^{-\alpha}, $$

where $\alpha$ and $\beta$ are positive constants and $\Lambda_{\gamma,\gamma'}^{-1}$ denotes, by abuse of notation, the Dirichlet-to-Neumann operator. However, for piecewise-constant scalar distributions of the conductivity over a bounded number of known disjoint Lipschitz domains but involving unknown real values, the previous estimate is significantly improved since the stability is now of Lipschitz-type [3]

$$ \|\gamma - \gamma'\|_{L^{\infty}(\Omega)} \le \beta \|\Lambda_{\gamma}^{-1} - \Lambda_{\gamma'}^{-1}\|_{H^{1/2}(\partial\Omega) \to H^{-1/2}(\partial\Omega)}. $$

Thus, this class of problems generally entails non-uniqueness, ill-conditioning or lack of stability towards 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 [4, 5, 6] and the references therein.

Over the past two decades, the above considerations led to the paradigm shift in mathematical theories of inverse problems that have, to a large degree, focused on the development of the so-called qualitative methods for non-iterative object reconstruction from remote measurements (see e.g. [7] 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 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 [8] one may mention the so-called factorization method [9], the probe method and the point source method [10, 11], the topological sensitivity approach
[12, 13] and the Multiple Signal Classification (MUSIC) algorithm [14, 15] among the most prominent examples.

The approach adopted in the present study is based on the non-iterative method initially proposed in [16] and [17, 18, 19] for the inverse conductivity problem, and which has been generalized in a variety of models [20, 7]. This method can be considered as a generalization of the MUSIC algorithm in [21] and it is strongly connected with the factorization method presented in [22, 23, 24, 25, 26] for EIT. Its interest lies in the fact that it avoids the issue of the non-linearity and it does not require any a priori information on the topology or the conductivity of the hidden object(s). Based on the resolution of a linear equation that features a singular solution to the diffusion equation, this approach allows to detect 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 [27] for inverse scattering and which has been justified in [28] as a generalization of the MUSIC algorithm to extended objects. A comparable approach has also been successfully implemented numerically in [29] 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 geometry like circular domains on which most of the previous studies on the subject have focused.

The article is organized as follows. An overlook to 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, the Section 5 presents a set of numerical results, that feature both synthetic and experimental data, to assess for the efficiency and robustness of a reconstruction scheme that employs the aforementioned criterion.
## 2. Inverse conductivity problem

### 2.1. Preliminaries

Let $\Omega \subset \mathbb{R}^d$, $d = 2, 3$, denote a bounded and connected background domain, with Lipschitz boundary $\partial\Omega$ and unit reference conductivity. Consider a finite union $\mathcal{I} = \bigcup_{j=1}^K \Omega_j$ of disjoint inclusions $\Omega_j \subset \Omega$ such that $\Omega \setminus \mathcal{I}$ is connected, with corresponding real-valued conductivities $\gamma_j \in L^\infty(\Omega, \mathbb{R})$, $j = 1, \dots, K$. One defines the scalar conductivity map as

$$ \gamma(\xi) = \begin{cases} 1 & \text{if } \xi \in \Omega \setminus \mathcal{I} \\ \gamma_j(\xi) & \text{if } \xi \in \Omega_j, j = 1, \dots, K. \end{cases} \qquad (1) $$

For simplicity we further assume that $0 < c \le \gamma_j < 1$. This assumption could be changed into $1 < \gamma_j \le C$ (compare Remark 1).

Applying a current distribution $f \in L^2(\partial\Omega)$ which verifies $\int_{\partial\Omega} f dS = 0$, the potential $u$ that arises in $\Omega$ is solution of the following problem with imposed Neumann condition

$$ \nabla \cdot (\gamma \nabla u) = 0 \text{ in } \Omega \qquad (2) $$

$$ (\gamma \nabla u) \cdot n = f \quad \text{on } \partial\Omega, \qquad (3) $$

where $n$ is the outward unit normal on the boundary $\partial\Omega$. On using the normalization $\int_{\partial\Omega} u dS = 0$ then the solution $u$ of the diffusion equation (3) is unique.

On introducing the Sobolev space $H_\diamond^1(\Omega) = \{\varphi \in H^1(\Omega) : \int_{\partial\Omega} \varphi dS = 0\}$, the boundary value problem (3) is now interpreted in a variational sense: We seek $u \in H_\diamond^1(\Omega)$ such that

$$ \int_{\Omega} \gamma \nabla u \cdot \nabla \varphi \, dV = \int_{\partial\Omega} f \varphi \, dS \quad \forall \varphi \in H_{\diamond}^{1}(\Omega). \qquad (4) $$

Existence and uniqueness of solution follow from a Poincaré inequality and the Lax-Milgram lemma. Then, on noting $L_\diamond^2(\partial\Omega) = \{\varphi \in L^2(\partial\Omega) : \int_{\partial\Omega} \varphi dS = 0\}$, one introduces the Neumann-to-Dirichlet (NtD) map as $\Lambda : L_\diamond^2(\partial\Omega) \to L_\diamond^2(\partial\Omega)$ such that $\Lambda f = u|_{\partial\Omega}$ where $u \in H_\diamond^1(\Omega)$ solves (4), together with its counterpart $\Lambda_1 : L_\diamond^2(\partial\Omega) \to L_\diamond^2(\partial\Omega)$ defined by $\Lambda_1 f = u_1|_{\partial\Omega}$ where $u_1 \in H_\diamond^1(\Omega)$ is solution of

$$ \int_{\Omega} \nabla u_1 \cdot \nabla \varphi \, dV = \int_{\partial\Omega} f \varphi \, dS \quad \forall \varphi \in H_{\diamond}^{1}(\Omega). \qquad (5) $$

which corresponds to the reference problem where the conductivity $\gamma$ is set to the unit background value everywhere in $\Omega$, i.e. with $\gamma_j = 1$ for $j = 1, \dots, K$. The NtD maps are bounded when acting on $H_\diamond^{-1/2}(\partial\Omega)$ into $H_\diamond^{1/2}(\partial\Omega)$, with the subscript $\diamond$ indicating a mean-free property over $\partial\Omega$, and they are compact operators on $L_\diamond^2(\partial\Omega)$. Finally, denote the *measurement* operator, or the relative NtD map, as $\Pi = \Lambda - \Lambda_1$.

The electrical impedance tomography consists in determining the conductivities $\gamma_j$ from $\Pi$. It has been proved in [30, 31, 32, 33] that it is possible to reconstruct exactly the distribution of conductivity $\gamma$ within the domain $\Omega$ considered from the knowledge of the *complete*, i.e. infinite dimensional, Neumann-to-Dirichlet operator. Experimentally,
on using $M \in \mathbb{N}$ currents densities $f_m \in L^2_{\diamond}(\partial\Omega)$ applied on $\partial\Omega$ with the corresponding $u_m \in H^1_{\diamond}(\Omega)$ solutions of the problem (4), and $u_{1m} \in H^1_{\diamond}(\Omega)$ in the homogeneous case (5), the data accessible to measurement are then the traces on $\partial\Omega$ of these potentials. From the corresponding measured voltage densities $u_m|_{\partial\Omega}$ and $u_{1m}|_{\partial\Omega}$ for $m = 1, \dots, M$, one can finally form the discretized relative NtD operator $\tilde{\Pi}_M$ which will be analyzed in detail in Section 4.

## 2.2. Non-iterative sampling approach

This paper investigates a non-iterative sampling approach which aims at reconstructing geometrically the unknown set $\mathcal{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 $\Omega$ which exhibits a singular behavior at a chosen sampling point $z \in \Omega$ as $z$ varies over the domain of interest.

The operator $\Pi$ is self-adjoint and compact [1], therefore there exists an eigensystem $\{\lambda_j, \psi_j\}$ with positive eigenvalues $\lambda_j$ sorted here in decreasing order and eigenfunctions $\psi_j \in L^2_{\diamond}(\partial\Omega)$ for $j \in \mathbb{N}$, such that for $f \in L^2_{\diamond}(\partial\Omega)$

$$ \Pi f = \sum_{j=1}^{\infty} \lambda_j(f, \psi_j)_{L^2(\partial\Omega)} \psi_j. \qquad (6) $$

### 2.2.1. Infinitesimal inhomogeneities
First, consider the case where the inclusions are characterized by the common scaling parameter $\varepsilon > 0$, i.e. $\Omega_j \equiv \Omega_j^{\varepsilon} = z_j + \varepsilon \hat{\Omega}_j$ with the centers $z_j$ of the inhomogeneities and their normalized shapes $\hat{\Omega}_j$ for $j = 1, \dots, K$. Let $\Lambda \equiv \Lambda^{\varepsilon}$ denote the corresponding Neumann-to-Dirichlet map. Then, in the limit $\varepsilon \to 0$, it has been proved in [34, 19] that the relative NtD operator $\Pi \equiv \Pi^{\varepsilon} = \Lambda^{\varepsilon} - \Lambda_1$ converges to a finite-rank operator $\hat{\Pi}$ in the operator norm of $\mathcal{L}(L^2_{\diamond}(\partial\Omega), L^2_{\diamond}(\partial\Omega))$ as

$$ \|\Pi^{\varepsilon} - \varepsilon^{d} \hat{\Pi}\|_{L^2_{\diamond}(\partial\Omega) \to L^2_{\diamond}(\partial\Omega)} = \mathcal{O}(\varepsilon^{d+\frac{1}{2}}). $$

Moreover, let $N(\cdot, z)$ denote the Green's function for the Laplace operator in $\Omega$ with respect to Neumann boundary conditions, i.e. the function $N(\cdot, z)$ such that $\int_{\partial\Omega} N(\cdot, z) dS = 0$ and which solves

$$ \nabla_{\xi} \cdot \nabla_{\xi} N(\xi, z) = -\delta(\xi - z) \quad \text{in } \Omega \qquad (7) $$

$$ \nabla_{\xi} N(\xi, z) \cdot n(\xi) = -\frac{1}{|\partial\Omega|} \quad \text{on } \partial\Omega, \qquad (8) $$

for a point $z \in \Omega$ and where $\delta$ is the Dirac delta function and the subscript $\xi$ indicates spatial derivatives w.r.t. this variable. Then the range of $\hat{\Pi}$ is given by

$$ \mathcal{R}(\hat{\Pi}) = \operatorname{span}\{e_k \cdot \nabla_z N(\cdot, z_j), k = 1, \dots, d; j = 1, \dots, K\}, \qquad (9) $$

where the vectors $e_k$, $k = 1, \dots, d$, constitute an orthonormal basis of $\mathbb{R}^d$. From the identity (9), the operator $\hat{\Pi}$ has maximal rank $dK$. On introducing an eigensystem
$\{\hat{\lambda}_j, \hat{\psi}_j\}$ of $\hat{\Pi}$ for $j = 1, \dots, dK$, and $\hat{\lambda}_j = 0$ for $j > dK$, the following convergence result holds (see [19]) in terms of eigenvalues

$$ \lambda_j \equiv \lambda_j^{\varepsilon} = \varepsilon^d \hat{\lambda}_j + \mathcal{O}(\varepsilon^{d+\frac{1}{2}}). $$

Given $m \le dK$, the orthogonal projection onto the set of the $m$ first eigenfunctions $\psi_j \equiv \psi_j^{\varepsilon}$ of $\hat{\Pi}^{\varepsilon}$ is also an approximation at the order $\varepsilon^{1/2}$ of the projector generated by $\hat{\psi}_j$ for $1 \le j \le m$.

The main result for imaging purposes revolves around the characterization of the point-like inhomogeneities at $z_j$ by establishing the relationship between the range of the operator $\hat{\Pi}$ and the traces on the boundary $\partial\Omega$ of dipoles located near the true inhomogeneities. In order to state such a relation, let $G_{z,d} = d \cdot \nabla_z N(\cdot, z)$ denote a test function expressed in terms of a given arbitrary unit vector $d \in \mathbb{R}^d$ and the Neumann function $N(\cdot, z)$. From the definition (8) of the latter, one finds that the field $G_{z,d}$ is harmonic in $\Omega \setminus \{z\}$ with homogeneous Neumann boundary condition on $\partial\Omega$ and it verifies $\int_{\partial\Omega} G_{z,d} dS = 0$. Moreover this solution has a singularity at the point $z$. Now, the characterization is given by the following proposition proved in [19].

**Theorem 1.** For any $d \in \mathbb{R}^d \setminus \{0\}$ and $z \in \Omega$, consider the trace $g_{z,d} = G_{z,d}|_{\partial\Omega}$ of the test function $G_{z,d} = d \cdot \nabla_z N(\cdot, z)$, then $z \in \{z_1, \dots, z_K\}$ if and only if $g_{z,d} \in \mathcal{R}(\hat{\Pi})$.

The Theorem 1 constitutes the key result for the justification of the MUSIC algorithm for the reconstruction of point-like inhomogeneities. In this approach, a discrete version of the operator $\hat{\Pi}$, the so-called multi-static response matrix, is computed from a finite number of imposed currents and measured voltages on the boundary of the domain. Then, from the singular value decomposition of this matrix, the vectorized 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 then expected to be small when $z \in \{z_1, \dots, z_K\}$ and large everywhere else, a feature which can then be used as an indicator function as $z$ varies over a domain of interest.

### 2.2.2. Extended inclusions

The previous approach can be transposed to the case of extended inhomogeneities using the method introduced in [18] and implemented numerically in [17]. In this reconstruction algorithm, a suitable factorization of the operator $\Pi$ plays a central role. Key results are stated hereafter and reference to [1] can be made for detailed proofs.

Defining the contrast $q < 0$ by $\gamma = 1 + q$, then for $f \in L^2(\partial\Omega)$, one has that $\Pi f = (\Lambda - \Lambda_1)f = (u - u_1)|_{\partial\Omega}$ with $u$ and $u_1$ verifying

$$ \int_{\Omega} (1+q) \nabla(u-u_1) \cdot \nabla \varphi \, dV = - \int_{\mathcal{I}} q \nabla u_1 \cdot \nabla \varphi \, dV \quad \forall \varphi \in H_{\diamond}^{1}(\Omega), (10) $$

owing to the equations (4) and (5). Moreover, let $L^2(\mathcal{I})^d$ denote the space of the vector-valued functions in $\mathcal{I}$ whose components are in $L^2(\mathcal{I})$, and define the operators
$A : L_{\diamond}^2(\partial\Omega) \to L^2(\mathcal{I})^d$ by $Af = \nabla u_{1|\mathcal{I}}$ where $u_1$ solves (5), together with $T : L^2(\mathcal{I})^d \to L^2(\mathcal{I})^d$ such that $Th = q(\nabla w - h)$ with $w \in H_{\diamond}^1(\Omega)$ solution of

$$ \int_{\Omega} (1+q) \nabla w \cdot \nabla \varphi \, dV = \int_{\mathcal{I}} qh \cdot \nabla \varphi \, dV \quad \forall \varphi \in H_{\diamond}^1(\Omega). \quad (11) $$

Notably, the adjoint operator of $A$ is given by $A^* : L^2(\mathcal{I})^d \to L_{\diamond}^2(\partial\Omega)$ with $A^*h = v|_{\partial\Omega}$ such that $v \in H_{\diamond}^1(\Omega)$ solves the equation

$$ \int_{\Omega} \nabla v \cdot \nabla \varphi \, dV = \int_{\mathcal{I}} h \cdot \nabla \varphi \, dV \quad \forall \varphi \in H_{\diamond}^1(\Omega). \quad (12) $$

**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 $\gamma_j > 1$ then $q = \gamma - 1 > 0$ and the operator $T$ has to be defined by $Th = q(h - \nabla w)$ with $w$ verifying (11) to recover the coercivity.

From (10) one obtains the following key relationship (see [1])

**Proposition 1.** Consider the operators $A$ and $T$ defined as above, then the relative Neumann-to-Dirichlet operator can be factorized as $\Pi = A^*TA$. Moreover, the ranges of $\Pi^{1/2}$ and $A^*$ coincide, i.e. $\mathcal{R}(\Pi^{1/2}) = \mathcal{R}(A^*)$.

From the characterization (12) of the adjoint operator $A^*$, the above proposition entails that the range of the operator $\Pi^{1/2}$ consists of these functions that are harmonic in $\Omega \setminus \mathcal{I}$ and have homogeneous Neumann boundary conditions on $\partial\Omega$. Then, since the test function $G_{z,d}$ is such a function if the sampling point $z$ belongs to the support $\mathcal{I}$ of the inhomogeneities, the associated current density $g_{z,d}$ on the domain boundary $\partial\Omega$ verifies the following theorem (see [18])

**Theorem 2.** Given $d \in \mathbb{R}^d \setminus \{0\}$ and $z \in \Omega$, then the trace $g_{z,d} = G_{z,d}|_{\partial\Omega}$ of the test function $G_{z,d} = d \cdot \nabla_z N(\cdot, z)$ is such that $z \in \mathcal{I}$ if and only if $g_{z,d} \in \mathcal{R}(\Pi^{1/2})$.

This is the key theorem on which the Factorization method in the context of EIT is based on, and it extends the results previously given in the context of point-like inhomogeneities.

From the spectral decomposition (6) of the compact and self-adjoint operator $\Pi$, one has for $f \in L_{\diamond}^2(\partial\Omega)$

$$ \Pi^{1/2} f = \sum_{j=1}^{\infty} \lambda_j^{1/2} (f, \psi_j)_{L^2(\partial\Omega)} \psi_j, $$

from which one can conclude that a given $g \in L_{\diamond}^2(\partial\Omega)$ verifies $g \in \mathcal{R}(\Pi^{1/2})$ if and only if the series $\sum_{j=1}^{\infty} \lambda_j^{-1/2} (g, \psi_j)_{L^2(\partial\Omega)} \psi_j$ converges in $L_{\diamond}^2(\partial\Omega)$. This property can be characterized by the well-known Picard criterion, that has been successfully employed in a number of studies to identify the geometrical support of the hidden inclusion(s). This criterion is summarized by the following corollary derived from the Theorem 2 which enables a practical reconstruction approach
**Corollary 1.** For $d \in \mathbb{R}^d \setminus \{0\}$ and $z \in \Omega$, then $z \in \mathcal{I}$ if and only if the series $\sum_{j=1}^{\infty} \frac{|(g_{z,d}, \psi_j)_{L^2(\partial\Omega)}|^2}{\lambda_j}$ converges.

The Theorem 2 and its Corollary 1 allow to reconstruct the topology and geometry of the inclusion set $\mathcal{I}$ by providing a point-by-point binary criterion if one is able to characterize the possible “blow-up” of the above series for the points $z$ that lie outside of the support of the inclusion, i.e. lying in the domain $\Omega \setminus \mathcal{I}$.

**3. 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 $\Pi$ 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 [27, 28]. 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 $\delta > 0$ and $M_\delta$ such that any $j > M_\delta$ verifies $\lambda_j \le \delta$, this space denoted as $S_\delta$ is defined by $S_\delta = \text{span}\{\psi_j, j = 1, \dots, M_\delta\}$. Then the noise subspace is the orthogonal complement $N_\delta = S_\delta^\perp$, i.e. $N_\delta = \text{span}\{\psi_j, j \ge M_\delta + 1\}$.

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 $\mathcal{I}$ and large when $z \in \Omega \setminus \mathcal{I}$. Let $M_*, M^* \in \mathbb{N}$ such that $M^* > M_* > 0$ and define the density $h_{z,d}^{M_*, M^*}$ by

$$h_{z,d}^{M_*, M^*} = \sum_{j=M_*}^{M^*} \frac{(g_{z,d}, \psi_j)_{L^2(\partial\Omega)}}{\lambda_j^{1/2}} \psi_j, \quad (13)$$

together with the function $\hat{h}_{z,d}^{M_*, M^*}$ such that

$$\hat{h}_{z,d}^{M_*, M^*} = \|h_{z,d}^{M_*, M^*}\|_{L^2(\partial\Omega)}^{-1} \sum_{j=M_*}^{M^*} \frac{(g_{z,d}, \psi_j)_{L^2(\partial\Omega)}}{\lambda_j} \psi_j. \quad (14)$$

**Remark 2.** On introducing $\mathcal{P}_{M_*, M^*}: L^2(\partial\Omega) \to L^2(\partial\Omega)$ the orthogonal projection onto the space spanned by the eigenfunctions $\psi_j$ for $j = M_*, \dots, M^*$,

$$\mathcal{P}_{M_*, M^*} f = \sum_{j=M_*}^{M^*} (f, \psi_j)_{L^2(\partial\Omega)} \psi_j, \quad f \in L^2(\partial\Omega),$$

then, from the definition (13), one has $\Pi^{1/2}h_{z,d}^{M_*, M^*} = \mathcal{P}_{M_*, M^*}g_{z,d}$.

Given a threshold value $\delta > 0$ representing the level above which the eigenvalues $\lambda_j$ can be measured with a satisfying precision, then if $M_*$ is sufficiently large one has

$$\text{span}\{\psi_j, j = M_*, \dots, M^*\} \subset N_\delta,$$
and thus, by construction, the function $\hat{h}_{z,d}^{M_*,M^*}$ (14) belongs to the noise subspace of the Neumann-to-Dirichlet operator.

Now, let $\mathcal{G}: L_{\diamond}^2(\partial\Omega) \to L^2(\Omega)$ be defined for $d \in \mathbb{R}^d \setminus \{0\}$ as

$$[\mathcal{G} h](\mathbf{x}) = (h, g_{\mathbf{x},d})_{L^2(\partial\Omega)}, \qquad (15)$$

with $\mathbf{x} \in \Omega$, then on noting $B(\mathbf{z}, \alpha)$ the open ball of center $\mathbf{z}$ and radius $\alpha > 0$, the following key result holds

**Theorem 3.** Let $\varepsilon > 0$ and $d \in \mathbb{R}^d \setminus \{0\}$.

(a) For $\mathbf{x} \in \mathcal{I}$, there exist $M_*, M^* \in \mathbb{N}$ with $M^* > M_* > 0$, such that

$$|[\mathcal{G} \hat{h}_{\mathbf{z},d}^{M_*,M^*}](\mathbf{x})| < \varepsilon \text{ if } \mathbf{z} \in \Omega.$$

(b) For $\mathbf{z} \in \Omega \setminus \mathcal{I}$, there exist $M_*, M^* \in \mathbb{N}$ with $M^* > M_* > 0$ and $\alpha > 0$, such that

$$|[\mathcal{G} \hat{h}_{\mathbf{z},d}^{M_*,M^*}](\mathbf{x})| > \frac{1}{\varepsilon} \text{ if } \mathbf{x} \in B(\mathbf{z}, \alpha).$$

If in (a) $\mathbf{x}$ belongs to a compact subset of $\mathcal{I}$, or if in (b) $\mathbf{z}$ belongs to a compact subset of $\Omega \setminus \mathcal{I}$, then one can choose $M_*, M^* \in \mathbb{N}$ with $M^* > M_* > 0$ uniformly with respect of $\mathbf{x}$ and $\mathbf{z}$, respectively.

*Proof.* From the definition (15) one has for $\mathbf{x}, \mathbf{z} \in \Omega$

$$|[\mathcal{G} \hat{h}_{\mathbf{z},d}^{M_*,M^*}](\mathbf{x})| \le \|h_{\mathbf{z},d}^{M_*,M^*}\|_{L^2(\partial\Omega)}^{-1} \sum_{j=M_*}^{M^*} \frac{1}{\lambda_j} |(g_{\mathbf{z},d}, \psi_j)_{L^2(\partial\Omega)} (\overline{g_{\mathbf{x},d}, \psi_j})_{L^2(\partial\Omega)}|, \qquad (16)$$

and thus the Cauchy-Schwarz inequality and the definition (13) yield

$$|[\mathcal{G} \hat{h}_{\mathbf{z},d}^{M_*,M^*}](\mathbf{x})| \le \left( \sum_{j=M_*}^{M^*} \frac{|(g_{\mathbf{x},d}, \psi_j)_{L^2(\partial\Omega)}|^2}{\lambda_j} \right)^{1/2} = \|h_{\mathbf{x},d}^{M_*,M^*}\|_{L^2(\partial\Omega)}.$$

(a) Let $h_{\mathbf{x},d}^{M_*,\infty} = \lim_{M_* \to \infty} h_{\mathbf{x},d}^{M_*,M^*}$, then the previous inequality together with the Corollary 1 entails that if $\mathbf{x} \in \mathcal{I}$ the sequence $M_* \mapsto \|h_{\mathbf{x},d}^{M_*,\infty}\|_{L^2(\partial\Omega)}$ is bounded and $\|h_{\mathbf{x},d}^{M_*,\infty}\|_{L^2(\partial\Omega)} \to 0$ when $M_* \to \infty$. Thus, $M_*$ sufficiently large yields the desired inequality.

(b) Alternatively, if $\mathbf{x} = \mathbf{z}$ then

$$|[\mathcal{G} \hat{h}_{\mathbf{z},d}^{M_*,M^*}](\mathbf{z})| = \|h_{\mathbf{z},d}^{M_*,M^*}\|_{L^2(\partial\Omega)}, \qquad (17)$$

and when $\mathbf{z} \in \Omega \setminus \mathcal{I}$, from the Corollary 1, one has that $\|h_{\mathbf{z},d}^{M_*,M^*}\|_{L^2(\partial\Omega)} \to \infty$ as $M^* \to \infty$. Then the second inequality can be verified for an appropriate choice of the parameter $M^*$ and owing to the continuity of $[\mathcal{G} \hat{h}_{\mathbf{z},d}^{M_*,M^*}](\mathbf{x})$ in $\mathbf{z}$ and $\mathbf{x}$.

In case (a) the convergence $\|h_{\mathbf{x},d}^{M_*,\infty}\|_{L^2(\partial\Omega)} \to 0$ when $M_* \to \infty$ is monotonic in $M_*$, and in case (b) the convergence $1/\|h_{\mathbf{z},d}^{M_*,M^*}\|_{L^2(\partial\Omega)} \to 0$ as $M^* \to \infty$ is monotonic in $M^*$. Hence, the two limits are uniform in $\mathbf{x}$ and $\mathbf{z}$ on compact sets due to Dini's theorem. In case (b) one can also choose the radius $\alpha$ uniformly in $\mathbf{z}$ on compact sets. $\square$
From the definition (8) of the Green's function $N(\cdot, \boldsymbol{x})$, the solution $u_1$ to the Laplace equation in $\Omega$ verifying (5) but with imposed Neumann boundary condition $\bar{h} \in L_{\diamond}^2(\partial\Omega)$, is given by

$$u_1(\boldsymbol{x}) = \int_{\partial\Omega} N(\boldsymbol{\xi}, \boldsymbol{x}) \bar{h}(\boldsymbol{\xi}) \, dS. \qquad (18)$$

Moreover, since the equation (15) can be recast in

$$[Gh](\boldsymbol{x}) = \int_{\partial\Omega} \boldsymbol{d} \cdot \nabla_{\boldsymbol{x}} \overline{N(\boldsymbol{\xi}, \boldsymbol{x})} h(\boldsymbol{\xi}) \, dS,$$

for any $\boldsymbol{d} \in \mathbb{R}^d \setminus \{\boldsymbol{0}\}$ and $\boldsymbol{z} \in \Omega$, one finally obtains that $[Gh](\boldsymbol{x}) = \boldsymbol{d} \cdot \nabla \bar{u}_1(\boldsymbol{x})$ where $u_1$ is the solution (18). Therefore Theorem 3 implies that, if the current density $\hat{h}_{\boldsymbol{z},\boldsymbol{d}}^{M_*,M^*}$ is applied on the boundary $\partial\Omega$ of the domain, then the perturbation due to the inclusions is negligible as $\hat{h}_{\boldsymbol{z},\boldsymbol{d}}^{M_*,M^*}$ belongs to the noise subspace of the operator $\Pi$. Moreover, the corresponding solution potential $u_1$ in the reference configuration is such that $|\boldsymbol{d} \cdot \nabla u_1(\boldsymbol{x})| < \varepsilon$ for the vector $\boldsymbol{d} \in \mathbb{R}^d \setminus \{\boldsymbol{0}\}$ and at any point $\boldsymbol{x} \in \mathcal{I}$, while this current density does not vanish in $\Omega \setminus \mathcal{I}$. The physical interpretation is that this potential in $\Omega$ is characterized by current streamlines nearly orthogonal to the chosen direction $\boldsymbol{d}$ in the inclusions $\Omega_j$. The idea of constructing non-interacting excitations that can avoid some objects while illuminating some other parts of the domain has emerged in [27] where the notion of non-scattering waves was introduced.

From the Theorem 3 and the identity (17), the norm of $h_{\boldsymbol{z},\boldsymbol{d}}^{M_*,M^*}$ given by

$$\|h_{\boldsymbol{z},\boldsymbol{d}}^{M_*,M^*}\|_{L^2(\partial\Omega)}^2 = \sum_{j=M_*}^{M^*} \frac{|(g_{\boldsymbol{z},\boldsymbol{d}}, \psi_j)_{L^2(\partial\Omega)}|^2}{\lambda_j}, \qquad (19)$$

immediately verifies the following characterization

**Corollary 2.** For $\varepsilon > 0$, $\boldsymbol{d} \in \mathbb{R}^d \setminus \{\boldsymbol{0}\}$ and compact subsets $\mathcal{K}_\mathcal{I} \subset \mathcal{I}$ and $\mathcal{K}_{\Omega \setminus \mathcal{I}} \subset \Omega \setminus \mathcal{I}$, there exist $M_*, M^* \in \mathbb{N}$ with $M^* > M_* > 0$, such that

$$\|h_{\boldsymbol{z},\boldsymbol{d}}^{M_*,M^*}\|_{L^2(\partial\Omega)}^2 < \varepsilon \text{ if } \boldsymbol{z} \in \mathcal{K}_\mathcal{I} \text{ and } \|h_{\boldsymbol{z},\boldsymbol{d}}^{M_*,M^*}\|_{L^2(\partial\Omega)}^2 > \frac{1}{\varepsilon} \text{ if } \boldsymbol{z} \in \mathcal{K}_{\Omega \setminus \mathcal{I}}.$$

# 4. Finite-dimensional approximations of NtD operators

## 4.1. Error estimates

In numerical examples and when dealing with real measurements, the full Neumann-to-Dirichlet operator $\Pi$ is never at hand. To this end we investigate in this section the approximation quality of certain finite-dimensional approximations of $\Pi$ within the framework of the continuous model (4).

Consider a set $\{f_m\}_{m=1}^M \subset L_{\diamond}^2(\partial\Omega)$ of $M \in \mathbb{N}$ linear independent current densities $f_m$ in $L_{\diamond}^2(\partial\Omega)$ that are applied on $\partial\Omega$ to generate the solutions $u_m \in H_{\diamond}^1(\Omega)$ and $u_{1m} \in H_{\diamond}^1(\Omega)$ of the problems (4) and (5) respectively. The trace $(u_m - u_{1m})|_{\partial\Omega}$ of these potentials is then measured using a projection operator onto a finite-dimensional space. Exemplarily, 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 modelling the measurements might be chosen as a projection onto a similar space spanned by indicator functions, modelling extended electrodes, or as a Lagrange interpolation projection, modeling point electrodes.

In the following, we consider the finite-dimensional spaces $F_M = \operatorname{span}\{f_m\}_{m=1}^M \subset L^2(\partial\Omega)$ together with associated projections

$$ \mathcal{P}_{FM}: L^2(\partial\Omega) \to F_M \subset L^2(\partial\Omega). $$

Further, we introduce finite-dimensional spaces $G_N = \operatorname{span}\{g_n\}_{n=1}^N \subset L^2(\partial\Omega)$ spanned by $N$ linearly independent functions $g_n \in L^2(\partial\Omega)$ and associated bounded projections

$$ \mathcal{P}_{G_N}: H^s(\partial\Omega) \to G_N \subset L^2(\partial\Omega) \quad \text{for some } s > 0. $$

Both projections $\mathcal{P}_{FM}$ and $\mathcal{P}_{G_N}$ are not assumed to be orthogonal projections, and $G_N$ is not required to consist of mean-free functions. Typically, $\mathcal{P}_{G_N}$ might an interpolation projection (for $s$ large enough to have the interpolation operation well-defined) or a $L^2(\partial\Omega)$-orthogonal projection onto $G_N$ (for $s = 0$). The $L^2(\partial\Omega)$-adjoint $\mathcal{P}_{G_N}^*$ is hence bounded from $L^2(\partial\Omega)$ into $H^s(\partial\Omega)$, and naturally, the adjoint is a projection onto the orthogonal complement of the nullspace of $\mathcal{P}_{G_N}$. In the following, we make use of an explicit representation of $\mathcal{P}_{G_N}$ via bounded linear forms $\mathfrak{g}_n: H^s(\partial\Omega) \to \mathbb{R}$,

$$ \mathcal{P}_{G_N} g = \sum_{n=1}^{N} \mathfrak{g}_n(g)g_n, \quad g \in H^s(\partial\Omega). $$

If $\mathcal{P}_{G_N}$ is a Lagrange interpolation projection, then $\mathfrak{g}_n(g)$ is the point evaluation of $g$ in one of the interpolation nodes; if $\mathcal{P}_{G_N}$ is an orthogonal projection onto $\operatorname{span}\{g_n\}_{n=1}^N$, then $\mathfrak{g}_n(g) = (g, g_n)_{L^2(\partial\Omega)}$. For the basis $\{g_n\}_{n=1}^N$ of $G_N$ we choose a dual basis $\{g_n^*\}_{n=1}^N \subset L^2(\partial\Omega)$ such that $(g_n, g_{n'}^*)_{L^2(\partial\Omega)} = \delta_{n,n'}$ for $n, n' = 1, \dots, 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} = \Pi f_m$ using the projection $\mathcal{P}_{G_N}$. The finite-dimensional linear operator corresponding to this model is

$$ \Pi_{NM} = \mathcal{P}_{G_N} \Pi \mathcal{P}_{FM}. \quad (20) $$

The operator $\Pi_{NM}$ is bounded from $L^2(\partial\Omega)$ into $L^2(\partial\Omega)$ if the boundary $\partial\Omega$ is sufficiently regular such that the projection $\mathcal{P}_{G_N}$ is well-behaved on the image space of $\Pi$ (see the subsequent Proposition 2). The finite-dimensional operator $\Pi_{NM}$ is characterized by the entries of the matrix $\tilde{\Pi}_{NM} \in \mathbb{R}^{N \times M}$,

$$ (\tilde{\Pi}_{NM})_{n,m} = \mathfrak{g}_n(u_m - u_{1m}), \quad n = 1, \dots, N, m = 1, \dots, M. \quad (21) $$

Indeed, using the dual basis $\{g_n^*\}_{n=1}^N$,

$$ (\Pi_{NM} f_m, g_n^*)_{L^2(\partial\Omega)} = (\mathcal{P}_{G_N} \Pi f_m, g_n^*)_{L^2(\partial\Omega)} \quad (22) $$

$$ = \sum_{j=1}^{N} \mathfrak{g}_j(u_m - u_{1m})(g_j, g_n^*)_{L^2(\partial\Omega)} = (\tilde{\Pi}_{NM})_{n,m} \quad (23) $$
for $1 \le n \le N$ and $1 \le m \le M$. Moreover, $(\Pi_{NM}f, g)_{L^2(\partial\Omega)} = 0$ if either $f$ or $g$ belong to the orthogonal complement of $F_M$ or $G_N$, respectively.

The next proposition investigates the approximation quality of $\Pi_{NM}$, relying on the fractional Sobolev spaces $H_{\diamond}^s(\partial\Omega)$. We concentrate on domains that are either smooth or convex polygons in $\mathbb{R}^2$, since our later computational examples use polygonal domains. It is not too difficult to state corresponding results for domains with intermediate smoothness, for non-convex polygons, and also for polyhedra in $\mathbb{R}^3$, using essentially the same technique and [35, Section 1]). However, technicalities would increase and we prefer to keep the presentation simple.

For a definition of the fractional Sobolev spaces $H_{\diamond}^s(\partial\Omega)$ we refer to, e.g., [36] for either smooth domains and to [37, Chapter 3] for polygons in case that $|s| > 1$. In the latter case, some references (e.g. [38]) denote these spaces as $H_{\text{pw}}^s(\partial\Omega)$.

**Proposition 2.** Assume that there is $s > 0$, $N_0 \in \mathbb{N}$, and $C > 0$ such that

$$ \|I - \mathcal{P}_{F_M}\|_{L_{\diamond}^2(\partial\Omega) \to H_{\diamond}^{-s}(\partial\Omega)} \le CM^{-s}, \quad \|I - \mathcal{P}_{G_N}\|_{H_{\diamond}^s(\partial\Omega) \to L^2(\partial\Omega)} \le CN^{-s} \quad (24) $$
for $N, M \ge N_0$. If $\Omega \subset \mathbb{R}^2$ or $\mathbb{R}^3$ is a smooth domain, then

$$ \|\Pi_{NM} - \Pi\|_{L_{\diamond}^2(\partial\Omega) \to L^2(\partial\Omega)} \le C(s)(N^{-s} + M^{-s}) \quad \text{for } N, M \ge N_0. $$

If $\Omega \subset \mathbb{R}^2$ is a convex polygon such that all interior angles are less than $\pi/(k-1)$ for $2 \le k \in \mathbb{N}$, then the latter estimate holds for $s \le k - 1/2$.

*Proof.* Obviously,

$$ \Pi_{NM} - \Pi = \mathcal{P}_{G_N} \Pi(\mathcal{P}_{F_M} - I) + (P_{G_N} - I)\Pi. \quad (25) $$

We start to estimate the operator norm of the second term. To this end, we recall the factorization $\Pi = A^*TA$ from Proposition 1. The product $TA$ is bounded from $L_{\diamond}^2(\partial\Omega)$ into $L^2(\mathcal{I})^d$ and the bounded linear operator $A^* : L^2(\mathcal{I})^d \to L_{\diamond}^2(\partial\Omega)$ maps $h \in L^2(\mathcal{I})^d$ to the solution $v|_{\partial\Omega}$ of (12), that is,

$$ \int_{\Omega} \nabla v \cdot \nabla \varphi \, dV = \int_{\mathcal{I}} h \cdot \nabla \varphi \, dV \quad \forall \varphi \in H_{\diamond}^{1}(\Omega). $$

If we assume that the boundary $\Omega$ 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 \mapsto v|_{\partial\Omega}$ is bounded from $L^2(\mathcal{I})^d$ into any Sobolev space $H_{\diamond}^t(\partial\Omega)$, $t \in \mathbb{R}$.

If $\Omega$ is a convex polygon such that all interior angles are less than $\pi/(k-1)$ for $2 \le k \in \mathbb{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(\partial\Omega)$ implies that $v|_{\partial\Omega} \in H^{k-1/2}(\partial\Omega)$. In this case, the subsequent estimates are correct if $s < k - 1/2$.

We exploit the latter regularity result, estimating that

$$ \begin{aligned} \|(P_{G_N} - I)\Pi\|_{L_{\diamond}^2(\partial\Omega) \to L^2(\partial\Omega)} &\le C(s)\|P_{G_N} - I\|_{H_{\diamond}^s(\partial\Omega) \to L^2(\partial\Omega)} \|\Pi\|_{L_{\diamond}^2(\partial\Omega) \to H_{\diamond}^s(\partial\Omega)} \\ &\le C(s)N^{-s}\|\Pi\|_{L_{\diamond}^2(\partial\Omega) \to H_{\diamond}^s(\partial\Omega)}, \end{aligned} $$
where we used (24). Next we consider the first term of (25) and start with

$$ \begin{aligned} \|\mathcal{P}_{G_N} \Pi(\mathcal{P}_{FM} - I)\|_{L^2(\partial\Omega) \to L^2(\partial\Omega)} &\le \|\mathcal{P}_{G_N}\|_{H^s(\partial\Omega) \to L^2(\partial\Omega)} \\ \|\Pi\|_{H^t(\partial\Omega) \to H^s(\partial\Omega)} \|\mathcal{P}_{FM} - I\|_{L^2(\partial\Omega) \to H^t(\partial\Omega)}, &\quad t \le 0. \end{aligned} $$

If we choose $t = -s$, then assumption (24) implies that

$$ \|\mathcal{P}_{FM} - I\|_{L^2(\partial\Omega) \to H^t(\partial\Omega)} \le CM^{-s} \quad M \ge N_0. $$

We use the factorization $\Pi = A^*TA$ and the boundedness of $A^* : L^2(\mathcal{I})^d \to H^t(\partial\Omega)$ for $0 \le t$ another time, to see by duality that $A$ is bounded from $H^{-t}(\partial\Omega)$ into $L^2(\mathcal{I})^d$ for all $t > 0$ in case that $\Omega$ is a smooth domain. If $\Omega$ is a convex polygon with interior angles less than $\pi/(k-1)$, then the last statement holds true for $k+1/2 > t \ge 0$.

Hence, $\Pi$ is bounded from $H^{-s}(\partial\Omega)$ into $H^s(\partial\Omega)$ and it remains to show that the operator norms $\|\mathcal{P}_{G_N}\|_{H^s(\partial\Omega) \to L^2(\partial\Omega)}$ are uniformly bounded. Again, we exploit the second estimate from (24), to obtain that $\|\mathcal{P}_{G_N}\|_{H^s(\partial\Omega) \to L^2(\partial\Omega)} \le 1 + C(s)N^{-s}$. This implies the claimed estimate for $\Pi_{NM} - \Pi$. $\square$

## 4.2. 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 in different settings.

**Example 1.** A simple, yet important, example is the special case when $\Omega \subset \mathbb{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.

In this setting, $F_M = \text{span}\{\sin(m\phi), \cos(m\phi), m = 1, \dots, M\}$ and the projection operator $\mathcal{P}_{F_N}$ is the orthogonal projection onto $F_M$. It is well known (see [39]) that $\mathcal{P}_{F_M}$ satisfies (24) for any $s \ge 0$. In consequence, the finite-dimensional approximation $\Pi_{MM}$, defined by $\Pi_{MM} = \mathcal{P}_{FM} \Pi_{FM}$, converges super-exponentially to $\Pi$ as $M \to \infty$.

**Example 2.** An important example of a non-orthogonal projection is a mean-value projection (first introduced in [40] in the context of the complete electrode model). Consider the case where one applies a constant current on certain parts of the boundary (crudely modelling electrodes), and where one only measures the mean-value of the potential on those boundary parts. We next construct the corresponding projections explicitly and prove the required error estimates.

For this example, $\Omega \subset \mathbb{R}^2$ can be a smooth domain or a convex polygon. Assume that the boundary $\partial\Omega$ is covered by $M$ disjoint, connected, and relatively open electrodes $E_m \subset \partial\Omega, m = 1, \dots, M$. Associate to each $E_m$ its indicator function $e_m \in L^2(\partial\Omega)$, that is, $e_m(x) = 0$ for $x \in E_m$ and $e_m(x) = 0$ else. Moreover, we denote by $\hat{E}_m, m = 1, \dots, M$ a connected and open superset of $E_m$ such that the closure of $\cup_{m=1}^M \hat{E}_m$
equals $\partial\Omega$ (so-called extended electrodes). The indicator function of $\hat{E}_m$ is denoted as $\hat{e}_m$. The space $F_M$ is then defined as $F_M = \text{span}\{e_m\}_{m=1}^M \cap L_{\diamond}^2(\partial\Omega)$. The projection

$$ \mathcal{P}_{F_M} g = \sum_{m=1}^{M} \frac{1}{|E_m|} \int_{\hat{E}_m} g \, dS \, e_m $$

is a non-orthogonal mean-value projection from $L_{\diamond}^2(\partial\Omega)$ into $F_M$. The $L^2$-adjoint operator $\mathcal{P}_{F_M}^*$ is given by

$$ \mathcal{P}_{F_M}^* f = \sum_{m=1}^{M} \frac{1}{|E_m|} \int_{E_m} f \, dS \, \hat{e}_m, \quad f \in L_{\diamond}^2(\partial\Omega), $$

because

$$ \int_{\partial\Omega} f \mathcal{P}_{F_M} g \, dS = \sum_{m=1}^{M} \frac{1}{|E_m|} \int_{E_m} f \, dS \int_{\hat{E}_m} g \, dS \\ = \int_{\partial\Omega} \sum_{m=1}^{M} \frac{1}{|E_m|} \int_{E_m} f \, dS \hat{e}_m g \, dS = \int_{\partial\Omega} \mathcal{P}_{F_M}^* f g \, dS. $$

A function $f \in H_{\diamond}^1(\partial\Omega)$ is continuous due to Sobolev's embedding lemma. By the mean-value theorem, $|E_m|^{-1} \int_{E_m} f \, dS = f(\mathbf{x}_0)$ for some $\mathbf{x}_0 \in E_m$, and

$$ \|f - \mathcal{P}_{F_M}^* f\|_{L_{\diamond}^2(\partial\Omega)}^2 = \sum_{m=1}^{M} \int_{\hat{E}_m} |f(\mathbf{x}) - f(\mathbf{x}_0)|^2 \, dS \le \max_{m=1,...,M} |\hat{E}_m|^2 \|f\|_{H_{\diamond}^1(\partial\Omega)}^2. $$

A duality argument finally yields that

$$ \|f - \mathcal{P}_{F_M} f\|_{L_{\diamond}^2(\partial\Omega) \to H_{\diamond}^{-1}(\partial\Omega)} \le \left( \max_{m=1,...,M} |\hat{E}_m|^2 \right)^{1/2} \|f\|_{L_{\diamond}^2(\partial\Omega)}. $$

The interesting point about this estimate is that it does not depend on the electrodes $\{E_m\}$ (especially not on their size), but only on the size of the extended electrodes $\hat{E}_m$.

Consider now a sequence of electrode configurations, leading to a sequence of spaces $\{F_M\}_{M \in \mathbb{N}}$, such that the associated extended electrodes $\{\hat{E}_m^M\}_{m=1}^M$ are quasi-uniform: there is a constant $C > 0$ independent of $M \ge M_0$ such that

$$ \max_{m=1,...,M} |\hat{E}_m^M| \le C |\hat{E}_m^M| \quad \text{for all } M \ge M_0, m = 1,...,M. $$

Summing up the last equation from $m = 1, ..., M$ and exploiting the assumption that the $\hat{E}_m$ cover $\partial\Omega$, we find that $\max_{m=1,...,M} |\hat{E}_m| \le C/M$ for all $M \ge M_0$. In consequence,

$$ \|f - \mathcal{P}_{F_M} f\|_{L_{\diamond}^2(\partial\Omega) \to H_{\diamond}^1(\partial\Omega)} \le CM^{-1}, \quad N \ge N_0. $$

The operator $\Pi_M = \mathcal{P}_{F_M}^* \Pi \mathcal{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 $\Pi_M$ converges linearly to $\Pi$ as $M \to \infty$.
**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, one of the standard techniques for numerical simulation in impedance tomography. Here, we assume for simplicity that $\Omega \subset \mathbb{R}^2$ is a convex polygon. One can extend the example to smooth domains and to dimension three using curved surface elements and 3D regularity results for the Neumann problem, respectively, but again we do not want to go into too many technicalities.

Denote by $V_{h_1}^p$ a discontinuous finite element space of piecewise polynomials of degree $p \in \mathbb{N}$ on a shape-regular triangulation $\mathcal{T}_1$ of $\Omega$ with mesh size $h_1 > 0$ (see [38, Chapter 4.1.3] for a construction and further details). The projector $P_{h_1}^p$ maps $f \in H_{\diamond}^s(\partial\Omega)$ into $V_{h_1}^p$ using local $L^2$-orthogonal projections on each triangle $\mathcal{T}$ of the triangulation: If $\mathcal{T} \in \mathcal{T}_1$ and if $L_{\mathcal{T}}^p$ is the orthogonal projection in $L^2(\mathcal{T})$ onto the polynomials of degree $p$ on $\mathcal{T}$, then

$$ (P_{h_1}^p f)|_{\mathcal{T}} = L_{\mathcal{T}}^p(f|_{\mathcal{T}}), \quad \mathcal{T} \in \mathcal{T}_1. $$

Note that the latter definition implies that $P_{h_1}^p$ maps $L_{\diamond}^2(\partial\Omega)$ into $L_{\diamond}^2(\partial\Omega)$, since

$$ \int_{\partial\Omega} P_{h_1}^p f \, dS = \int_{\partial\Omega} (f - P_{h_1}^p f) \, dS = \sum_{\mathcal{T} \in \mathcal{T}_1} \int_{\mathcal{T}} (f - L_{\mathcal{T}}^p(f|_{\mathcal{T}})) \, dS = 0. $$

It is well known (see [38, Theorem 4.3.19]) that for $s \ge 0$ the estimate

$$ \|f - P_{h_1}^p f\|_{L_{\diamond}^2(\partial\Omega)} \le C h_1^{\min(p+1,s)} \|f\|_{H_{\diamond}^s(\partial\Omega)} $$

holds. A duality argument and the local $L^2$-orthogonality of the projection $P_{h_1}^p$ show that

$$ \begin{aligned} \|f - P_{h_1}^p f\|_{H_{\diamond}^{-s}(\partial\Omega)} &\le \sup_{\phi \in H_{\diamond}^s(\partial\Omega)} \frac{|\langle f - P_{h_1}^p f, \phi \rangle|}{\|\phi\|_{H_{\diamond}^s(\partial\Omega)}} = \sup_{\phi \in H_{\diamond}^s(\partial\Omega)} \frac{|\langle f - P_{h_1}^p f, \phi - P_{h_1}^p \phi \rangle|}{\|\phi\|_{H_{\diamond}^s(\partial\Omega)}} \\ &\le \|f - P_{h_1}^p f\|_{L_{\diamond}^2(\partial\Omega)} \sup_{\phi \in H_{\diamond}^s(\partial\Omega)} \frac{\|\phi - P_{h_1}^p \phi\|_{L_{\diamond}^2(\partial\Omega)}}{\|\phi\|_{H_{\diamond}^s(\partial\Omega)}} \le C h_1^{\min(p+1,s)} \|f\|_{L_{\diamond}^2(\partial\Omega)}. \end{aligned} $$

Assume that we discretize the input currents $f$ of (4) using piecewise polynomials in $V_{h_1}^p$, and that we measure point values of the resulting potentials on the nodes of a possibly different shape-regular grid $\mathcal{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_{h_2}^1$, the measurements are hence $Q_{h_2}^1(u|_{\partial\Omega})$ where $u$ is the potential corresponding to an input current in $V_{h_1}^p$. It is well-known that $\|g - Q_{h_2}^1 g\|_{L^2(\partial\Omega)} \le C h_2^s \|g\|_{H^s(\partial\Omega)}$ for $1/2 < s \le 2$, see, e.g., [38, Theorem 4.3.20]. Note here that the interpolation $Q_{h_2}^1 g$ does not necessarily integrate to zero.

The operator $Q_{h_2}^1 \Pi P_{h_1}^p$ corresponds to a finite-dimensional current-to-voltage matrix, and the Proposition 2 implies that if all interior angles of $\Omega$ are less than $\pi/(k-1)$ for $2 \le k \in \mathbb{N}$, then

$$ \|\Pi - Q_{h_2}^1 \Pi P_{h_1}^p\|_{L_{\diamond}^2(\partial\Omega) \to L^2(\partial\Omega)} \le C(h_1^{p+1} + h_2^{\min(k-1/2,2)}) \quad \text{as } h_{1,2} \to 0. $$

In consequence, if one uses piecewise constant functions as input currents and if one measures point values of the resulting potential, then the approximation error coming
from the current discretization is dominant. If one uses piecewise linear discontinuous elements to discretize the currents, and if $h_1 = h_2$, then the approximation error tends to zero quadratically if all interior angles are less than $\pi/2$.

### 4.3. Approximation of the spectrum

For implementations of noise subspace projections, it is important to know how the spectrum of the matrix $\tilde{\Pi}_{NM}$ is related to the spectrum of the operator $\Pi$. Basic perturbation theory [41] implies that if two operators are close, then their spectra are also close. Recall that the eigenvalues of $\Pi$, sorted in decreasing order according to their multiplicity, are denoted as $\lambda_j, j \in \mathbb{N}$. We write $\sigma(\Pi_{NM}) = \{\lambda_j^{NM}, j \in \mathbb{N}\}$ where the eigenvalues $\lambda_j^{NM}$ are sorted in decreasing order according to their multiplicity. The spectrum of the finite dimensional operator $\Pi_{NM}$, defined in (20), contains zero as an eigenvalue with infinite multiplicity, and a finite number of further eigenvalues.

**Proposition 3.** Under the assumptions of Proposition 2, there is $N_0 \in \mathbb{N}$ such that

$$ \mathrm{dist}(\sigma(\Pi), \sigma(\Pi_{NM})) \leq C(s)(N^{-s} + M^{-s}) \quad N, M \geq N_0, $$

where the metric on the left is the Hausdorff distance. Further, for each $J \in \mathbb{N}$,

$$ |\lambda_j - \lambda_j^{NM}| \leq C(J, s)(N^{-s} + M^{-s}) \quad j \leq J \text{ and } N, M \geq N_0. $$

*Proof.* This is due to the estimate $\mathrm{dist}(\sigma(\Pi), \sigma(\Pi_{NM})) \leq \|\Pi - \Pi_{NM}\|_{L^2(\partial\Omega) \to L^2(\partial\Omega)}$ that holds for the spectra of bounded linear operators, see [41], and the ordering of the eigenvalues. $\square$

To investigate the relation between the spectrum of the operator $\Pi_{NM}$ and the matrix $\tilde{\Pi}_{NM}$, we assume in the following that $N = M$, such that $\tilde{\Pi}_{MM}$ is a square matrix, and write $\tilde{\Pi}_M \equiv \tilde{\Pi}_{MM}$, as well as $\Pi_M \equiv \Pi_{MM}$.

**Proposition 4.** Assume that the basis function of $F_M$ and $G_M$ equal each other for some fixed $M \geq N_0$, that is, $f_m = g_m$ for $m = 1, \dots, M$.

(a) If $(\lambda, \phi) \in \mathbb{R} \times L^2(\partial\Omega)$ is an eigenpair of the operator $\Pi_M$, then $(\lambda, \tilde{\phi})$, where $\tilde{\phi} \in \mathbb{R}^M$ has entries $\tilde{\phi}_m = g_m(\phi)$, is an eigenpair of the matrix $\tilde{\Pi}_M$.

(b) Then, if $(\lambda, \tilde{\phi}) \in \mathbb{R} \times \mathbb{R}^M$ is an eigenpair of $\tilde{\Pi}_M$, then $(\lambda, \phi)$ with $\phi \in L^2(\partial\Omega)$ such that $\phi = \sum_{m=1}^M \tilde{\phi}_m f_m$ is an eigenpair of $\Pi_M$.

*Proof.* (a) Since $\mathcal{P}_{G_M} \Pi \mathcal{P}_{F_M} \phi = \lambda \phi$, the eigenfunction $\phi = \sum_{m=1}^M \tilde{\phi}_m g_m = \sum_{m=1}^M \tilde{\phi}_m f_m$ belongs to $G_M$ (and thus, by assumption, to $F_M$). Due to Equation (23), we conclude that, for $n = 1, \dots, M$

$$ \sum_{m=1}^{M} \mathfrak{g}_n(u_m) \tilde{\phi}_m = (\Pi_N(\sum_{m=1}^{M} f_m \tilde{\phi}_m), g_n^*)_{L^2(\partial\Omega)} = \lambda(\phi, g_n^*)_{L^2(\partial\Omega)} = \lambda \tilde{\phi}_n. $$
(b) Since $(\lambda, \tilde{\phi})$ is an eigenpair of $\tilde{\Pi}_M$, it holds that $\sum_{m=1}^M \mathfrak{g}_n(u_m)\tilde{\phi}_m = \lambda\tilde{\phi}_n$ for $n = 1, \dots, M$. From (23) we conclude that

$$ \mathcal{P}_{G_M} \Pi \left( \sum_{m=1}^{M} \tilde{\phi}_m f_m \right) = \lambda \sum_{m=1}^{M} \tilde{\phi}_m g_m. $$

As above, our assumption that $f_m = g_m$ yields the claim. $\square$

If we use a sufficiently large set of basisfunctions $\{f_m\}_{m=1}^M$, then the error $\|\Pi - \Pi_M\|_{L_2^2(\partial\Omega)\to L_2^2(\partial\Omega)}$ is small and the eigenvalues and eigenvectors of the operator $\Pi$ are satisfactorily approximated by those of the matrix $\tilde{\Pi}_M$.

**Remark 3.** The matrix $\tilde{\Pi}_M$ does in general not possess an eigenvalue decomposition, especially when it is perturbed by noise. Numerically, one has hence to resort to the singular value decomposition instead of using an eigenvalue decomposition, and use corresponding perturbation results (see, e.g., [42] or [28, Lemma 5.1.]).

# 5. Numerical implementation

## 5.1. Indicator function

Consider the discretized version $\tilde{\Pi}_M \in \mathbb{R}^{M \times M}$ of the operator $\Pi$ using an orthonormal basis of $M \in \mathbb{N}$ current densities $\{f_m\}_{m=1}^M$ applied on $\partial\Omega$, $M = N$, and $F_M = G_N$ (compare (21)). Since $\tilde{\Pi}_M$ might no longer be normal, its singular value decomposition is introduced as

$$ \tilde{\Pi}_M = \mathcal{U} \Sigma \mathcal{V}^T, \qquad (26) $$

where the superscript $\tau$ denotes the transpose operation. In (26), $\Sigma$ is a $M \times M$ diagonal matrix composed of the real-valued and non-negative singular values $\{\sigma_m\}$, in decreasing order $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_M$, and the columns of the orthogonal matrices $\mathcal{U}$ and $\mathcal{V}$ are respectively the left and right singular vectors $u_m$ and $v_m$ in $\mathbb{R}^M$, $m = 1, \dots, M$, satisfying

$$ \tilde{\Pi}_M v_m = \sigma_m u_m, \quad \tilde{\Pi}_M^\tau u_m = \sigma_m v_m. $$

Given the expected measurement error $\delta > 0$ and on using the notations introduced in Section 3, the response matrix $\tilde{\Pi}_M$ possesses significant singular values $\sigma_m > \delta$ for $m \le M_\delta$ while the other eigenvalues for $m = M_\delta + 1, \dots, M$ verify $\sigma_m \le \delta$. Furthermore, the $M_\delta$ left singular vectors $u_m$ associated with significant eigenvalues form an orthonormal basis of the subspace $\mathcal{S}_\delta \subset \mathbb{R}^N$, the so-called signal subspace

$$ \mathcal{S}_\delta = \operatorname{span}\{u_m, m = 1, \dots, M_\delta\}, $$

while the right singular vectors $v_m$ for $m = M_\delta + 1, \dots, M$ form an orthonormal basis of the orthogonal complement $\mathcal{N}_\delta = \mathcal{S}_\delta^\perp$ to this space. In the context of this study and for a later numerical implementation of the method, one introduces the parameter $M_* > M_\delta$ as well as $M^*$ such that $M_* < M^* \le M$ in order to prevent taking into
account the eigenvalues of high rank which are likely to be polluted by errors in the simulated data. A companion noise subspace $\mathcal{N}_*$ can then be defined as

$$ \mathcal{N}_* = \operatorname{span}\{v_m, m = M_*, \dots, M^*\}. $$

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 $\Omega$. 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 $\mathbb{R}^d$. For a given point $z \in \Omega$ and a unit vector $d \in \mathbb{R}^d$ consider the dipole potential $\Phi_{z,d}(\xi) = -d \cdot \nabla_{\xi} N(\xi, z)$ which is the harmonic function in $\mathbb{R}^d \setminus \{z\}$ given by

$$ \Phi_{z,d}(\xi) = \frac{1}{\omega_d} \frac{(z - \xi) \cdot d}{|z - \xi|^d}, \quad \text{for } \xi \neq z, \qquad (27) $$

where $\omega_2 = 2\pi$ and $\omega_3 = 4\pi$. Owing to the property of the Neumann function (8), one can conclude that the difference $G_{z,d} - \Phi_{z,d}$ is harmonic in $\Omega \setminus \{z\}$ and solves a boundary value problem with the imposed Neumann condition $-\nabla \Phi_{z,d} \cdot n$ on $\partial\Omega$ which, on using the definition of the Neumann-to-Dirichlet map, finally entails that the test function is given by

$$ g_{z,d} = \Phi_{z,d}|_{\partial\Omega} - \Lambda_1(\nabla\Phi_{z,d} \cdot n) + c, \qquad (28) $$

where the constant $c$ is to be adjusted to ensure the condition $\int_{\partial\Omega} g_{z,d} dS = 0$.

Replacing the test function $g_{z,d}$ given by (28) by a suitable discretized version $g_{z,d} \in \mathbb{R}^M$, the Picard criterion from Corollary 1 is replaced by the truncated finite sum

$$ \sum_{m=1}^{M'} \frac{|g_{z,d}^T u_m|^2}{|\sigma_m|}, \quad 1 \le M' \le M. \qquad (29) $$

The theoretical results in [43] on a regularization technique for such discretizations of Picard series suggest a rule how to choose $M' = M'(M) \le M$ such that the latter sum converges (as the discretization parameter $M$ tends to $\infty$) if and only if $z \in \mathcal{I}$. Still, the problem is to decide for a fixed but sufficiently large discretization parameter $M$, whether sampling points $z$ are inside or outside the inclusion $\mathcal{I}$. To do so, different approaches have been used in the literature and among the most widely used ones one can distinguish the following techniques.

**Linear regression.** Some of the previous studies that have been concerned by a numerical implementation of the method (see e.g. [8]) have been based on the expectation that both the numerator $\{|g_{z,d}^T u_m|^2\}$ and the denominator $\{\sigma_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 $\{\sigma_m^{-1} |g_{z,d}^T u_m|^2\}$ as a function of the sampling point $z$. The inclusion set $\mathcal{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. [23, 26]) is directly based on the equation (29), 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 $\Pi$, this partial sum is large outside $\mathcal{I}$ and small inside.

However, our own experience have 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 [44]). 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 picked 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 $\gamma_I = 0.01$. The 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 (Fig. 1(a)) or inside (Fig. 1(b)). The reader may refer to Section 5.3 for details about the numerical settings.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Picard coefficients, numerators and denominators in log scale as functions of *m*</figcaption></figure>

The approach proposed in this study has been partly motivated by these
observations. Interestingly, similar conclusions were given in [29] 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 [27, 28] 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 \mapsto \sigma_m$.

As a consequence of the previously discussed difficulties plaguing in the characterization of the signal subspace projection, one introduce the following indicator function

**Definition 1.** For a given $d \in \mathbb{R}^d \setminus \{0\}$ and $M_*, M^* \in \mathbb{N}$ with $0 < M_* < M^* \le M$, let $\mathcal{I}_{M_*, M^*}$ be defined on $\Omega$ by

$$ \mathcal{I}_{M_*, M^*}(z) = \left( \sum_{m=M_*}^{M^*} \frac{|g_{z,d}^T v_m|^2}{\sigma_m} \right)^{-1} \quad (30) $$

**Remark 4.** From the identity (19), the Corollary 2 and the results of Section 4, our conclusion (that we do not formalize) is as follows: If the dimension of the matrix $\tilde{\Pi}_M \in \mathbb{R}^{M \times M}$ is large enough to feature a noise subspace $N_*$ with $M_*, M^* \in \mathbb{N}$, then $z \mapsto \mathcal{I}_{M_*, M^*}(z)$ takes large values in a region inside the inclusion $\mathcal{I}$, whereas this function is considerably smaller in $\Omega \setminus \bar{\mathcal{I}}$, and it smoothly changes from large to small values over the boundary of $\mathcal{I}$.

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. It is expected that the indicator function $\mathcal{I}_{M_*, M^*}$ introduced provides a good identification of the interior points, while an accurate representation of the geometry of the unknown objects will be achieved by the recourse to a thresholded version of the function.

## 5.2. 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 (30) rather than with the truncated Picard series (29). 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 [28] 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 $\delta > 0$, the original measurement matrix $\tilde{\Pi}_M$ is replaced by a noisy counterpart $\tilde{\Pi}_M^\delta$ such that $\|\tilde{\Pi}_M^\delta - \tilde{\Pi}_M\|_2 = \delta\|\tilde{\Pi}_M\|_2$ in the spectral norm. As a special noisy configuration we consider that

$$ \tilde{\Pi}_M^\delta = \tilde{\Pi}_M + \frac{1}{\sqrt{M}}\Lambda \qquad (31) $$

where $\Lambda$ is a real-valued random $M \times M$ matrix with independent and identically distributed entries described by Gaussian statistics with zero mean and standard deviation $\sigma_n$. In particular, with such a choice of perturbation, then for large $M$, one obtains at the first order $\delta = 2\sigma_n/\sigma_1$ in probability (see [45]). To simplify the estimates below, we assume moreover that the noise subspace is not modified by the noise,

$$ \mathrm{span}\{v_m^\delta, m = M_*, \dots, M\} = \mathrm{span}\{v_m, m = M_*, \dots, M\}. $$

This assumption is motivated by the fact that the signal subspace is only slightly perturbed by the noise, and that the noise subspace is by construction of the singular value decomposition part of the orthogonal complement of the signal subspace. Of course, in general one cannot exclude the possibility that noise subspace does change under the influence of noise, but this situation is beyond the scope of this paper.

Let us denote by $\sigma_m^\delta$, $v_m^\delta$ the singular values and vectors of the perturbed matrix $\tilde{\Pi}_M^\delta$. Under the above assumption the perturbed indicator function $\mathcal{I}_{M_*,M}^\delta$ verifies for all $z \in \Omega$ the inequalities

$$ \frac{1}{\max_{m=M_*,\dots,M}} \sigma_m^\delta \sum_{m=M_*}^M |g_{z,d}^\top v_m^\delta|^2 \le \mathcal{I}_{M_*,M}^\delta(z)^{-1} \le \frac{1}{\min_{m=M_*,\dots,M}} \sigma_m^\delta \sum_{m=M_*}^M |g_{z,d}^\top v_m^\delta|^2, $$

and one finally obtains the estimates

$$ \frac{\min_{m=M_*,\dots,M} \sigma_m^\delta}{\max_{m=M_*,\dots,M} \sigma_m} \mathcal{I}_{M_*,M}(z) \le \mathcal{I}_{M_*,M}^\delta(z) \le \frac{\max_{m=M_*,\dots,M} \sigma_m^\delta}{\min_{m=M_*,\dots,M} \sigma_m} \mathcal{I}_{M_*,M}(z). $$

Assuming the form (31), then for large $M$ and due to the rapidly decaying behavior of the singular values of the Neumann-to-Dirichlet operator, their noisy counterparts are such that $\max_{m=M_*,\dots,M} \sigma_m^\delta \approx \delta$, and on using that the eigenvalues are sorted in decreasing order, one has $\max_{m=M_*,\dots,M} \sigma_m = \sigma_{M_*} = \sigma_{M_\delta+1} \le \delta$, the above estimates can thus be rewritten as

$$ \frac{\sigma_M^\delta}{\delta} \mathcal{I}_{M_*,M}(z) \le \mathcal{I}_{M_*,M}^\delta(z) \le \frac{\delta}{\sigma_M} \mathcal{I}_{M_*,M}(z). \qquad (32) $$

The set of inequalities (32) provides a rough estimate since the term $\sigma_M^\delta$ has not been evaluated. The known literature concerning the characterization of the set of perturbed singular values and of the difference between the perturbed and the unperturbed right singular subspaces, surveyed in [46], deals with matrices of maximal rank strictly smaller
than $M$. Estimates are then known in such a case for deterministic or stochastic [45] perturbations. When the matrix of finite rank involves small singular values then higher order estimates have to be derived [47]. Therefore, a perturbation analysis in the present case of a matrix of rank $M$ with asymptotically small singular values remains a challenging problem and an open question.

If $\min_{m=M_*,...,M^*} \sigma_m^\delta \propto \delta$ with $M^* < M$, which can be conjectured from [45] when $1 \ll \text{rank}(\tilde{\Pi}_M) \ll M$ and in the asymptotic regime $M \gg 1$, then (32) shows that the indicator function (30) behaves linearly with respect to additive noise which makes the present approach particularly suited for the applications dealing with simulated or experimental data. Interestingly, this reconstruction strategy circumvents the issues related to the “inversion” of the compact Neumann-to-Dirichlet operator which is employed in the linear sampling method [7] and which necessitates the recourse to an adapted regularization scheme.

Following the interpretation of the indicator function (30) that is synthesized in Remark 4, given compact subsets $K_I \subset I$ and $K_{\Omega \setminus I} \subset \Omega \setminus I$ and on employing the inequalities (32) for an exterior point $z_e \in K_{\Omega \setminus I}$ and an interior point $z_i \in K_I$ one has that, if the reconstruction in the noise-free configuration is such that

$$ \frac{\delta^2}{\sigma_M \sigma_M^\delta} \mathcal{I}_{M_*,M}(z_e) \le \mathcal{I}_{M_*,M}(z_i), $$

then after perturbation, the set $I$ of inclusions can still be distinguished from the background since $\mathcal{I}_{M_*,M}^\delta(z_e) \le \mathcal{I}_{M_*,M}^\delta(z_i)$. Conversely, if the contrast provided by the unperturbed indicator function is initially too low, i.e. in the case where

$$ \frac{\delta^2}{\sigma_M \sigma_M^\delta} \mathcal{I}_{M_*,M}(z_i) \le \mathcal{I}_{M_*,M}(z_e), $$

then the noisy configuration leads to situations where one can obtain $\mathcal{I}_{M_*,M}^\delta(z_i) \le \mathcal{I}_{M_*,M}^\delta(z_e)$.

Two remarks can be made about these conclusions. On the one hand, given the properties of the measurement operator, the term $\delta^2/\sigma_M \sigma_M^\delta$ can be relatively large. Moreover the noise level is generally not known and obviously the unperturbed indicator function cannot be computed. On the other hand, from the Corollary 2 and the arguments developed in Section 4.3 regarding the approximation of the spectrum of $\Pi$, then in the absence of measurement noise and computation errors the indicator function is expected to be such that $\mathcal{I}_{M_*,M}(z_i)/\mathcal{I}_{M_*,M}(z_e) \gg 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 $\delta^2/\sigma_M \sigma_M^\delta$, then the reconstruction from data corrupted by the additive noise $\delta$ cannot be misinterpreted.

### 5.3. 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 the 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 $\Omega_h = [0; 1] \times [0; 1]$ with unit conductivity and which may contain an homogeneous inclusion (or a set thereof) of conductivity $\gamma_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 \in \partial\Omega_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 dependence to the mesh size $h$, the discrete version $\tilde{\Pi}_M$ of the relative Neumann-to-Dirichlet operator, is computed according to (21) where the entries $g_n(u_m - u_{1m})$ are consistently substituted by the orthonormal projection $f_n^T(u_m - u_{1m})$.

The indicator function (30) is implemented from the singular value decomposition of the matrix $\tilde{\Pi}_M$, the discrete version $g_{z,d}$ of the test function featured in Theorem 2 and a choice of dipole directions on a half-circle, i.e. $d_k = (\cos \theta_k; \sin \theta_k)$ where $\theta_k = (k-1)\pi/8$ and $k = 1, \dots, 8$. To circumvent the limitations due to the difficulties plaguing in the computation of the Green's function on the geometry considered, the formula (28) is employed. This latter expression requires the evaluation of the scalar products of the analytical solution (27) 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 $\Lambda_1$ which is given by the entries $f_n^T u_{1m}$. The sampling point $z$ varies on a regular sampling grid $\Omega_h^{prob}$ of size $n_{prob} \times n_{prob}$ points where $n_{prob} = 40$ provides a sufficient resolution. The contribution of each underlying dipole direction considered is summed and the optimal values of the indicator function are sought as $M_* = 20, \dots, 40$ and $M^* = 130, \dots, 140$, i.e.

$$ \hat{\mathcal{I}}'(z) = \max_{M_* = 20, \dots, 40} \max_{M^* = 130, \dots, 140} \sum_{k=1}^{8} \frac{\mathcal{I}_{M_*, M^*}(z)}{\max_{x \in \Omega_h^{prob}} \mathcal{I}_{M_*, M^*}(x)}. $$

Then the following graphs (Fig. 2 to 5) show, for clarity, the unitary indicator function

$$ \hat{\mathcal{I}}(z) = \frac{\hat{\mathcal{I}}'(z)}{\max_{x \in \Omega_h^{prob}} \hat{\mathcal{I}}'(x)}. $$

*Single inclusion.* The 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 $\gamma_I = 0.01$.
<figure><img src="image_3.png" /><figcaption>Figure 2. Identification of a single inclusion</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3. Identification of multiple inclusions</figcaption></figure>

*Non-convex inclusion.* A single L-shaped inclusion of conductivity $\gamma_I = 0.01$ is identified on the Figure 4.

<figure><img src="image_5.png" /><figcaption>Figure 4. Identification of L-shaped inclusion</figcaption></figure>
*Noisy data.* The noise-free configuration presented on the 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 (31) and the standard deviation $\sigma_n$ of the noisy matrix $\Lambda$ is chosen in order to achieve the desired value of $\delta$.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Identification of K=4 inclusions from noisy data</figcaption></figure>

*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, the figures 2 and 3 shows that a relatively small set of inclusions, characterized by different sizes or conductivities, can be located with a satisfying precision by the maxima of the function (30). 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 the Figure 3(c), the reconstruction provided is less contrasted between interior and exterior points, however four spikes clearly visible permit 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. The 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 value of the noise level, as shown on the Figure 5(a), and the reconstruction is still informative for the larger value $\delta = 0.05$ (Fig. 5(b)). Thus, the use of polluted spectral information has enabled a qualitatively satisfactory and stable identification. As expected, its quality decreases for larger value of the noise level (Fig 5(c)).

It is worth underlying as a final remark that the approach developed in this article is based on the *relative* Neumann-to-Dirichlet operator, which necessitates the 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 [26].

### 5.4. Experimental result

This paragraph presents two preliminary experimental results provided by V. Choquet and J. Alaterre in [48]. The experimental domain was a 20 × 28cm carbon-paper sheet (from www.pasco.com) with currents and potentials measured from standard laboratory sources and meters as proposed in [49]. The inclusion was a circular cut in the carbon-paper, centered at $x_0 = (20, 7)$ with radius $r = 3cm$, and the discretized version of the operator $\tilde{\Pi}_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 (30) the numerical approach described in Section 5.3 was used to estimate the discrete version of $g_{z,d}$.

<figure><img src="image_7.png" /><figcaption>Figure 6. Experiment 1: $M_* = 8$ and $M^* = 13$</figcaption></figure>

The figures 6 and 7 correspond to two different positionings of the electrodes that are represented by black diamonds on the right panels. Figures 6(a) and 7(a) plot the distribution of the singular values in the two experiments with the dashed red lines indicating the particular choice of the parameters $M_*$ and $M^*$. Based on the set of eight equidistributed dipole directions $d_k$, the following function is computed

$$ \hat{\mathcal{I}}'_{M_*, M^*}(z) = \sum_{k=1}^{8} \frac{\mathcal{I}_{M_*, M^*}(z)}{\max_{x \in \Omega_h^{prob}} \mathcal{I}_{M_*, M^*}(x)}, $$

and then the corresponding unitary indicator function

$$ \hat{\mathcal{I}}_{M_*, M^*}(z) = \frac{\mathcal{I}'_{M_*, M^*}(z)}{\max_{x \in \Omega_h^{prob}} \mathcal{I}'_{M_*, M^*}(x)}, $$

is shown on the figures 6(b) and 7(b). The parameters $M_*$ and $M^*$ are manually chosen from the evaluation of the signal and noise trends in the displayed distribution of the singular values. These satisfying initial results involving a single large inclusion prove
the feasibility of the method from an experimental point of view. More measurements
points and precision would contribute to increase the accuracy of the reconstruction.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Experiment 2: $M_* = 7$ and $M^* = 12$</figcaption></figure>

**6. 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-to-measurement operator and which commonly span its signal subspace, the approach developed in this article 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 article 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 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.

Acknowledgments

This work was initiated by Thomas Coquet and Thomas Jaravel during the course "Inverse Problems" EA MEC574 in the Department of Mechanics at the Ecole Polytechnique, Palaiseau, France. It has been pursued by Vincent Choquet and Julien Alaterre who are kindly acknowledged for providing us with the experimental results of Section 5.4. The first author would like to thank Guillaume Bal and François Monard for fruitful interactions and helpful discussions.

References

[1] A. Kirsch. *An introduction to the mathematical theory of inverse problems*. Springer, 2nd edition, 2011.

[2] G. Alessandrini. Stable determination of conductivity by boundary measurements. *Appl. Anal.*, 27:153–172, 1988.

[3] G. Alessandrini and S. Vessella. Lipschitz stability for the inverse conductivity problem. *Adv. in Appl. Math.*, 35:207–241, 2005.

[4] M. Cheney, D. Isaacson, and J. C. Newell. Electrical impedance tomography. *SIAM Review*, 41:85–101, 1999.

[5] L. Borcea. Electrical impedance tomography. *Inverse Problems*, 18:R99-R136, 2002.

[6] W. R. B. Lionheart. Developments in EIT reconstruction algorithms: pitfalls, challenges and recent developments. *Physiol. Meas.*, 25:125–142, 2004.

[7] F. Cakoni and D. Colton. *Qualitative methods in inverse scattering theory*. Springer-Verlag, Berlin, 2006.

[8] M. Hanke and M. Brühl. Recent progress in electrical impedance tomography. *Inverse Problems*, 19:S65–S90, 2003.

[9] A. Kirsch and N. Grinberg. *The Factorization Method for Inverse Problems*. Oxford University Press, New York, 2008.

[10] R. Potthast. A survey on sampling and probe methods for inverse problems. *Inverse Problems*, 22:R1–R47, 2006.

[11] M. Ikehata. Reconstruction of inclusion from boundary measurements. *J. Inv. Ill-Posed Problems*, 10:37–65, 2002.

[12] L. Jackowska-Strumillo, J. Sokolowski, and A. Zochowski. Topological optimization and inverse problems. *Computer Assisted Mechanics and Engineering Sciences*, 10(2):163–176, 2002.

[13] M. Bonnet and B. B. Guzina. Sounding of finite solid bodies by way of topological derivative. *Int. J. Num. Meth. in Eng.*, 61:2344–2373, 2004.

[14] A. J. Devaney. Super-resolution processing of multi-static data using time reversal and music. *Northeastern University Preprint*, 2000.

[15] A. J. Devaney. Time reversal imaging of obscured targets from multistatic data. *IEEE Trans. Antennas Propag.*, 53:1600–1610, 2005.

[16] P. Hähner. An inverse problem in electrostatics. *Inverse Problems*, 15:961–975, 1999.
[17] M. Brühl and M. Hanke. Numerical implementation of two noniterative methods for locating inclusions by impedance tomography. *Inverse Problems*, 16:1029–1042, 2000.

[18] M. Brühl. Explicit characterization of inclusions in electrical impedance tomography. *SIAM J. Math. Anal.*, 32:1327–1341, 2001.

[19] M. Brühl, M. Hanke, and M. S. Vogelius. A direct impedance tomography algorithm for locating small inhomogeneities. *Numer. Math.*, 93:635–654, 2003.

[20] R. Kress and L. Kühn. Linear sampling methods for inverse boundary value problems in potential theory. *Appl. Numer. Math.*, 43:161–173, 2002.

[21] M. Cheney. The linear sampling method and the MUSIC algorithm. *Inverse Problems*, 17:591–595, 2001.

[22] M. Hanke and B. Schappel. The factorization method for electrical impedance tomography in the half-space. *SIAM J. Appl. Math.*, 68:907–924, 2008.

[23] B. Gebauer and N. Hyvönen. Factorization method and irregular inclusions in electrical impedance tomography. *Inverse Problems*, 23:2159–2170, 2007.

[24] N. Hyvönen, H. Hakula, and S. Pursiainen. Numerical implementation of the factorization method within the complete electrode model of electrical impedance tomography. *Inverse Probl. Imaging*, 1:299–317, 2007.

[25] A. Lechleiter, N. Hyvönen, and H. Hakula. The factorization method applied to the complete electrode model of impedance tomography. *SIAM J. Appl. Math.*, 68:1097–1121, 2008.

[26] B. Harrach and J. K. Seo. Detecting inclusions in electrical impedance tomography without reference measurements. *SIAM J. Appl. Math.*, 69:1662–1681, 2009.

[27] D. R. Luke and A. J. Devaney. Identifying scattering obstacles by the construction of nonscattering waves. *SIAM J. Appl. Math.*, 68:271–291, 2007.

[28] T. Arens, A. Lechleiter, and D. R. Luke. MUSIC for extended scatterers as an instance of the factorization method. *SIAM J. Appl. Math.*, 70:1283–1304, 2009.

[29] M.'B. Fares, S. Gratton, and P. L. Toint. SVD-tail: a new linear-sampling reconstruction method for inverse scattering problems. *Inverse Problems*, 25:095013, 2009.

[30] A. P. Calderón. On an inverse boundary value problem. *Seminar on Numerical Analysis and its Application to Continuum Physics (Soc. Brasileira de Matemática, Rio de Janeiro)*, pages 65–73, 1980.

[31] J. Sylvester and G. Uhlmann. A global uniqueness theorem for an inverse boundary value problem. *Ann. Math.*, 125:153–169, 1987.

[32] V. Isakov. *Inverse problems for partial differential equations*. Springer, 1998.

[33] K. Astala and L. Päivärinta. Calderón's inverse conductivity problem in the plane. *Ann. Math.*, 163:265–299, 2006.

[34] D. J. Cedio-Fengya, S. Moskow, and M. S. Vogelius. Identification of conductivity imperfections of small diameter by boundary measurements. Continuous dependence and computational reconstruction. *Inverse Problems*, 14:553–595, 1998.

[35] Martin Costabel, Monique Dauge, and Serge Nicaise. Analytic Regularity for Linear Elliptic Systems in Polygons and Polyhedra. 54 pages.

[36] W. McLean. *Strongly Elliptic Systems and Boundary Integral Operators*. Cambridge University Press, Cambridge, UK, 2000.

[37] P. Monk. *Finite Element Methods for Maxwell's Equations*. Oxford Science Publications, Oxford, 2003.

[38] S. Sauter and C. Schwab. *Boundary element methods*. Springer, 1. edition, 2011.

[39] J. Saranen and G. Vainikko. *Periodic integral and pseudodifferential equations with numerical approximation*. Springer, 2002.

[40] N. Hyvönen. Approximating idealized boundary data of electric impedance tomography by electrode measurements. *Mathematical Models and Methods in Applied Sciences*, 19:1185–1202, 2009.

[41] T. Kato. *Perturbation theory for linear operators*. Springer, repr. of the 1980 edition, 1995.
[42] G. W. Stewart. Perturbation theory for the singular value decomposition. *SVD and Signal Processing, II: Algorithms, Analysis and Applications*, pages 99–109, 1991.

[43] A. Lechleiter. A regularization technique for the factorization method. *Inverse Problems*, 22:1605–1625, 2006.

[44] M. A. Fischler and R. C. Bolles. Random sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography. *Communications of the ACM*, 24:381–395, 1981.

[45] J. Garnier. Use of random matrix theory for target detection, localization, and reconstruction. *Contemporary Mathematics*, 548:1–19, 2011.

[46] G. W. Stewart and J.-G. Sun. *Matrix perturbation theory*. Academic Press, New York, 1990.

[47] M. Stewart. Perturbation of the SVD in the presence of small singular values. *Linear Algebra and its Applications*, 419:53–77, 2006.

[48] V. Choquet and J. Alaterre. Localisation d'inclusions par tomographie d'impedance électrique. Rapport EA MEC574, Ecole Polytechnique, 2011.

[49] R. Lazarovitch, D. Rittel, and I. Bucher. Experimental crack identification using electrical impedance tomography. *NDT & E International*, 35:301–316, 2002.