Houssem Haddar

Giovanni Migliorati
email: giovanni.migliorati@gmail.com

Numerical analysis of the Factorization Method for EIT with piece-wise constant uncertain background

Keywords: AMS classification scheme numbers: 35R60, 35R30, 65M32 Inverse Problems, EIT, Factorization Method, uncertain background

We extend the Factorization Method for Electrical Impedance Tomography to the case of background featuring uncertainty. We first describe the algorithm for known but inhomogeneous backgrounds and indicate expected accuracy from the inversion method through some numerical tests. Then we develop three methodologies to apply the Factorization Method to the more difficult case of piece-wise constant but uncertain background. The first one is based on a recovery of the background through an optimization scheme and is well adapted to relatively low dimensional random variables describing the background. The second one is based on a weighted combination of the indicator functions provided by the Factorization Method for different realizations of the random variables describing the uncertain background. We show through numerical experiments that this procedure is well suited to the case where many realizations of the measurement operators are available. The third strategy is a variant of the previous one when measurements for the inclusion-free background are available. In that case, a single pair of measurements is sufficient to achieve comparable accuracy to the deterministic case.

Introduction

The problem of Electrical Impedance Tomography (EIT) arises in many applied contexts. It leads to well established operative procedures, e.g. in geophysics, nondestructive testing or imaging applications, while in other fields its usage is still experimental e.g. medicine, [START_REF] Muller | Linear and Nonlinear Inverse Problems with Practical Applications[END_REF][START_REF] Borcea | Electrical impedance tomography[END_REF][START_REF] Cheney | Electrical impedance tomography[END_REF][START_REF] Seppänen | State estimation with fluid dynamical evolution models in process tomography-an application to impedance tomography[END_REF][START_REF] Guo | A review of electrical impedance techniques for breast cancer detection[END_REF].

For some key theoretical results related to uniqueness and conditional stability of the EIT inverse problem, we refer to [START_REF] Astala | Calderón's inverse conductivity problem in the plane[END_REF][START_REF] Brown | Uniqueness in the inverse conductivity problem for conductivities with 3/2 derivatives in L p , p > 2n[END_REF][START_REF] Alessandrini | Stable determination of conductivity by boundary measurements[END_REF][START_REF] Lee | Determining anisotropic real-analytic conductivities by boundary measurements[END_REF] and references therein. On the numerical side, many mathematical models and inversion algorithms for EIT were proposed in the literature, trying to take into account as many physical phenomena as possible, see [START_REF] Muller | Linear and Nonlinear Inverse Problems with Practical Applications[END_REF] and the references therein for a complete overview. In this work we address the so-called Factorization Method (FM) as introduced in [START_REF] Kirsch | The factorization method for inverse problems[END_REF] applied to the Continuous Model (CM) (in contrast with the so-called complete electrode model [START_REF] Lechleiter | The factorization method applied to the complete electrode model of impedance tomography[END_REF][START_REF] Hyvönen | Numerical implementation of the factorization method within the complete electrode model of electrical impedance tomography[END_REF]) in the context of EIT featuring an inhomogeneous isotropic background. We shall in particular address the case where the inhomogeneous background cannot be considered as a small perturbation of the homogeneous one (as in [START_REF] Harrach | Detecting inclusions in electrical impedance tomography without reference measurements[END_REF]) and more importantly consider the case where this background is piecewise constant but with uncertain values on the conductivity (and not on the spatial distribution). The latter configuration is motivated by medical imaging experiments since the conductivity support may be known from the use of other imaging modalities, while the physical electrical parameters cannot be exactly determined.

Most of the works on the FM in the literature treat the case of homogeneous background, e.g. [START_REF] Azzouz | The factorization method for electrical impedance tomography data from a new planar device[END_REF][START_REF] Hanke | The factorization method for electrical impedance tomography in the half-space[END_REF]. Our work is mostly related to numerical issues associated with EIT in inhomogeneous and uncertain background. To begin with, we propose a numerical scheme to solve the dipole-like Neumann boundary-value problem when the background coefficient is inhomogeneous and deterministic, and use this scheme to design an efficient implementation of the FM algorithm for inhomogeneous but deterministic background. We then discuss the case where the background is piecewise constant with known spatial distribution but unknown parameter values. We propose three variants of the FM to cope with this configuration. In the first algorithm, by means of an optimization scheme motivated by the structure of the sampling operator and the outcome of the FM, we simultaneously recover the background parameters and the shape of the inclusion. This algorithm is well suited for low dimensional configurations of the parameter space, and can be applied to moderately high dimensional configurations by means of acceleration with polynomial approximation techniques [START_REF] Babuška | Galerkin finite element approximations of stochastic elliptic partial differential equations[END_REF][START_REF] Cohen | Analytic regularity and polynomial approximation of parametric and stochastic elliptic PDE's[END_REF][START_REF] Migliorati | Analysis of the discrete l 2 projection on polynomial spaces with random evaluations[END_REF]. In addition, an alternative approach is proposed, where the optimization scheme is replaced by a weighting of the FM indicator function with a misfit indicator for the background. This procedure requires an extensive sampling of the parameter space which would be very expensive. However, in the case where many independent realizations of the measurement operator are available, we describe how this strategy can be made efficient. The third algorithm deals with the case where paired measurements are available, namely measurements for the inclusion free background and measurements for the medium with inclusion. In this case, a uniform weight for the FM associated with random sampling of the parameter spaces provides an effective indicator function.

The outline of the paper is the following: in Section 2.1 we present the problem of EIT and the notation. In Section 2.2 we introduce the FM in deterministic inhomogeneous background and the numerical scheme proposed to solve the dipolelike singular problem. Then in Section 2.3 we describe the regularization technique employed. In Section 3 we introduce the problem of EIT in random media, starting from the parametrization of the background coefficient by random variables. A justification of how randomness affects the spectrum of the measurement operator is provided. Moreover, two situations are identified depending on how the measurements are collected, i.e. arbitrary or paired measurements. In Section 3.1 two variants of the FM are proposed, in the case of arbitrary measurements. In Section 3.2 a variant of the FM is proposed in the case of paired measurements. Sections 4, 5, 6 illustrate the performance of the variants of the FM proposed through several testcases. Finally, in Section 7 we draw some conclusions.

The EIT problem and the Factorization Method

Mathematical formulation of the EIT problem

Consider a bounded Lipschitzian domain B ⊂ R 2 and its subset D ⊂ B. We assume that D is a union of possibly disjoint Lipschitzian domains, each one with positive measure, and that B \ D is connected. The domain B represents the background medium, featured by the physical coefficient σ B . The domain D represents an inclusion, displaced somewhere inside B. The inclusion is characterized by unknown shape and unknown value of its physical coefficient σ D . We shall assume that σ B and σ D are real valued functions. Now consider the coefficient σ(x) ∈ L ∞ (B) defined as

σ(x) = σ B (x), in B \ D, σ D (x), in D, (1)
and define the functional spaces L2 (∂B), H1 (B), H 1 2 (∂B) and H-1 2 (∂B) as the subspaces of respectively L 2 (∂B), H 1 (B), H 1 2 (∂B) and H -1 2 (∂B) of functions with zero (boundary) mean value on ∂B.

Given g ∈ H-

∇ • (σ(x)∇u) = 0, in B, σ(x)∇u • ν = g, on ∂B, (2)
where ν denotes the outward normal to ∂B. It is well known that problem (2) has a unique solution if σ ∈ L ∞ (B) and is positive definite. In the case of the continuous EIT model, one prescibes a current pattern g and measures the generated potential f = u| ∂B on the boundary ∂B. The final aim is to retrieve some information about D, e.g. its location inside B and its shape. The operator Λ that maps currents g into potentials f is the NtDm associated to problem (2) with coefficient [START_REF] Muller | Linear and Nonlinear Inverse Problems with Practical Applications[END_REF]. It is a continuos operator from H-1 2 (∂B) to H 1 2 (∂B).

Denote also by Λ 0 the NtDm for the problem (2) with the coefficient

σ(x) = σ B (x), in B, (3)
and by u 0 the solution of the same problem. In this way u 0 represents the potential generated by the incoming current g in the domain B, when the inclusion D is not present. Let f 0 be the corresponding measured potential f 0 = u 0 | ∂B = Λ 0 g. Both the operators Λ and Λ 0 are weakly compact, when restricted on L2 (∂B), i.e. Λ : L2 (∂B) → L2 (∂B) and Λ 0 : L2 (∂B) → L2 (∂B). Moreover, we define the operator Λ

Λ := Λ -Λ 0 : L2 (∂B) → L2 (∂B), (4)
that has the remarkable property of being strongly compact (eigenvalues decay exponentially fast: see for instance (32) ).

We will need also the Green function N (•, ξ) ∈ L2 (B), which is a solution of the following Neumann boundary-value problem with a Dirac source δ ξ centered in ξ ∈ B:

∇ 1 • (σ B (x)∇ 1 N (x, ξ)) = -δ ξ (x), x ∈ B, ν • σ B (x)∇ 1 N (x, ξ) = -1 |B| , x on ∂B. (5)
Note that ( 5) does not embed any information about the inclusion D. Denote by p a two-dimensional versor (i.e. |p| = 1). In the sequel we will often need the test function

ψ(x, ξ, p) := p • ∇ 2 N (x, ξ), x ∈ B, (6)
as well as its restriction on ∂B

l p ξ = ψ(x, ξ, p) x on ∂B . (7)
Moreover, denoting with {e 1 , e 2 } an orthonormal basis of R 2 , we define

l k ξ = l e k ξ , k = 1, 2.
Remark 1 When the domain B is a circle with radius R there is an explicit formula for the solution N (x, ξ) of the Neumann problem [START_REF] Guo | A review of electrical impedance techniques for breast cancer detection[END_REF] with σ B ≡ 1 (see e.g. [START_REF] Brühl | Explicit characterization of inclusions in electrical impedance tomography[END_REF]):

N (x, ξ) = - 1 2π log |x -ξ| + log R |x| x - |x| R ξ + log R π . (8)
In this case we have an explicit formula for the evaluation of ψ defined in [START_REF] Astala | Calderón's inverse conductivity problem in the plane[END_REF] on the boundary of the unitary circle :

ψ(x, ξ, p) ∂B = 1 π p • (ξ -x) |ξ -x| 2 , for |x| = 1. (9)

The Factorization Method with deterministic inhomogeneous background

In this section we outline the basis of the FM applied to EIT and some numerical issues related to the case of inhomogeneous background. The main result that motivates the FM algorithm is Theorem [START_REF] Muller | Linear and Nonlinear Inverse Problems with Practical Applications[END_REF]. It provides a necessary and sufficient criterion to localize a point inside the support of the inclusions. Remark 2 shows how this criterion becomes operative. This theorem can be proved following the lines of the proof of [10, Theorem 6.8] that treats the case of constant σ B (see also [START_REF] Brühl | Explicit characterization of inclusions in electrical impedance tomography[END_REF]). The only point that requires attention is the use of unique continuation arguments for equations of type (2), and this is why a regularity assumption on σ B is required.

Assumption 1 Assume σ and σ B are bounded and positive definite functions on B and that σ B is piecewise Lipschitz continuous and with Lipschitz discontinuity interfaces. Moreover, denoting by

σ m B = min x∈D σ B (x) and σ M B = max x∈D σ B (x),
then σ B and σ D satisfy

σ D (x) < σ m B -δ or σ M B < σ D (x) -δ a.
e. in D, for some positive constant δ.

Theorem 1 Assume that assumption (1) holds. Then the operator Λ : L2 (∂B) → L2 (∂B) is self adjoint and injective. Moreover, for any versor p

l p ξ ∈ Range(| Λ| 1/2 ) ⇐⇒ ξ ∈ D.
Remark 2 Denote by {λ j , w j : j ∈ J} the spectral system of the self adjoint and compact operator Λ. Then

w ∈ Range(| Λ| 1/2 ) ⇐⇒ ∞ j=1 |(w, w j ) L 2 (∂B) | 2 |λ j | < ∞. (10)
The main difficulty in implementing the FM for inhomogeneous background is the precise evaluation of l p ξ that can only be done in the general case by numerically solving problem [START_REF] Guo | A review of electrical impedance techniques for breast cancer detection[END_REF]. To overcome the singularity in the forcing term, we resort to the fundamental solution

φ(•, ξ) = - 1 2πσ B (ξ) log • -ξ , x ∈ R 2 , (11)
of the problem

-∇ • σ B (ξ)∇φ(•, ξ) = δ ξ , in R 2 . ( 12
)
Since the singularity at ξ in problem ( 5) is of the same kind as in problem [START_REF] Hyvönen | Numerical implementation of the factorization method within the complete electrode model of electrical impedance tomography[END_REF], we can restrict problem (5) in a small neighbourhood O(ξ) ⊂ R 2 of ξ where

σ B (x) ≈ σ B (ξ), ∀ x ∈ O(ξ),
and approximate the solution N (x, ξ) of problem (5) near the singularity in ξ as

N (x, ξ) ≈ φ(x, ξ), ∀ x ∈ O(ξ). ( 13
)
Then it is possible to write a nonsingular problem for the difference

ϕ N (•, ξ) = N (•, ξ) -φ(•, ξ), (14)
plugging ( 14) in [START_REF] Guo | A review of electrical impedance techniques for breast cancer detection[END_REF]. The quantity ψ(•, ξ, p) defined in ( 6) is a solution of problem (5) for any p. Moreover, in ψ the dependece of N (•, •) on the second argument is smooth.

So we can exploit this regularity to derive a numerical scheme that given ξ ∈ B directly computes ∇ 2 N (•, ξ), that is needed to obtain l p ξ through ψ and apply Theorem 1. To this aim we differentiate [START_REF] Azzouz | The factorization method for electrical impedance tomography data from a new planar device[END_REF] with respect to the second argument, and obtain

∇ 2 ϕ N (•, ξ) = ∇ 2 N (•, ξ) -∇ 2 φ(•, ξ). ( 15
)
The derivatives of φ w.r.t. ξ are ∂φ(•, ξ)

∂ξ k = 1 2πσ B (ξ) 2 ∂σ B (ξ) ∂ξ k log |x -ξ| + x k -ξ k 2πσ B (ξ)|x -ξ| 2 , k = 1, 2,
and since φ is analytical when x = ξ then the mixed derivatives coincide:

∂ 2 φ(x, ξ) ∂x k ∂ξ k = 1 2πσ B (ξ) 2 ∂σ B (ξ) ∂ξ 1 x k -ξ k |x -ξ| 2 + |x -ξ| 2 -2(x k -ξ k ) 2 2πσ B (ξ)|x -ξ| 4 , ∂ 2 φ(x, ξ) ∂x ℓ ∂ξ k = 1 2πσ B (ξ) 2 ∂σ B (ξ) ∂ξ 1 x ℓ -ξ ℓ |x -ξ| 2 - (x k -ξ k ) (x ℓ -ξ ℓ ) πσ B (ξ)|x -ξ| 4
, for k, ℓ = 1, 2 and ℓ = k. Denoting by superscript the partial derivative, we have two problems for the unknowns

ϕ k N = ∂ϕ N /∂ξ k , k = 1, 2:    -∇ x • σ B (x)∇ 1 ϕ k N (x, ξ) = -∇ x • σ B (ξ) -σ B (x) Φ k (x, ξ) , x in B, σ B (x) ∇ 1 ϕ N (x, ξ) • ν = -σ B (x) Φ k (x, ξ) • ν -1 |∂B| , x on ∂B, (16)
with

Φ k (x, ξ) = ∂ 2 φ(x, ξ) ∂x 1 ∂ξ k , ∂ 2 φ(x, ξ) ∂x 2 ∂ξ k , ∀ (x, ξ) ∈ B × B : x = ξ.
All the second derivatives of φ are singular in x = ξ, and behave like 1/|x -ξ| 2 . Therefore if σ B (•) is piecewise Lipschitz then the solution of ( 16) is in H 1 (B). Hence the variational formulation of ( 16) is suited for standard finite element discretizations.

Once the solution of ( 16) has been computed, then ∇ 2 N (•, ξ) can be recovered from [START_REF] Hanke | The factorization method for electrical impedance tomography in the half-space[END_REF]. Note that when σ B is piece-wise constant and ξ does not fall where it jumps, the function ζ(x, ξ) = σ B (ξ) -σ B (x) vanishes as x approaches ξ.

Description of the algorithm for the deterministic setting

Now we present our operative criterion to implement the range test. It is based on the spectral decomposition of the operator | Λ| 1/2 and on the use of an appropriate regularization. We focus mainly on Tikhonov Regularization (TR), although also a straightforward application of the Picard Criterion (23) can give good results when the background is homogeneous (see Remark 3 below). To simplify the notation we set M = | Λ| 1/2 . To check whether a given l k ξ belongs to Range(M ) we have to solve the problem

M g k ξ = l k ξ , (17)
which demands for regularization. The keypoint is

ξ ∈ D ⇐⇒ g k ξ L 2 (∂B) -1 = 0.
Numerical analysis of the Factorization Method for EIT with uncertain background 7

The TR of [START_REF] Cohen | Analytic regularity and polynomial approximation of parametric and stochastic elliptic PDE's[END_REF] reads

α + M * M g k ξ = M * l k ξ . (18)
We choose the Fourier basis to discretize the operator M . Since we discretize l k ξ and g k ξ in [START_REF] Cohen | Analytic regularity and polynomial approximation of parametric and stochastic elliptic PDE's[END_REF] over the same orthogonal basis, then the discretization of M is a square matrix.

In the sequel we denote by σ i the singular values of M and by u i and v i the corresponding left and right singular vectors. Both the sets {u i } ∞ i=1 and {v i } ∞ i=1 form an orthonormal basis of L 2 (∂B). The regularized solution of problem [START_REF] Migliorati | Analysis of the discrete l 2 projection on polynomial spaces with random evaluations[END_REF] is

g k ξ = ∞ i=1 σ i α + σ 2 i (l k ξ , u i )v i . (19)
We choose the regularization parameter α using Morozov principle, i.e. imposing

M g k ξ -l k ξ L 2 (∂B) = δ, (20)
with the parameter δ > 0 to be in some way related to the accuracy of the measurements. We pick

δ = γσ 1 , (21)
being σ 1 the largest singular value of M , and γ a given threshold. The terms l k ξ depends on ξ, so to ensure that ( 20) is uniformly satisfied we normalize l k ξ as lk ξ = l k ξ / l k ξ L 2 (∂B) . The optimal value of α is computed plugging [START_REF] Brühl | Explicit characterization of inclusions in electrical impedance tomography[END_REF] into the nonlinear equation [START_REF] Gebauer | Factorization method and irregular inclusions in electrical impedance tomography[END_REF], and solving with respect to α. To show the inversion results, we display the isolines of the indicator function

C(ξ) = log v 1 (ξ) + v 2 (ξ) -1 , (22)
where v k (ξ) := g k ξ L 2 (∂B) . We can summarize the algorithm by the following:

Algorithm 1 The Factorization Method in inhomogeneous deterministic background Sample the region of B to be probed with a set of points P = {ξ j } P j=1 for ξ in the set P do solve problem (5) to find its solution N (•, ξ), compute l k ξ from N (•, ξ) using ( 7), use l k ξ to compute the indicator C(ξ) in [START_REF] Haddar | Conformal mapping and impedance tomography[END_REF], plot ξ → C(ξ). end for Remark 3 (Picard criterion) Another criterion to test the range condition is analyzed in [START_REF] Gebauer | Factorization method and irregular inclusions in electrical impedance tomography[END_REF]. It employs the indicator function:

I k (ξ) = 1 l k ξ 2 L 2 (∂B) ∞ i=1 (l k ξ , u i ) 2 L 2 (∂B) 2 |σ i | , ( 23
)
and it reads ξ ∈ D ⇐⇒ I k (ξ) < ∞. Note that the normalization of l k ξ is already embedded in [START_REF] Siltanen | An implementation of the reconstruction algorithm of A Nachman for the 2D inverse conductivity problem[END_REF]. We refer to this criterion as Picard Criterion (PC). It does not require to solve problem [START_REF] Cohen | Analytic regularity and polynomial approximation of parametric and stochastic elliptic PDE's[END_REF] to find its regularized solution. However, as discussed in [START_REF] Haddar | Numerical analysis of the factorization method for electrical impedance tomography in inhomogeneous medium[END_REF], the regularization technique TR is more accurate than PC when the background is not homogeneous, above all in presence of many inclusions. Hence, in this paper we present only the results obtained using the TR technique. See also [START_REF] Haddar | Numerical analysis of the factorization method for electrical impedance tomography in inhomogeneous medium[END_REF] for a comparison between TR and PC for several numerical examples.

We consider backgrounds that are piece-wise constant in r regions {R i } r i=1 that are a partition of unity of the computational domain B = r i=1 R i , see e.g. Fig. 3. Denote by m the r-dimensional vector of positive numbers that specifies the coefficient value in each region R i of the domain B. The background coefficient σ B is defined as

σ B (x) = r i=1 m i I R i (x), x ∈ B, (24)
where I X (x) denotes the indicator function of the subdomain X ⊆ B. If m is the unitary vector then we obtain a homogeneous deterministic coefficient. Since σ B is a diffusion coefficient, it has to be positive: therefore we assume that m i > 0, ∀ i. More general deterministic inhomogeneous coefficients σ B , nonlinearly depending on x, have been analyzed in [START_REF] Haddar | Numerical analysis of the factorization method for electrical impedance tomography in inhomogeneous medium[END_REF]. The effect of noise has been analyzed in [START_REF] Haddar | Numerical analysis of the factorization method for electrical impedance tomography in inhomogeneous medium[END_REF] as well.

The problem of EIT in random background

A natural issue that arises in the mathematical modeling of physical phenomena concerns the uncertainty in the input data. There are mainly two interpretations on how the presence of uncertainty in the model could be characterized: the epistemic and the aleatory interpretation. We introduce them directly in the context of EIT and the applications that we have in mind.

In the EIT framework the most natural choice to incorporate in the model the uncertainty affecting the background media is to consider the background coefficient σ B as a random variable. In the epistemic uncertainty, we can imagine that the coefficient that models the background material is deterministic but unkwnown, because of the lack of knowledge or because of inaccurate measure instruments. Therefore we can model it as a random variable, although there are also other ways to cope with epistemic uncertainty. In the aleatory uncertainty, the coefficient is really random, according to a prescribed probability density. We could think of its realizations coming from a stationary or from a nonstationary stochastic process.

In general, the treatment of uncertainty through random variables distributed according to given probability densities is known as probabilistic approach. For any positive integer N , we introduce the N -dimensional parameter space Γ ⊆ R N , the probability density ρ : Γ → R and the Hilbert space L 2 (Γ, dρ) equipped with the inner

product u, v = Γ u(y) v(y) dρ(y).
Denote by Y ∼ ρ a multidimensional random variable, and by y ∈ Γ any of its realizations. We denote by L 2 (Γ, L ∞ (B)) the Bochner space of functions that take values into Γ and are L ∞ (B)-valued. Then we define the random diffusion coefficient

σ B ∈ L 2 (Γ, L ∞ (B)) as σ B (•, y) ∈ L ∞ (B), ∀ y ∈ Γ. ( 25
)
Moreover, we assume a smooth dependence of σ B on y:

Assumption 2 ∀ y ∈ Γ the map y → σ(•, y) is assumed to be smooth (analytical).

We will consider two parametrizations of the random coefficient σ B introduced in (25). Again, the spatial background is piece-wise constant in r regions {R i } r i=1 that are a partition of unity of the computational domain B = r i=1 R i . The randomness concerns the value of the coefficient in each region R i . We associate each component of the Ndimensional random variable Y to a region, hence N = r.

The first parametrization is suited for uncertainty analysis, and is an extension of (24) to the case of random background:

σ B (x, Y) = r i=1 m i + d i • Y i I R i (x), x ∈ B, Y ∈ Γ = [-1, 1] r . (26)
The vectors m and d quantify the mean and dispersion of the coefficient in each region, and

m 1 = 1, d 1 = 0. The range of variation in each region is [m i -d i , m i + d i ],
and we have to assume that m i -d i > 0, ∀ i, as we did for [START_REF] Lagarias | Convergence properties of the Nelder-Mead simplex method in low dimensions[END_REF].

The second parametrization is suited for large variations of the background coefficient σ B , and we parametrize it with the exponential model

σ B (x, Y) = r i=1 10 Y i I R i (x), x ∈ B, Y ∈ Γ = [-1, 1] r , (27)
so that the coefficient jumps by two orders of magnitude in each one of the r regions.

Because of the randomness in σ B , the diffusion coefficient (1) of problem (2) becomes random as well, and the solution u = u(x, Y) of (2) belongs to the space H1 (B) ⊗ L 2 (Γ, dρ) [START_REF] Babuška | Galerkin finite element approximations of stochastic elliptic partial differential equations[END_REF][START_REF] Cohen | Analytic regularity and polynomial approximation of parametric and stochastic elliptic PDE's[END_REF][START_REF] Migliorati | Analysis of the discrete l 2 projection on polynomial spaces with random evaluations[END_REF]. In the same way the Neumann-to-Dirichlet maps [START_REF] Guo | A review of electrical impedance techniques for breast cancer detection[END_REF] with the random coefficient (25). For the same reason, the quantities (7) depend now on any realization y ∈ Γ of Y,

Λ : L2 (∂B) ⊗ L 2 (Γ, dρ) → L2 (∂B) ⊗ L 2 (Γ, dρ) and Λ 0 : L2 (∂B) ⊗ L 2 (Γ, dρ) → L2 (∂B)⊗L 2 (Γ, dρ) become random themselves. Also the Green function N = (x, ξ, y) ∈ D ⊗ D ⊗ Γ → N (x, ξ, y) : D ⊗ D ⊗ Γ → L2 (B) ⊗ L 2 (Γ, dρ) depends on the random variable Y because it solves problem
l p ξ = l p ξ (y), ∀ y ∈ Γ. ( 28
)
We define a measurement of the random operator Λ (or Λ 0 ) its observation for a single realization of the random variable Y. In practice the measurements of Λ come from experiments, while the measurements of Λ 0 could be simulated numerically. We will always assume that it is not possible to recover the value of the realization of Y that generated the measurement of Λ, because this would reduce the inverse problem in random background to an inverse problem with deterministic background. Moreover, we need also the following assumption: Assumption 3 We assume to be able to evaluate the operator Λ for any realization of the random variable Y. To evaluate the operator Λ means that we are able to collect all the voltages corresponding to the current patterns in the Fourier basis, up to a certain accuracy.

This assumption is completely fine in the epistemic uncertainty or in the stationary aleatory uncertainty, but could not be the case in the nonstationary aleatory uncertainty. Now we discuss the crucial point of how to extend to random background the operator Λ defined in (4) for inhomogeneous but deterministic background. The main question is: how do we collect the measurements of the random operators Λ and Λ 0 ? Denote by {y m } and { y t } the sets of unknown independent realizations of the random variable Y that generated the corresponding measurements of the random operator Λ and Λ 0 , respectively. Even if the values of the realizations are unknown, at least two different situations can be identified, depending on how the measurements are collected:

• arbitrary measurements: the measurements of Λ and Λ 0 are collected separately, i.e. the sets {y m } m and { y t } t have nothing to do each other. In the sequel we need only to assume that the realizations in each set are independent only among those belonging to the same set.

• paired measurements: the measurements of Λ and Λ 0 are collected in pair, i.e. the mth measurement of Λ and Λ 0 is collected with the same (but unknown) realization y m = y m of the random variable Y .

The case of arbitrary measurements is the general one, and it arises in both epistemic and aleatory uncertainty. The case of paired measurements arises in epistemic uncertainty. For example, in crack testing: several mechanical objects are built together with the same internal but unknown properties. The measurements of Λ are collected using the mechanical object to test, and the measurements of Λ 0 are collected using a reference object, where the absence of cracks has been checked by means of other more expensive techniques.

When the measurements are paired, we can define the operator Λ in the case of random background as in (4),

Λ(y) = Λ(y) -Λ 0 (y), y ∈ Γ, ( 29
)
where now the operator is parametrized by the realization y ∈ Γ of the random variable

Y . Therefore Λ : L2 (∂B) ⊗ L 2 (Γ, dρ) → L2 (∂B) ⊗ L 2 (Γ, dρ),
and Λ(y) : L2 (∂B) → L2 (∂B) is strongly compact for any value y ∈ Γ. At this point, the only difficulty that remains is that the value of y that generated the measurement is unknown. We will see in Section 3.2 how to cope with this issue. When the measurements are arbitrary, the operator Λ depends on y and y,

Λ(y, y) = Λ(y) -Λ 0 ( y), y, y ∈ Γ, ( 30
)
and the property of strongly compactness of the operator Λ(y, y) : L2 (∂B) → L2 (∂B) is lost in general. However, we know that the singular values of a strongly compact operator decay exponentially fast. So we can look at the singular values of Λ as defined in (30), and choose the value of y that yields the fastest decay of the spectrum of Λ.

The singular values of several strongly compact operators are computed in [START_REF] Haddar | Numerical analysis of the factorization method for electrical impedance tomography in inhomogeneous medium[END_REF]. We introduce the quantity K defined as

K(y, y) = supp x∈B |σ(x, y) -σ B (x, y)| = 0 , y, y ∈ Γ, (31)
that gives the size of the support of the contrast between σ and σ B , for any realization y and y in Γ. Of course

|D| = K(y, y) ≤ K(y, y) ≤ |B|, ∀ (y, y) ∈ Γ × Γ.
In general, the singular values of Λ increase with respect to the value of K. The next remark proves this statement in a particular situation where explicit calculations are made possible.

Remark 4 Suppose that B is a circular domain with unitary radius, and D is a concentric circular inclusion with radius ρ < 1. The diffusion coefficient is assumed to be homogeneous both in the background and in the inclusion. We define µ = σ B /σ D and observe that K defined in (31

) satisfies K = |supp x∈B {σ -σ B = 0}| = |D|. Of course we have ρ 2 ∝ K.
|λ DN i | = i µ + 1 -(µ -1)ρ 2i µ + 1 + (µ -1)ρ 2i , i = 1, 2, . . .
and therefore the eigenvalues of the dual Neumann-to-Dirichlet map decay as

|λ N D i | = 1 i µ + 1 + (µ -1)ρ 2i µ + 1 -(µ -1)ρ 2i , i = 1, 2, . . .
From this we can derive the decay of the eigenvalues of the operator Λ:

|λ i ( Λ)| = 1 i µ + 1 + (µ -1)ρ 2i µ + 1 -(µ -1)ρ 2i -1 = 1 i 2 ρ -2i (2(µ -1) -1 + 1) -1 , i = 1, 2, . . . ( 32
)
and clearly the smallest ρ < 1 the fastest the decay, as shown in Fig. 2. See also the numerical section where we tested inclusions with radius ρ = 0.1 and ρ = 0.0275. From (32), if we suppose µ > 1 we derive

|λ N D i | ∼ K i i µ -1 µ + 1 , i = 1, 2, . . .
The dependency on µ of the eigenvalues is very small, and it becomes unperceivable already from the third eigenvalue λ 3 on (see Fig. 2). Looking at the first eigenvalue we see a slight dependence of λ 1 on µ: the higher µ (that is equivalent to say the higher the contrast σ B -σ D ) the slower the decay. The case ρ = 0.1 and µ = 10000 has been tested in the numerical section, since we have set σ D = 10 -4 and σ B ≡ 1. Remark 4 and Remark 5 justify analytically that the decay of the eigenvalues is slower, the larger is K(y, y) with respect to ρ. Formally K(y, y) can attain the value |D|, but when the random coefficient is associated to a continuous probability density a contrast is always present with probability one. The same occurs also in the applications: there is always a contrast somewhere in the domain between the background coefficient of Λ and Λ 0 , because of modeling, measurement and numerical errors. As a consequence, the quantity K and the value of the contrast |σ(x, y) -σ B (x, y)| have both influence on the decay of the eigenvalues. Since the background coefficients σ and σ B are continuous w.r.t. the random variables (because of Assumption (2)), the magnitude of their contrast |σ(x, y) -σ B (x, y)| satisfies

sup x∈B |σ(x, y) -σ B (x, y)| ≤ C y -y ∞ , ∀ (y, y) ∈ Γ × Γ.
In conclusion, we proved that in a circular domain with a circular concentric inclusion and piece-wise constant background the decay of the eigenvalues of Λ is governed by yy ∞ . Reasonably the same conclusion can be extended to more general geometries and partitions of the domain into regions with arbitrary shape.

Consequently, when M = 1 and Λ(y) is the only measurement generated by the (unknown) realization y, we propose to solve the inverse problem

y =argmin y * ∈Γ Λ(y, y * ) k , (33)
to obtain y that satisfies

y -y ∞ ≤ ε. ( 34
)
The norm • k is defined as

• k = k i=1 σ i (•)
, and is known as Ky Fan k-norm. For any continuous linear operator in Hilbert spaces, the Ky Fan 1-norm corresponds to the usual operator norm equal to the largest singular value of the operator. The value of ε has to be small enough, depending on the size of the inclusion |D| and on the jumps of the background coefficient. At the discretization level, one can use also the Frobenius norm, that takes into account all the singular values including the smallest ones.

The problem ( 33) is an ill-posed inverse problem over the parameter space Γ, which dimension can be very large. If more than one measurement of the operator Λ are available, i.e. M > 1, it is possible to exploit this further information solving

y =argmin y * ∈Γ min 1≤m≤M Λ(y m , y * ) k , (35)
so that the objective function contains more points of minima displaced in Γ, and the optimization algorithm can detect more easily one of them which satisfies (34), being y one among the M realizations that generated the M measurements. We will discuss in the numerical section how to efficiently solve problems (33) and (35).

The Factorization Method in the case of arbitrary measurements

In the case of random background with arbitrary measurements, we propose two versions of the FM. The first version of the FM that we present is an optimized variant (Algorithm 2) where the background is reconstructed by means of an optimization approach. This optimization is mandatory whenever only one measurement or few measurements of the operator Λ are available. The additional operations w.r.t. Algorithm (1) concern the approximation of the solution y of the optimization problem (33) when only one measurement of Λ is available, or (35) when more than one measurement are available. The second version is named pure variant of the FM (Algorithm 3), because the final imaging reconstruction is a weighted linear combination of intermediate imaging reconstructions, each one obtained without any attempt to reconstruct the background. The approximation of the solution of problem (35) satisfying the convergence criterion (34) within a tolerance ε arbitrarily small is found by means of a sampling approach, rather than by optimization. This variant becomes attractive only when many measurements are available.

Algorithm 2

The optimized Factorization Method in random background (arbitrary measurements)

Sample the region of B to be probed with a set of points P = {ξ j } P j=1 , Choose the initial guess y (i) ∈ Γ, i = 0, repeat Compute y (i+1) performing one or several steps in the optimization procedure for solving problem (33) with initial guess y (i) . for ξ in the set P do solve problem [START_REF] Guo | A review of electrical impedance techniques for breast cancer detection[END_REF] with

σ B = σ B (•, y (i) ) to find its solution N (•, ξ, y (i) ), compute l k ξ ( y (i) ) from N (•, ξ, y (i) ) using (7), use l k ξ ( y (i) ) to compute the indicator C(ξ) in (22), plot ξ → C(ξ), end for i ← i + 1 until the background is NOT reconstructed in C(ξ) OR i ≥ maxit .
To set up an easy notation, we consider only the parametrization (26) of the coefficient σ B with the vector d. The same idea can be applied with any parametrization of σ B . Define the l ∞ norm of an N -dimensional vector v, weighted w.r.t. the Ndimensional vector d as v ∞,d = max{d 1 v 1 , . . . , d N v N }. Given the values of N , M , ε and the vector d, then the number T of realizations in the sampling procedure is chosen such that

P min 1≤m≤M 1≤t≤T y m -y t ∞,d ≤ ε ≥ 0.99. ( 36
)
Note that ε in (36) plays the same role as in (34). If the dispersion vector d has components all equal to d > 0, so that the dispersion is the same for all the components of the random variable Y, then (36) can be rewritten with the usual l ∞ norm as

P min 1≤m≤M 1≤t≤T y m -y t ∞ ≤ ε d -1 ≥ 0.99. ( 37
)
In the pure FM T is choosen such that (36) is satisfied, and this allows to obtain with 99% probability an approximation of one of the realizations that generated the measurements, within a tolerance ε. In this way the optimization of problem (35) can be avoided. Note that, given M , N , d and ε, condition (37) states how to choose T such that the pure FM detects the inclusions with at least 99% probability. In Section 6.2 the calibration of ε is discussed and the value of T such that condition (37) holds is calculated for some values of M, N, d, ε.

Looking at Algorithms 2 and 3, it is clear that the most computationally expensive part of these two variants of the FM consists in the repetitive evaluation of Λ 0 , for many realizations of the random variable Y. In Algorithm 2 this operation is hidden inside ) to find its solution N (•, ξ, y (i) ), compute l k ξ ( y (i) ) from N (•, ξ, y (i) ) using ( 7), use l k ξ ( y (i) ) and the SVD decomposition of Λ(y

m i , y i ) to compute g k ξ (y m i , y i ) L 2 (∂B) , use g k ξ (y m i , y i ) L 2 (∂B) to compute the indicator C i (ξ) in (22), update C(ξ) ← C(ξ) + C i (ξ) exp - (min 1≤j≤i F m i j ) F m i i , plot ξ → C(ξ), end for i ← i + 1 until the background is NOT reconstructed in C(ξ) OR i ≥ T .
the optimization step, where point-wise evaluations of Λ 0 are needed. The number of evaluations of Λ 0 to achieve an accurate approximation of the solution of problems (33) and (35) can be considerably reduced by means of polynomial approximation [START_REF] Babuška | Galerkin finite element approximations of stochastic elliptic partial differential equations[END_REF][START_REF] Cohen | Analytic regularity and polynomial approximation of parametric and stochastic elliptic PDE's[END_REF] of problem (2) with the random coefficient (25), which defines the operator Λ 0 .

The Factorization Method in the case of paired measurements

When a pair of measurements of Λ and Λ 0 is available, i.e. the measurements are generated from the same realization y of the random variable Y , then we can build straightforwardly the operator Λ from (29). Since the value of the realization y that generated the measurements is unknown, we have to work with the sampling procedure to span the parameter space Γ looking for a good approximation of y.

Let {s t } T t=1 be T independent and identically distributed realizations of the random variable Y that are used to explore the parameter space Γ. We propose the following indicator function, that substitutes the one defined in [START_REF] Haddar | Conformal mapping and impedance tomography[END_REF]:

C(ξ) = T t=1 1 log 10 v 1 (ξ, s t ) + v 2 (ξ, s t ) , (38)
In (38) the functions v k : (ξ, y)

∈ (D × Γ) → v k (ξ, y) : D × Γ → R are defined as v k (ξ, y) := g k ξ (y) L 2 (∂B)
and g k ξ (y) is again given by ( 19) but with l p ξ that depends on y, as in (28). The indicator (38) is based on the sample mean estimator, which is very sensitive to the outliers that correspond to the points where it diverges due to the presence of the inclusion. This variant of the FM that exploits the indicator function (38) is described in Algorithm (4). The number T of realizations to be evaluated with the FM has to be tuned according to the dimension N of the parameter space Γ.

Algorithm 4

The Factorization Method in random background (paired measurements) Sample the region of B to be probed with a set of points P = {ξ j } P j=1 , Sample T realizations s 1 , . . . , s T iid ∼ ρ of the random variable Y, for ξ in the set P do for t = 1, . . . T do solve problem [START_REF] Guo | A review of electrical impedance techniques for breast cancer detection[END_REF] with σ B = σ B (•, s t ) to find its solution N (•, ξ, s t ), compute l k ξ (s t ) from N (•, ξ, s t ) using ( 7), update the indicator C(ξ) in (38) with l k ξ (s t ), end for plot ξ → C(ξ). end for

Numerical tests featuring deterministic background

In this section we present several numerical tests to illustrate the capabilities of the Factorization Method in homogeneous and inhomogeneous but deterministic background. To show the results, we display the isolines of the indicator function C(ξ) given by ( 22) A crucial issue is the tuning of the scale and the choice of the isovalue that represents the inclusion. In practice, this requires additional information on the value of the coefficient in the inclusions. In all the numerical tests we plot the indicator C(ξ) in the domain B choosing an uniform coloring scale, such that the black color is associated to the zero value and the white color is associated to max ξ∈B C(ξ). Throughout the paper the space L2 (∂B) is discretized using the Fourier basis with 128 terms.

The domain B is always the unitary ball centered in (0, 0). The sampling points consist of a 50 × 50 uniform grid over the square [-1, 1] × [-1, 1]. In the FM we use only the points that fall 0.05 farther from ∂B. In the testcases with small inclusions the resolution is increased to 100 × 100. The dashed purple line always marks the exact geometry of the inclusion(s). The value of the coefficient that identifies the inclusion is always set to σ D = 0.001, although any value that satisfies Assumption 1 is allowed.

We focus on two configurations of the background σ B :

• (E) the union of two semicircular regions,

• (U) the union of five anular regions.

The inclusion(s) can be configured as

• (N) a circular inclusion with radius 0.1 centered in (0, 0.33),

• (S) four small circular inclusions with radius 0.0275.

Combining each configuration of the background with each configuration of the inclusions we obtain the four configurations EN, ES, UN, US depicted in Fig. 3.

-

R 2 R 6 R 1 R 4 R 5 R 3

The homogeneous case

We begin to apply the Factorization Method to the EIT problem with an homogeneous deterministic coefficient σ B ≡ 1. Many numerical results in the literature show the potentiality of the method. The testcases are named as

• testcase hN: inclusion configuration N with m i = 1, i = 1, . . . , r, in (24),

• testcase hS: inclusion configuration S with m i = 1, i = 1, . . . , r, in [START_REF] Lagarias | Convergence properties of the Nelder-Mead simplex method in low dimensions[END_REF].

With this choice of m the background configurations E and U become identical. Fig. 4 shows the results obtained with the inclusion configurations N and S, displaying the coloring of the indicator C(ξ) in the domain B. In both cases, the reconstructions are very accurate and the shape and location of the inclusions are clearly detected.

The inhomogeneous case

We now proceed to present the results concerning the inhomogeneous deterministic coefficient σ B . In this case the Factorization Method is described in Algorithm 1. As Numerical analysis of the Factorization Method for EIT with uncertain background 18 will be explained in the next section, we pick the background coefficient equal to one in the region R 1 that is accessible to the boundary. In Fig. 5 and Fig. 6 we present the results obtained when applying the Factorization Method in the following testcases,

• testcase i1EN: configuration EN with m = [1, 5, 0.5] in ( 24),

• testcase i1ES: configuration ES with m = [1, 5, 0.5] in [START_REF] Lagarias | Convergence properties of the Nelder-Mead simplex method in low dimensions[END_REF], to achieve one order of magnitude jumps between different regions, or chosen as

• testcase i2EN: configuration EN with m = [1, 10, 0.1] in [START_REF] Lagarias | Convergence properties of the Nelder-Mead simplex method in low dimensions[END_REF],

• testcase i2ES: configuration ES with m = [1, 0.1, 10] in [START_REF] Lagarias | Convergence properties of the Nelder-Mead simplex method in low dimensions[END_REF], to achieve two orders of magnitude jumps. Next, in Fig. 7 and Fig. 8 we present the results of the following testcases,

• testcase i1UN: configuration UN with m = [1, 5, 1, 5, 0.5, 0.5] in [START_REF] Lagarias | Convergence properties of the Nelder-Mead simplex method in low dimensions[END_REF],

• testcase i1US: configuration US with m = [1, 5, 1, 0.5, 1, 0.5] in [START_REF] Lagarias | Convergence properties of the Nelder-Mead simplex method in low dimensions[END_REF], to achieve one order jumps between different regions, or set as order of magnitude in Fig. 5 and Fig. 7 is the same as the quality in Fig. 4 with the homogeneous background. See also many other numerical tests in [START_REF] Haddar | Numerical analysis of the factorization method for electrical impedance tomography in inhomogeneous medium[END_REF], also with more general nonlinear inhomogeneous background coefficients. When the order of magnitude of the jumps increases to two, the effect of the jump in the background coefficient becomes visible in the reconstruction, as in Fig. 8-left. Nonetheless the inclusion is still accurately detected. In Fig. 6-right we observe that the inclusions on the interface where the coefficient jumps are still detected, albeit with a fainter intensity of the indicator C(ξ). The same occurs in Fig. 8-right with the most inner inclusion. All the four small inclusions are detected, but the method would hardly inspect the presence of additional inclusions inside their convex envelope. See also [START_REF] Haddar | Numerical analysis of the factorization method for electrical impedance tomography in inhomogeneous medium[END_REF], where a test with ten small inclusions is performed.

•
We can conclude that the Factorization Method is able to accurately detect the shape and location of all the inclusions, when the jumps in the background are up to two orders of magnitude.

Numerical tests with one measurement featuring random background

In this section we present some numerical tests with random background, in the arbitrary measurements and paired measurements cases. We are interested in large variations of the background coefficient σ B , and therefore we parametrize it with the exponential model ( 27) so that the coefficient jumps by two orders of magnitude in each one of the r regions, or

σ B (x, Y) = r i=1 0.5 • 10 Y i I R i (x), x ∈ B, Y ∈ Γ = [-0.5, 0.5] r , (39)
to have one order jumps. Usually in the problem of EIT the boundary is accessible. This would allow to directly recover the value of the realization of the random variable Y 1 that generated the voltage and current, corresponding to the observation of Λ. As a consequence, the value of σ B in R 1 could be retrieved, and the starting problem could be reduced to a problem with deterministic coefficient in R 1 . For this reason, we will choose σ B equal to 1 in R 1 and random only in Ω \ R 1 . In this way, the parameter space Γ in the parametrizations (26), (39), (27) has to be dimensionally reduced by one, so that its dimension in practice becomes N = r -1.

We proceed now to present several numerical tests performed with the optimized FM (Algorithm 2) in the case of one arbitrary measurement, and with the FM of Algorithm (4) in the case of one paired measurement.

The optimized Factorization Method with one arbitrary measurement.

To show the capabilities of the proposed approach we set up three testcases, with an increasing difficulty in the optimization operations. We choose the following testcases featuring a random background:

• testcase aEN2: configuration EN with two orders jumps in the background coefficient defined in (27). The parameter space is Γ = [0.1, 10] 2 . The measurement of Λ(y) is generated by the realization y = (8.3567, 0.3558).

• testcase aES2: configuration ES with two orders jumps in the background coefficient defined in (27). The parameter space is Γ = [0.1, 10] 2 . The measurement of Λ(y) is generated by the realization y = (4.5, 0.75).

• testcase aUN1: configuration UN with one order jump in the background coefficient defined in (39). The parameter space is Γ = [0.5, 5] 5 . The measurement of Λ(y) is generated by the realization y = (4.4950, 1.3450, 0.9450, 3.3450, 0.6240).

The testcase aEN2 is a standard configuration, with a not too small inclusion and a twodimensional parameter space Γ. The testcase aES2 is more challenging than testcase aAS2, because the inclusions have a smaller diameter and their presence is more easily obfuscated by the variations of the random background. Lastly, in testcase aUN1 we treat an ill-posed optimization problem over a five-dimensional parameter space, and the limits of the optimization approach are highlighted. From a computational standpoint, the optimization of the Ki Fan k-norm • k can be equivalently replaced by the optimization of the real-valued operator q k (•) = k i=1 log σ i (•), which has the same points of minima as • k . In Figs. 9, 10 and 11 we display the function q k evaluated over the parameter space Γ in the testcases aEN2 and aES2, for some values of k. As explained in section 3, the solution of problem (33) that corresponds to the point of global minima is the realization y that generated the measurement of Λ(y). The higher the value of k, and the steeper the objective function becomes, because more singular values are taken into account. Moreover, when k increases the objective function shows a product structure that suggests the use of Alternating Minimization. When Γ is two-dimensional, we present the numerical results only using a global optimization approach, although we observed that the Alternating Minimization allows to reduce the overall computational cost.

When Γ = [0.5, 5] 5 , the global optimization approach converges to the correct solution in the subspace of the parameter space associated to the three outer regions, but then gets stuck into local minima when exploring the two most inner regions. Alternating Minimization converges even more easily, but still gets stuck into local minima whenever the initial data of the two most inner regions is not close enough to the exact solution. One possibility is to initialize with different initial data the optimization procedure, but since the parameter space is very large this turns out to be very costly. Therefore, we think that the optimization procedure could still be employed effectively but on problem (35) with more than one measurement of Λ and resizing the range of variation of the random variables of the most inner regions. For example, when Y 5 , Y 6 deviate only 25% from their mean values the alternating minimization always easily finds one of the points of minima of problem (35). Some test with many measurements are presented in Section 6. To solve problem (33) we choose a derivative-free approach, since ∂ y * Λ(y, y * ) k is not available. We employ the well-known Melder-Mead method [START_REF] Lagarias | Convergence properties of the Nelder-Mead simplex method in low dimensions[END_REF] to perform a constrained optimization of the objective function q k over the whole parameter space Γ.

In testcase aEN2 the method converged in less than 100 iterations to the solution y of problem (33) within an accuracy ε = 10 -4 in (34), for several different choices of the initial point in the parameter space. The number of evaluations of Λ 0 has always been less than three times the number of iterations. In testcase aES2 the method converged in less than 100 iterations to the solution y of problem (33) with an accuracy ε = 10 -3 , for several different choices of the initial point in the parameter space. The evaluations of Λ 0 are more costly, because the finite element mesh has been chosen finer.

-1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1
Fig. 12. Coloring of C(ξ), aEN2, γ = 10 -2 . Top-left: using the solution y of (33) with ε = 0.00014. Top-right: using the solution y of (33) with ε = 0.001. Bottom-left: using the solution y of (33) with ε = 0.01. Bottom-right: using the point y = (5.05, 5.05).

Once we obtain an approximation y of the solution of problem (33), we can apply the FM with the realization y. In Fig. 12 we show the reconstructions obtained in the testcase aEN2, with different values of ε. We solve problem (33) with the Ky Fan 5-norm with initial point (0.1, 10), and after 76 iterations and 223 evaluations of Λ 0 the method converges to the solution y = (8.356653, 0.355660) such that the convergence criterion (34) is satisfied with ε = 0.00014. The corresponding reconstruction obtained is depicted in Fig. 12-top-left. The other Figs. 12-top-right and 12-bottom-left show the effect of using a less accurate realization y for different values of ε. Fig. 12-bottom-right is obtained using y = (5.05, 5.05), that is the expected value of the random variable Y over Γ, and yields a reconstruction where the inclusion is completely obfuscated by the contrast of the background. As justified in section 3, when applying the FM with a realization y that is not enough accurate with respect to the criterion (34), the background begins to be detected in the reconstruction. This fact can be used to check a posteriori that the solution of problem (33) is accurate enough: if the region R i of the background is sharply detected, then the ith component of the solution y of problem (33) is too inaccurate, and its accuracy needs to be improved. To achieve an improvement one can either utilize Alternating Minimization, or increase k or choose other types of optimization methods, e.g. Simulated Annealing or Global Pattern Search.

In Fig 13-left we show the reconstruction obtained in testcase aES2 using the solution y of problem (33) that satisfies (34) with ε = 10 -3 . This solution was computed solving problem (33) with the Ky Fan 32-norm with initial point (0.1, 10), and required 131 iterations and 312 evaluations of Λ 0 .

Finally we present the testcase aUN1, that shows the limit of the optimization procedure with large variations of the coefficients. In Fig. 14 we show the results obtained for different values of ε: on top-left we show the reconstruction obtained when solving accurately problem (33) with ε = 5 • 10 -4 , on bottom-right we display the reconstruction obtained when the optimization of problem (33) is accurate with ε = 10 -2 only in the three most outer regions, but converged only up to ε = 0.2 in the subspace of Γ associated to the two most inner regions. As a result, the two most inner regions are detected as an inclusion.

Numerical analysis of the Factorization Method for EIT with uncertain background 25

-1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1
Fig. 14. Coloring of C(ξ), aUN1, γ = 10 -2 . Top-left: using the solution y of (33) with ε = 0.0005. Top-right: using the solution y of (33) with ε = 0.001. Bottom-left: using the solution y of (33) with ε = 0.01. Bottom-right: using the solution y of (33) with ε = 0.01 in R 2 , R 3 , R 4 and ε = 0.2 in R 5 , R 6 .

In general, we observed that the global optimization over the whole five-dimensional parameter space is always successfully, when applied to testcases featuring small variations of the coefficients in the most inner regions, e.g. up to 25%. On the other side, when the starting point is not sufficiently close to the exact solution and the variations of the coefficients are of one order of magnitude in the two most inner regions, the global optimization over the five-dimensional parameter space got stuck into local minima. Clearly, this is a limitation of the proposed approach, and the availability of more than one measurement to facilitate the optimization seems required to overcome it.

The Factorization Method with one paired measurement

In the case of paired measurements, with the configurations EN and UN the inclusion is always easily detected, also when the background coefficient jumps by two orders of magnitude. Therefore, we will not present these results, but focus only on the more challenging configurations ES and US, again with two orders of magnitude jumps in the background coefficient. We consider the following testcases, to test the FM in the paired measurements case:

• testcase pES: configuration pES with σ B defined in (27),

• testcase pUS: configuration pUS with σ B defined in (27).

Numerical analysis of the Factorization Method for EIT with uncertain background 26

The results obtained are presented in Fig. 15 for the testcase pES, and in Fig. 16 for the testcase pUS. We use the indicator C(ξ) in (38) with the number T of realizations chosen according to the number of regions r. We choose T = 10 when r = 3 and T = 20 when r = 6. Then we apply the FM when observing only one measurement of the random operator Λ, generated from the unknown realization y. The results are presented when "common" realizations are observed, e.g. σ B = (1, 5, 0.5) or σ B = (1, 0.26, 4.17, 0.89, 7.05, 0.16), and when "extremal" realizations are observed as well. We distinguish between "common" and "extremal" realizations in the following sense: an extremal realization lies on the boundary of the support of the random variable, e.g. σ B = (1, 10, 0.1) and σ B = (1, 10, 0.1, 0.1, 10, 10). If a realization is not extremal, than it is a "common" realization. Extremal realizations provide a worst-case benchmark for the FM. With all the realizations observed, the FM detects at least three inclusions over four. The inclusion closest to the center of the domain is the hardest to detect. The reconstructions rely on the unique realization of the random operator Λ, and this justifies its slight dependence on the particular realization observed. In spite of this, at least three inclusions are always sharply detected, but one is occasionally undetected. 6. Numerical tests with many measurement featuring random background 6.1. The optimized Factorization Method with many arbitrary measurements.

-1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0
As described previously, the optimization in problem (33) can be facilitated, optimizing problem (35) where additional measurements of Λ(Y ) are included. Fig. 17 shows the evaluation over Γ of the objective function of problem (35) with • k replaced by q k (•), in the testcases aEN2 and when four additional measurements are added to the one measurement of Λ available.

6.2. The pure Factorization Method with many arbitrary measurements.

Here we discuss some details about the application of the pure Factorization Method. The critical issue concerns the choice of ε. In Figs. 12-top-right, 13-left, 14-top-right several accurated reconstructions are obtained with ε = 10 -3 and much larger variations in the coefficient. Based on this observation, and on many other numerical experiments that we do not report, an accuracy ε = 10 -3 is always enough to clearly detect inclusions of normal size, e.g. with radius 0.1, when there are up to r = 4 regions, d is not larger than 0.125 and the vector m does not produce jumps in the coefficient larger than one order of magnitude. When interested in smaller inclusions, with more regions r > 4, or d and m that produce larger variations in σ B , then also the value of ε has to be furtherly decreased.

Table 1 reports the number T of realizations needed to satisfy (36), for some values of M , N , d, ε = 10 -3 , σ B parametrized by (26) and d up to 0.125, yielding 25% variations in the coefficient. Note that this variation is much smaller than those treated in Figs. 12-top-right, 13-left, 14-top-right.

In Table 2 we report the value of T such that a condition analogous to (36) holds, but with the exponential parametrizations (39) and (27). The ranges of variations achieved with these parametrizations are the same as in Figs. 12-top-right, 13-left, 14-top-right. We remark that the value of T directly corresponds to the number of evaluations of the operator Λ 0 in points of the parameter space Γ. This evaluation can be performed very efficiently by means of polynomial approximation techniques.

Conclusions

In this paper we provided a numerical scheme to solve the Neumann dipole-like boundary value problem in inhomogeneous deterministic background. Then we employed this scheme to apply the FM to the problem of EIT in inhomogeneous deterministic and random background. With any type of background the FM is able to detect the presence and location of one or many inclusions, and to supply information on their geometry. The quality of the reconstructions obtained in the inhomogeneous case is comparable to the quality of those in the homogeneous case, with jumps in the background up to two orders of magnitude.

When the background is random and the jumps are of two orders of magnitude, the inclusions are reconstructed at a lower accuracy. We provided three variants of the FM, based on the distinction between arbitrary and paired measurements of the random operator Λ. In the paired measurement case the FM always detects a normal or small inclusion. In presence of many inclusions those harder to detect could remain undetected. In the case of arbitrary measurements, the quality of the reconstruction depends on the accuracy ε of the criterion (34) when solving problem (33). We focused mainly on large variations of the coefficients, and showed that when the parameter space is low-dimensional the solution of problem (33) can be easily approximated up to a tolerance ε such that the FM is successful. When the dimension of the parameter space increases the optimization becomes challenging, in particular with large variations of the background coefficient. Another potentiality of the proposed approach, is the possibility to easily incorporate many measurements of the operator Λ in problem (35), which optimization can be a lot quicker than for problem (33) with only one measurement. Lastly we proposed a pure variant of the FM, that is particularly suited for uncertainty analysis and does not require the optimization step.

Fig. 1 .

 1 Fig. 1. An obstacle D which lies in a medium B.

Fig. 2 .Remark 5

 25 Fig. 2. The decay of the eigenvalues λ i ( Λ) in (32) for different choices of ρ and µ.

Algorithm 3

 3 The pure Factorization Method in random background (arbitrary measurements) Collect M measurements Λ(y 1 ), . . . , Λ(y M ), Sample the region of B to be probed with a set of points P = {ξ j } P j=1 , Sample T realizations y 1 , . . . , y T iid ∼ ρ of the random variable Y, i = 1, repeat Compute F mi = Λ(y m ) -Λ 0 ( y i ) k , ∀ m and m i = argmin m=1,...,M F mi , for ξ in the set P do solve problem (5) with σ B = σ B (•, y (i)

Fig. 3 .

 3 Fig. 3. Geometries of the regions in the background and displacement of the inclusions (in dashed line). Top-left: configuration EN. Top-right: configuration ES. Bottom-left: configuration UN. Bottom-right: configuration US.

Fig. 4 .

 4 Fig. 4. Coloring of C(ξ) . Left: hN, γ = 10 -2 . Right: hS, γ = 5 • 10 -3 .

Fig. 5 .

 5 Fig. 5. Coloring of C(ξ). Left: i1EN, γ = 10 -2 . Right: i1ES, γ = 7 • 10 -3 .

Fig. 6 .Fig. 7 .Fig. 8 .

 678 Fig. 6. Coloring of C(ξ). Left: i2EN, γ = 10 -2 . Right: i2ES, γ = 10 -2 .

Fig. 9 .

 9 Fig. 9. Testcase aEN2. Function q k over the parameter space Γ in log-log scale. Top-left: k = 1. Top-right: k = 5. Bottom-left: k = 32. Bottom-right: k = 64. The realization y = (8.3567, 0.3558) corresponds to the point of global minima (0.9220, -0.4488) in log-log scale.

Fig. 10 .

 10 Fig. 10. Testcase aES2. Function q k over the parameter space Γ in log-log scale. Left: k = 32. Right: k = 64. The realization y = (4.5, 0.75) corresponds to the point of global minima (0.6532, -0.1249) in log-log scale.

Fig. 13 .

 13 Fig.[START_REF] Harrach | Detecting inclusions in electrical impedance tomography without reference measurements[END_REF]. Coloring of C(ξ), aES2, γ = 10 -2 . Left: using the solution y of (33) with ε = 10 -3 . Right: using the solution y of (33) with ε = 10 -2 .

1 - 1 0 1 Fig. 15 .

 11115 Fig. 15. Coloring of C(ξ) with T = 10, pES, γ = 5 • 10 -2 . Top-left: σ B = (1, 0.5, 5). Top-right: σ B = (1, 5, 0.5). Bottom-left: σ B = (1, 0.1, 10). Bottom-right: σ B = (1, 10, 0.1).

Fig. 16 .

 16 Fig. 16. Coloring of C(ξ) with T = 20, pUS, γ = 5 • 10 -2 . Top-left: σ B = (1, 0.26, 4.17, 0.89, 7.05, 0.16). Top-right: σ B = (1, 10, 0.1, 0.1, 10, 10). Bottom-left: σ B = (1, 10, 1, 0.1, 1, 0.1). Bottom-right: σ B = (1, 0.1, 0.1, 10, 0.1, 10).

Fig. 17 .

 17 Fig. 17. Testcase aEN2 with 5 measurements of Λ(Y ). Evaluation over Γ in loglog scale of the objective function of problem (35), with • k replaced by q k . The realizations y 1 , y 2 , y 3 , y 4 , y 5 correspond to the five points of minima. Left: k = 1. Right: k = 32.

dTable 1 .

 1 = 0.025 d = 0.05 d = 0.1 d = 0.125 M = 3, N = 2 (r = 3Uncertainty Analysis. The value of T such that P[min mt y my t ∞ ≤ ε d -1 ] ≥ 0.99, given ε = 10 -3 and some values of N , M , d. Parametrization of σ B by (26), where y m , y t ∈ [-1, 1] N . The dash advises for a larger value of M .

Table 2 .

 2 Large variations. The value of T such that P[min mt y my t ∞ ≤ ε] ≥ 0.99, given ε = 10 -3 and some values of N , M . Left column: one order of magnitude variations with the parametrization (39) of σ B , where y m , y t ∈ [0.5, 5] N . Right column: two orders of magnitude variations with the parametrization (27) of σ B , where y m , y t ∈ [0.1, 10] N .

  Testcase aES2. Same as Fig.10but zoomed in a neighbourhood of the point of global minima (0.6532, -0.1249). Function q k over the parameter space Γ in log-log scale. Left: k = 32. Right: k = 64.

										-82
	-0.12										-0.12									-212
	-0.121									-84	-0.121									-214
	-0.122									-86	-0.122									-216
	-0.123										-0.123
										-88										-218
	-0.124										-0.124
										-90										-220
	-0.125										-0.125
																				-222
	-0.126									-92	-0.126
	-0.127										-0.127									-224
										-94
	-0.128										-0.128									-226
										-96
	-0.129										-0.129									-228
	-0.13									-98	-0.13
																				-230
	0.644	0.646	0.648	0.65	0.652	0.654	0.656	0.658	0.66	0.662	0.644	0.646	0.648	0.65	0.652	0.654	0.656	0.658	0.66	0.662
	Fig. 11.


Jordi Sanchez-Riera
email: jordi.sanchez-riera@inria.fr

Xavier Alameda-Pineda
email: xavier.alameda-pineda@inria.fr

Radu Horaud
email: radu.horaud@inria.fr

Audio-Visual Robot Command Recognition D-META'12 Grand Challenge

Keywords: H.4 [Information Systems]: Applications Algorithms, Datasets, Evaluation Audio-Visual Categorization, Multimodal Learning

This paper addresses the problem of audio-visual command recognition in the framework of the D-META Grand Challenge 1 . Temporal and non-temporal learning models are trained on visual and auditory descriptors. In order to set a proper baseline, the methods are tested on the "Robot Gestures" scenario of the publicly available RAVEL data set, following the leave-one-out cross-validation strategy. The classification-level audio-visual fusion strategy allows for compensating the errors of the unimodal (audio or vision) classifiers. The obtained results (an average audio-visual recognition rate of almost 80%) encourage us to investigate on how to further develop and improve the methodology described in this paper.

INTRODUCTION

Humans use several sensory modalities to gather information into order to build an internal representation of the real world, e.g, visual, auditory, tactile, or olfactory sensing. Among these modalities, vision and audition are the most suitable to be used by a robot due to the wide availability of associated sensors, namely cameras and microphones. Nevertheless, the data collected by these sensors are often corrupted by occlusions, reverberations, and bad recording conditions such as poor lighting, competing sound sources, background noise, etc. Hence unimodal interpretation leads to ambiguities. It seems natural to investigate how one can obtain a better understanding of a human action if several modalities are used simultaneously. Indeed, the fusion of information coming from different modalities is likely reinforce the true hypothesis. In addition, several challenging questions arise when dealing with different modalities: At which level of abstraction one should fuse the information coming from physical different sensors? How can one build the concept of multimodal features? how can one fuse features representing different aspects (time, position, appearance, spectrum, etc) of the events to be recognized or categorized?

This work is presented in the framework of the D-META ICMI Grand Challenge aiming at setting up the basis for comparison, analysis, and further improvement of multimodal data annotations and multimodal interactive systems. The main goal of D-META is to foster research and development in multimodal communication and to further elaborate algorithms and techniques for building various multimodal applications. Held by two coupled pillars, method benchmarking and annotation evaluation, the D-META challenge envisions a starting point for transparent and publicly available application and annotation evaluation on multimodal data sets. Five tasks were proposed in D-META: AVRGR, recognize gestures addressed to the robot by means of vision and the audio, AVSR, detect, localize and track multiple speakers using audio-visual information, CEP, estimate the level of engagement in a video-mediated communication, AVCGR, recognize conversational gestures in first encounter dialogues and AVFGR, recognize feedback gestures in first encounter dialogues.

In this paper we set the baseline for the task audio-visual recognition of human commands. The proposed method combines auditory and visual information for recognition using the publicly available RAVEL data set [START_REF] Alameda-Pineda | Ravel: An annotated corpus for training robots with audio visual abilities[END_REF], which contains annotated binocular-binaural recordings.

The remainder of the paper is organized as follows. Section 2 describes the previous work done on the topic and the main contributions of the paper. The learning approach as well as the auditory/visual features are described in section 3. In section 4 are explained the different experiments carried out to validate the fusion methods and in section 5 conclusions and directions for future work are briefly presented. Finally, section 6 is devoted to the discussion of some Challenge organizational aspects.

RELATED WORK & CONTRIBUTIONS

As mentioned before, a few choices have to be done when combining different sensory modalities for classification. Indeed, the level of abstraction at which the fusion should be done depends on the targeted application. Most of the existing audio-visual categorization methods perform the fusion either at the feature level [START_REF] Jiang | Short-term audio-visual atoms for generic video concept classification[END_REF][START_REF] Liu | An audio-visual fusion framework with joint dimensionality reduction[END_REF][START_REF] Lopes | Audio and video feature fusion for activity recognition in unconstrained videos[END_REF] or at the classification level [START_REF] Lacheze | Audio/video fusion for objects recognition[END_REF][START_REF] Saenko | Object category recognition using probabilistic fusion of speech and image classifiers[END_REF][START_REF] Xiong | Audio-visual sports highlights extraction using coupled hidden markov models[END_REF]. The reader is referred to [START_REF] Luo | Object category detection using audio-visual cues[END_REF] for a study using support vector machines (SVM) that compares feature-level fusion techniques to classification-level fusion techniques.

An example of a feature-level fusion method is [START_REF] Jiang | Short-term audio-visual atoms for generic video concept classification[END_REF]. The idea is to track short-term visual features. These features are associated to the auditory signal in order to construct a joint audio-visual codebook. Used to describe a part of a sequence, the codebook feeds a multiple-instance learning (MIL) paradigm, training a semantic concept detector. Another way to perform audiovisual fusion at a feature level is to map both modalities into the same "audio-visual object" space and to perform clustering [START_REF] Vasil Khalidov | Conjugate mixture models for clustering multimodal data[END_REF][START_REF] Alameda-Pineda | Finding audio-visual events in informal social gatherings[END_REF].

The authors in [START_REF] Liu | An audio-visual fusion framework with joint dimensionality reduction[END_REF] target a speech detection application, and perform the audio-visual fusion at a feature level as well. Principal component analysis (PCA) features are taken from the face images and Mel Frequency Cepstral Coefficients (MFCC) are the auditory features. Both types of features are then projected in a joint subspace using canonical correlation analysis (CCA). A Gaussian mixture model (GMM) is used for classification.

The authors in [START_REF] Lopes | Audio and video feature fusion for activity recognition in unconstrained videos[END_REF] target general activity recognition. High dimensional features (around 3000) for video and audio are collected and then reduced using the sequential forward floating selection (SFFS) algorithm. The SFFS is a pruning algorithm that selects most relevant features. With this the dimension is reduced (to 40). Finally a kNN algorithm is used as classifier.

Among the classification-level fusion methods we remark [START_REF] Lacheze | Audio/video fusion for objects recognition[END_REF], in which the authors experiment different combination strategies for object detection. Visual features are based on texture description and entropy-based variable-size patches. Auditory features correspond to the energy of the signal's gammatone filter bank decomposition. Monocular video and monaural audio are used and there is a strong need of uniform visual background.

Object recognition based SVMs is used in [START_REF] Saenko | Object category recognition using probabilistic fusion of speech and image classifiers[END_REF] where a probabilistic method combining posterior class probability output by each classifier is proposed; Basically this means that each modality is trained separately and then combined. SIFT descriptors are used as visual features and a commercial speech recognizer is used to classify the incoming audio signal.

A different approach from the ones mentioned so far is described [START_REF] Xiong | Audio-visual sports highlights extraction using coupled hidden markov models[END_REF] which finds sport highlights using a coupled hidden Markov model (CHMM). Several video features are used such as quantization average motion vectors and color. On the auditory side, the authors chose to use MFCC features. Both these features train a CHMM to perform the classification.

In this paper we introduce a method that compares different learning schemes. The first scheme is based on a bagof-words (BoW) approach while the second one incorporates temporal structure by means of HMMs. Scene-flow [START_REF] Sanchez-Riera | Action recognition robust to background clutter by using stereo vision[END_REF] and STIPs [START_REF] Laptev | On space-time interest points[END_REF] are used as visual features and MFCC as auditory features. The fusion is done at the classification stage through a modality-weighting scheme, i.e., pooling. Experimental validation is done using a publicly available data set, in which we set the performance baseline.

AUDIO-VISUAL CATEGORIZATION

This section is devoted to the proposed audio-visual command learning approach, which performs classification-level fusion. By means of the scene flow from one side and STIPs on the other side, we are able to describe the visual information (see section 3.1). The auditory information is characterized by standard features used in speech recognition (section 3.2). The learning is performed through two different supervised techniques: without any temporal model (section 3.3) and with a temporal model (section 3.4). Finally, the procedure to combine the output of the unimodal classifiers is described (section 3.5).

The Visual Descriptor

We used two different visual descriptors. The first one was proposed in [START_REF] Sanchez-Riera | Action recognition robust to background clutter by using stereo vision[END_REF] and it is based on the scene flow, which is the 3D equivalent of the optical flow [START_REF] Čech | Scene flow estimation by growing correspondence seeds[END_REF]. The scene flow is represented by the optical flow plus the depth at each image position. Together with the camera calibration, this is equivalent to a vector field of 3D position and associated 3D velocities. This intrinsic representation is potentially less sensitive to changes of texture and illumination than the intensity images. Moreover, the notion of depth allows to focus on the actor, while discarding any activity from the background. We assume that the actor of interest is the person closest to the camera. This is a reasonable assumption, since it holds in most of the human-robot-interaction applications and on movies. The final descriptor consists on the position and disparity relatives to the actor's face plus the optical flow (see Figure 1 for a detailed example).

The second descriptor, STIPs, was proposed in [START_REF] Laptev | On space-time interest points[END_REF]. This descriptor consists on the histogram of gradients concatenated to the histogram of optical flow (HOG-HOF), applied at Harris 3D interest points. Notice that, while the first descriptor uses stereo-vision and has low dimension (five), the second uses monocular vision and is high dimensional (200). Furthermore, while the former has a semi-dense nature, the latter is sparse in the spatiotemporal domain.

The Speech Descriptor

The auditory stream is represented by the Mel Frequency Cepstral Coefficients (MFCC). Widely used for speech and sound recognition (see [START_REF] Ramasubramanian | Continuous audio analytics by HMM and viterbi decoding[END_REF][START_REF] Lawrence | Theory and Applications of Digital Speech Processing[END_REF]), the MFCC are computed following the three steps: (i) perform the short-time Fourier transform (STFT), (ii) map the power spectrum onto the Mel scale and (iii) take the discrete cosine transform of these mapped powers. The are three main parameters associated with MFCC features. First, the frame size defines the length of the STFT (denoted by W ). Second, the frame shift (F ) determines the time between two consecutive STFT windows. Third, the amount of cepstral coefficients (D), that sets the dimension of the output MFCC representation.

(a) (b) (c) (d) (e) (f)

Bag-of-Words Gesture Learning

The Bag-of-Words (BoW) paradigm has shown to have very high performance in model learning from low-level features. Following [START_REF] Laptev | On space-time interest points[END_REF], it requires to: 1. Collect a set of local descriptors (possibly associated with interest points) of all training gesture instances.

2. Group these descriptors into K clusters.

3. Quantize all the descriptors assigning the label corresponding to its nearest cluster centroid. This provides for the "words".

4. Represent a gesture instance as a K-bin histogram of the quantized descriptor ("bag of words").

5. Train a classifier with these histograms in order to construct the gesture model. In steps 1-3, the word vocabulary (or the codebook) is constructed. The amount of local descriptors is typically high. Clustering these descriptors allows to represent an instance as a histogram, which in turn provides for a compact representation with fixed (chosen) length K. During steps 4-5, these compact representations together with their annotated labels are used to train a classifier. The BoW representation encodes the relative frequency of occurrences of the quantized descriptors, which discriminates among command classes.

We use the BoW paradigm to build auditory and visual models for each of the commands. Hence we have both a visual and an auditory classifier. When an instance of an unknown command class has to be recognized, the auditory and visual representations are computed and sent to their respective classifiers. The outputs of the two classifiers are fused to perform audio-visual gesture recognition (as explained in section 3.5).

Temporal Gesture Learning

The standard BoW approach is notoriously lacking a temporal model. Indeed, the order of the extracted features is not taken into account when learning the command model. A powerful statistical learning tool to model time series are hidden Markov models (HMMs). These models assume the existence of a discrete-valued hidden variable that describes the state of the sequence. At each of these states, the model is assumed to generate observations following different probability distributions (called emission distributions).

To fully specify such an HMM, two choices have to be done:

1. The amount of hidden states and the possible transitions among them.

2. The emission probability distributions at each of the states.

An expectation-maximization (EM) algorithm is derived in order to estimate the posterior distributions of the hidden state sequence, the transition probabilities, and the parameters of the emission distributions (see [START_REF] Lawrence | Theory and Applications of Digital Speech Processing[END_REF][START_REF] Bishop | Pattern Recognition and Machine Learning (Information Science and Statistics)[END_REF] for more details).

Audio-Visual Command Recognition

Once the auditory and visual classifiers are trained, we are ready to fuse the information coming from both classifiers. Let ac(g) and vc(g) denote the score of the command instance g to belong to class c given by the auditory and visual classifiers respectively. In order to combine the information from both classifiers we train a combined classifier consisting on (i) whitening the training data (unimodal classifier scores) and (ii) apply a weighting function.

The whitening procedure consists on computing the mean (µa) and the standard deviation (σa) of the auditory classification scores {ac(gn)} N,C n=1,c=1 , being N the number of command instances. A new auditory score is computed as ãc(g) = ac(g)-µa σa . The same procedure is applied to the visual classification scores.

Finally, the combined score is the result of a convex combination of the two whitened scores:

m l c (g) = lṽc(g) + (1 -l)ãc(g).
The value of l determines the trust we put on each modality. Actually, some cases deserve a special mention: l = 0 is equivalent to audio-based classification. l = 0.5 the auditory and visual scores stand on equal foot, l = 1 is equivalent to vision-based classification, and In general, l > 0.5 means that we put more trust on the visual classification score, whereas l < 0.5 means that we do it with the auditory score. This way of combining the two classifiers allows us to evaluate the relative trust we put on the modalities. The final classification is:

c * = arg max c m l c (g).

EXPERIMENTAL VALIDATION

In this section we describe the experiments done using different classification methods, such as HMMs and SVMs. The details are given below.

Implementation

The visual descriptors were computed using for STIPs the default vaules given in the binary code 2 and for Sceneflow using a disparity margin around 5 pixels with the intention to include most of the body of the actor and exclude the rest of the scene (see [START_REF] Sanchez-Riera | Action recognition robust to background clutter by using stereo vision[END_REF] for other details). The parameters to compute the MFCC features were set to the standard ones in speech recognition: W = 21.3 ms, f = W/2 and D = 13.

We tried the two different learning schemes (non-temporal and temporal) with audio and video features.

sSVM The classifier of the histogram of the scene-flow visual descriptors is a SVM. All the visual descriptors are grouped into K = 500 clusters. For each command clip, an histogram of centroid occurrences is computed and normalized. This is used to train the SVM.

lSVM We use the same method as before but on features from [START_REF] Laptev | On space-time interest points[END_REF].

2 http://www.di.ens.fr/ laptev/ aSVM The classifier of the histogram of auditory features is a SVM. We trained the SVM using histograms of size K = 500 in that case. Since the speech information is sparse in time, a voice activity detector is needed, we chose to use the state-of-the-art voice activity detector available in VOICEBOX [START_REF] Brookes | VOICEBOX: Speech processing toolbox for MATLAB[END_REF].

sHMM The per-frame histogram of visual features train temporal models (HMM). Five states per commands and mixture of five Gaussian components were trained. The size of the histogram was K = 20.

aHMM A temporal model is added to the auditory features. In that case eight states per command and Gaussian emission probabilities were trained.

lHMM This options was not explored because of the nature of the descriptor. As explained before, the HOG-HOF descriptor is sparse. Actually there are many image frames without any descriptor. Hence this descriptor is not well suited for a per-frame representation.

The Ravel Data Set

The experimental validation is performed on the "Robot Gestures" scenario of the Ravel data set [START_REF] Alameda-Pineda | Ravel: An annotated corpus for training robots with audio visual abilities[END_REF]. We use the eight sequences proposed for the Challenge, each one containing three instances of nine command categories. The set of voice-and-gesture commands are the followings: (i) wave ("Hello!"), (ii) walk towards the robot ("I am coming."), (iii) walk away from the robot ("Bye!"), (iv) stop handwave ("Stop!"), (v) turn around gesture ("Turn around."), (vi) come here gesture ("Come here."), (vii) point action ("Look!"), (viii) head motion for yes ("Yes") and (ix) head motion for no ("No"). In all cases, the human gesture is accompanied by speech corresponding to the gesture, shown above in brackets. Notice that all the actors in the data set are non-native English speakers of five different nationalities, hence there is a large variability in the speech commands.

Performance Results

In order to properly validate the method, the leave-oneout strategy is applied. This strategy is a cross-validation strategy that consists on selecting one of the sequences as test data and perform the training with the rest. Applying this to all the sequences makes the presented results statistically significant.

Evaluating multicategory classifiers means providing the confusion matrices. The ij-th entry of such matrix contains how many instances of the i-th class have been classified as class j. By averaging the elements of the diagonal, one obtains the average recognition rate (ARR) of the classifier.

Visual Categorization

Figures 3a, 3b and3c show the confusion matrices for the sSVM, the sHMM and the lSVM classifiers. Notice that sSVM performs very well in five out of nine gestures (Hello, I am coming, Look, No and Bye), well in two gestures (Yes and Come here) and poorly in two gestures (Stop and Turn around ). However, the sHMM classifier performs poorly compared to sSVM. Indeed, it gets good results for most of the actions, very good results for just two of them (Hello and Bye) and poor results in three gestures ( around). Notice also that there is no much difference between the sSVM and the lSVM classifiers. Actually, they mainly differ in two of the gestures No and Turn around.

Auditory Categorization

The confusion matrix of aSVM and aHMM can be seen in figures 4a and 4b respectively. The SVM-based classifier performs very good in six out of nine actions, good in two of them (Bye and Stop) and poorly just in one command (Turn around ). However the aHMM classifier has very good performance everywhere except for the Bye and Turn around commands.

Audio-Visual Categorization

It is worth to notice that, for instance, some actions that are difficult to recognize by sSVM as Stop or Come here are easily classified by aSVM, and viceversa with the actions Bye or Hello. This supports the idea that the combined classifier should outperform both unimodal classifiers. Figure 5 show the ARR of the combined classifier as a function of the weighting parameter l. Please remark that when using the same underlying model (curves aSVM-sSVM, aSVM-lSVM and aHMM-sHMM) the maximum performance of the combined classifier is achieved for values of l around 0.5. However, when the temporal classifier is combined with a non-temporal classifier, the maximum performance of the combined tends to shift towards the modality with temporal modeling. Furthermore, the temporally modeled classifier performs much better when some global information (coming from the non-temporal classifier) is taken into account. In that sense, it is worth to notice that, except for one case (left side of aSVM-sHMM), the combined classifiers outperform the unimodal ones when a small weight is given to the other modality's classifier.

Finally, Figures 6a-6f show the confusion matrix of all the multimodal classifiers for the optimal weighting parameter l. Generally speaking, we notice that the performance of these combined classifiers improved with respect to the unimodal ones. Most of them obtain outstanding results for some of the commands and very good results for most of the commands. We need to remark that the aHMM-sSVM and the aHMM-lSVM classifiers achieve an ARR of 77%, that represents a considerable increment respect the ARR computed on each modality independently.

CONCLUSIONS & FUTURE WORK

In this work we presented a method to recognize human commands from auditory and visual input signals. Based on a high-performance and solid learning technique, the method uses binocular/monocular visual features and monaural auditory features. Audio-visual recognition is performed at a later stage in which the classification scores from audio and video are combined to get the final audio-visual score, yielding to important increasings on the ARR. The method is tested in leaving-one-out fashion on a publicly available data set. Results show the importance of using both modalities for recognition and leave some room for improving.

More tests can be done trying different type of descriptors for visual and auditory data, as well as, try to fuse these data at different stages.

Acknoledgments

This work was supported by the European project HUMAVIPS FP7-ICT-2009-247525, http://humavips.eu/

CHALLENGE ORGANIZATION

D-META (Data sets for Multimodal Evaluation of Tasks and Annotations

) is the answer that a group of researchers gave to the call for Grand Challenges at the International Conference on Multimodal Interaction (ICMI) 2012. The aim of this challenge is to group researchers working in the same multimodal applications and to encourage them to test their algorithms on a common data set. Obviously, this will set a baseline on the quality of the current methods, the usability of such data sets as well as some hints about the quality of the provided data and annotations.

We would like to remark that the spirit of D-META is far from pursuing a standardization or universal agreement on how to process, present or annotate multimodal data. Instead, we would like to have a deep understanding of the underlying problems related to the use of multimodal data tar-geting different tasks. Where do the algorithms fail? Why? Which algorithm is best suited for each task?

In order to implement those thoughts we selected five different tasks for which we had associated data sets. The tasks were defined by a few researchers working on the field, and they covered a variety of topics: (i) recognition of gestures addressed to the robot by means of the vision and the audio, (ii) detection, localization and tracking of multiple speakers using audio-visual information, (iii) estimate of the level of engagement in a video-mediated communication, (iv) recognition of conversational gestures in first encounter dialogues and (v) recognition of feedback gestures in first encounter dialogues. A call-for-papers was enthusiastically written accordingly to the spirit of the challenge, delineating the proposed tasks and their evaluation procedures, referring to the associated data sets and fitting the quality and logistic criteria of the holding conference, ICMI.

Once the first boost was given and the structure was set, the data set needed to be made public. The data and annotation formats had to be specified too. In addition, the access to this data set had to be free, such that any researcher in the world could participate to the Challenge. Notice that a Golden Standard (the annotation of the Ground Truth) had to be provided for evaluation and transparency purposes.

In parallel, a crucial task for the success of such Challenge needed to be done: the advertising. We informed the community about what we were proposing. This is a difficult task because the communication channels are many and not always easy to find and use. Several mailing lists and the ICMI'12 website were used to make the event visible.

We encountered several problems through the organization of this Challenge. First of all, the process of publishing a data set is very resource-and time-consuming. Lots of efforts were made to set a proper access to the data set and their annotations. Secondly, we had also troubles advertising the Challenge. Even if we used all our communication channels, we had almost no external submissions. We think that several factors may have influenced this sad outcome:

• People are not always willing to test their methods to

Figure 1 :

 1 Figure 1: Construction of the proposed descriptor. The actor's face is detected from the left input image (a). The raw disparity map (b) is segmented, such that all pixels having the lower disparity than the actor's face are discarded (c). The descriptor is then computed for all remaining pixels undergoing non-zero motion, such that it consists of the pixel's position relative to the face, it's disparity (d), and horizontal (e) and vertical (f) components of optical flow.

Figure 2 :

 2 Figure 2: Mel Frequency Cepstral Coefficients for one voicecommand instance. From the raw signal (top) the voice activity is detected (middle) and used to mask the extracted MFCC (bottom).

Figure 3 :Figure 5 :

 35 Figure 3: Confusion matrix of the three visual classifiers: (a) sSVM, (b) sHMM and (c) lSVM.

Figure 4 :

 4 Figure 4: Confusion matrix of the two auditory classifiers: (a) aSVM and (b) aHMM.

  Yes, Stop and Turn

					Confusion Matrix (sSVM)
	hello	70.83	0.00	4.17	0.00	0.00	20.83	0.00	4.17	0.00
	i am coming	4.17	66.67	0.00	4.17	12.50	4.17	0.00	0.00	8.33
	look	0.00	0.00	70.83	0.00	0.00	4.17	12.50	0.00	12.50
	yes	0.00	8.33	0.00	54.17	0.00	8.33	0.00	0.00	29.17
	no	4.17	8.33	0.00	0.00	79.17	0.00	4.17	0.00	4.17
	stop	29.17	0.00	8.33	4.17	0.00	33.33	12.50	0.00	12.50
	turn around	0.00	12.50	20.83	0.00	4.17	8.33	33.33	4.17	16.67
	bye	0.00	0.00	0.00	8.33	0.00	0.00	12.50	75.00	4.17
	come here	0.00	12.50	0.00	4.17	4.17	8.33	8.33	8.33	54.17
		hello	i am coming	look	yes	no	stop	turn around	bye	come here
					(a) sSVM
		75.00	0.00	0.00	0.00	0.00	16.67	0.00	4.17	4.17
		0.00	58.33	8.33	12.50	4.17	0.00	4.17	4.17	8.33
		4.17	16.67	45.83	4.17	8.33	0.00	4.17	0.00	16.67
		0.00	25.00	8.33	29.17	4.17	0.00	0.00	4.17	29.17
		0.00	20.83	0.00	4.17	50.00	4.17	4.17	12.50	4.17
		33.33	12.50	12.50	0.00	0.00	16.67	12.50	0.00	12.50
		0.00	12.50	20.83	4.17	0.00	0.00	29.17	4.17	29.17
		0.00	12.50	0.00	4.17	4.17	4.17	0.00	66.67	8.33
		0.00	29.17	8.33	8.33	4.17	4.17	4.17	4.17
		hello	i am coming	look	yes	no	stop	turn around	bye	come here

Datasets for Multimodal Evaluation of Tasks and Annotations: http://d-meta.inrialpes.fr Organizers: Xavier Alameda-Pineda (PhD candidate at Perception Team, at INRIA Rhône-Alpes), Kristiina Jokinen (adjunct Professor of language technology at the University of Helsinki) Dirk Heylen (Professor at University of Twente), Roman Bednarik (interaction technology researcher at the University of Eastern Finland) and Michal Hradiš (PhD candidate at Faculty of Information Technology, Brno University of Technology)

other data sets/data type without having any guarantees of good performance or a possible award.

• Adapt the software/method to a new data set is not always straightforward, and that may be an obstacle for some researchers.

• Since the Challenge just started (it's its first edition), some people may not find this attractive enough or not even a good way to publish their research, even if we think it is worth.

• We should also consider the possibility that the community is simply not interested in what we propose.

For next editions of the Challenge, we should take into account two important points. On one hand, be aware that the processing and formatting of a data set is a long and tedious task, even for small data sets. The format has to be chosen carefully, since it may constrain the access to the data. On the other hand, take advance on the advertising campaign. The larger the amount of people knows about it, the higher the number of participants and the richer the discussion. We hope to be able to overcome these difficulties for the oncoming editions of D-META.