Jerome Brachat

Pierre Comon
email: pcomon@unice.fr

Bernard Mourrain

Elias Tsigaridas

SYMMETRIC TENSOR DECOMPOSITION

We present an algorithm for decomposing a symmetric tensor of dimension n and order d as a sum of of rank-1 symmetric tensors, extending the algorithm of Sylvester devised in 1886 for symmetric tensors of dimension 2. We exploit the known fact that every symmetric tensor is equivalently represented by a homogeneous polynomial in n variables of total degree d. Thus the decomposition corresponds to a sum of powers of linear forms.

The impact of this contribution is two-fold. First it permits an efficient computation of the decomposition of any tensor of sub-generic rank, as opposed to widely used iterative algorithms with unproved convergence (e.g. Alternate Least Squares or gradient descents). Second, it gives tools for understanding uniqueness conditions, and for detecting the tensor rank.

INTRODUCTION

Symmetric tensors show up in applications mainly as high-order derivatives of multivariate functions. For instance in Statistics, cumulant tensors are derivatives of the second characteristic function [START_REF] Mccullagh | Tensor Methods in Statistics, Monographs on Statistics and Applied Probability[END_REF]. Tensors have been widely utilized in Electrical Engineering since the nineties, because of the use of High-Order Statistics [START_REF] Swami | Multichannel ARMA processes[END_REF] [3] [4] [START_REF] Chevalier | Optimal separation of independent narrow-band sources -concept and performance[END_REF] [6] [START_REF] De Lathauwer | Tensor-based techniques for the blind separation of ds-cdma signals[END_REF]. Even earlier in the seventies, tensors have been used in Chemometrics [START_REF] Bro | Parafac, tutorial and applications[END_REF] or psychometrics [START_REF] Kiers | An efficient algorithm for Parafac of three-way data with large numbers of observation units[END_REF]. Another important application field is Data Analysis. For instance, Independent Component Analysis was originally introduced for symmetric tensors whose rank did not exceed dimension [START_REF] Comon | Independent Component Analysis[END_REF] [START_REF] Cardoso | Blind signal separation: statistical principles[END_REF]. Now, it has become possible to estimate more factors than the dimension. Further references may be found in [START_REF] Comon | Symmetric tensor decomposition[END_REF] [START_REF] Brachat | Symmetric tensor decomposition[END_REF], and numerous applications of tensor decompositions may be found in [START_REF] Cichocki | Adaptive Blind Signal and Image Processing[END_REF] [START_REF] Smilde | Multi-Way Analysis[END_REF].

The goal of this paper is to devise an algebraic technique able to decompose a symmetric tensor of arbitrary order and dimension in an essentially unique manner (i.e. up to scale and permutation) into a sum of rank-one terms. Of course, reaching such a goal requires some conditions, in particular related to its rank, which must be sub-generic. Our algorithm could be seen as an extension of the SVD algorithm from matrices to n-way arrays. We exploit the strong connection of symmetric tensors and homogeneous polynomials. This approach allows us to use effective algebraic geometry techniques, and to tackle the problem of decomposition using Veronese varieties, duality of vector spaces, and algorithms for polynomial system solving. To the best of our knowledge this is the first time that a decomposition algorithm for symmetric tensors is presented.

The rest of the paper is structured as follows: In the remaining of the section we present some historical remarks and we shed light to the connection of symmetric tensors and homogeneous polynomials. Sec. 2 presents Sylvester's approach for the binary case. In Sec. 3 we exploit many different, albeit equivalent, algebraic formulations of the decomposition problem, as well as the necessary algebraic tools. In Sec. 4 we present the algorithm and illustrate it with an example.

Historical remarks

Despite their obvious practical interest, numerical algorithms presently used in most scientific communities are suboptimal, in the sense that they either do not fully exploit symmetries [START_REF] Albera | Blind identification of overcomplete mixtures of sources (BIOME)[END_REF], minimize different successive criteria sequentially [START_REF] Yeredor | Non-orthogonal joint diagoinalization in the LS sense with application in blind source separation[END_REF] [START_REF] De Lathauwer | Fourth-order cumulant-based identification of underdetermined mixtures[END_REF], or are iterative and lack a guarantee of global convergence [START_REF] Harshman | PARAFAC: Parallel factor analysis[END_REF] [START_REF] Paatero | The multilinear engine: A table-driven, least squares program for solving multilinear problems, including the n-way parallel factor analysis model[END_REF]. In addition, they often request the rank to be much smaller than generic [START_REF] Jiang | Kruskal's permutation lemma and the identification of CANDE-COMP/PARAFAC and bilinear models[END_REF]. For estimating the mixing matrix in independent component analysis, based on fourth-order cumulant tensor, we refer the reader to [START_REF] De Lathauwer | Fourth-order cumulant-based identification of underdetermined mixtures[END_REF].

On the other hand, the algorithm based on Sylvester's theorem [START_REF] Kung | The invariant theory of binary forms[END_REF], recalled in section 2, provides a complete answer to the questions of uniqueness and computation, for any order [START_REF] Comon | Decomposition of quantics in sums of powers of linear forms[END_REF]. However, the latter is devoted to 2-dimensional symmetric tensors, and techniques based on pairwise processing have a very limited range of use when the rank exceeds the dimension.

The algorithm proposed in this paper is inspired from Sylvester's theorem, and extends its principle to larger dimensions. In addition, it fully exploits symmetry, and when the solution is essentially unique, it provides the decomposition for any sub-generic rank.

Tensors and Polynomials

Any symmetric tensor of dimension n, i.e. the range of each index, and order d, i.e. the number of indices, can be associated with a homogeneous polynomial in n variables of degree d. For instance, a third order tensor T ijk can be associated with the polynomial ijk c(i, j, k) T ijk x i x j x k , where c(•) denotes some fixed symmetric function. See e.g. [START_REF] Comon | Decomposition of quantics in sums of powers of linear forms[END_REF] for further details. We consider a homogeneous polynomial F h (x)

F h (x) = j0+j1+•••+jn=d a j0,j1,...,jn x j0 0 x j1 1 • • • x jn n . ( 1
)
Our goal is to compute a decomposition of F h as a sum of dth powers of linear forms,

F h (x) = r i=1 λ i (k i,0 x 0 + • • • + k i,n x n ) d = λ 1 k 1 (x) d + λ 2 k 2 (x) d + • • • + λ r k r (x) d , where λ i ∈ (
C, such that r is the smallest possible. This smallest r is often referred to as the tensor rank, or sometimes the polynomial width.

Let's see now how the decomposition is made possible in the case of homogeneous polynomials in two variables only.

THE BINARY CASE

Let p(x 1 , x 2 ) = d i=0 d i c i x i 1 x d-i 2 be a homogeneous polynomial of degree d in 2 variables. Denote H[r] the Hankel matrix of dimensions (d -r + 1) × (r + 1) with entries H[r] ij = c i+j-2 : H[r] =   c 0 c 1 • • • c r . . . . . . c d-r • • • c d-1 c d   .
Then we have:

Sylvester, 1886 p(x 1 , x 2 ) can be written as a sum of d th powers of r distinct linear forms in ( C as:

p(x 1 , x 2 ) = r j=1 λ j (α j x 1 + β j x 2 ) d , (2)
if and only if (i) there exists a vector q of dimension r + 1, with components q ℓ , such that

H[r] q = 0. (3)
and (ii) the polynomial q(x 1 , x 2 ) = r ℓ=0 q ℓ x ℓ 1 x r-ℓ 2 admits r distinct roots, i.e. it can be written as q(x 1 , x 2 ) = r j=1 (β * j x 1α * j x 2 ). It turns out that the proof of this theorem is constructive [START_REF] Sylvester | Sur une extension d'un théorème de Clebsch relatif aux courbes du quatrième degré[END_REF] [23] [START_REF] Comon | Symmetric tensors and symmetric tensor rank[END_REF] and yields the algorithm below.

1. Initialize r = 0 2. Increment r ← r + 1 3. If the row rank of H[r] is full, then go to step 2 4. Else compute a basis {k 1 , . . . , k l } of the right kernel of H[r].

Specialization (pick a random vector in the kernel):

• Take a generic vector q in the kernel, e.g.

q = i µ i k i • Compute the roots of the associated polyno- mial q(x 1 , x 2 ) = r ℓ=0 q ℓ x ℓ 1 x d-ℓ 2 . Denote them (β j , -α j ), where |α j | 2 + |β j | 2 = 1.
• If the roots are not distinct in P 2 , try another specialization. If distinct roots cannot be obtained, go to step 2. • Else if q(x 1 , x 2 ) admits r distinct roots then compute the coefficients λ j , 1 ≤ j ≤ r, by solving the linear system below

      α d 1 . . . α d r α d-1 1 β 1 . . . α d-1 r β r α d-2 1 β 2 1 . . . α d-1 r β 2 r : : : β d 1 . . . β d r       λ =      c 0 c 1 c 2 : c d      . 6. The decomposition is p(x 1 , x 2 ) = r j=1 λ j k j (x) d , where k j (x) = (α j x 1 + β j x 2 ).
Note that step 5 is a specialization only if the dimension of the right kernel is larger than one, which will not occur for ranks smaller than generic.

The goal is now to extend this kind of numerical algorithm to polynomials in more variables. This problem was open until now.

PROBLEM FORMULATIONS

Notation. If a = (a 1 , . . . , a n ) is a vector in IN n , then |a| is the sum of its elements, i.e. |a| = n i=a a i . By x a will denote the monomial x a1 1 • • • x an n . Let R be the ring of polynomials ( C[x 1 , . . . ,
x n ], while R d will denote the vector space of polynomials of (total) degree at most d. The set elements. The corresponding basis of the dual space R d , that is the set of linear forms that compute the coefficients of a polynomial in the primal basis, is the set {d a } |a|≤d , where d a : R d → ( C and

{x a } |a|≤d = {x a1 1 • • • x an n } a1+•••
d a (f ) = 1 Q n i=1 ai! d a 1 dx1 • • • d an dxn f (0)
. The superscript h denotes the homogeneous version of the polynmomial. Let S be the set of homogeneous polynomials in n + 1 variables. S d represents the homogeneous polynomials of degree d, and P(S d ), the corresponding projective space. Similarly interpretations hold for the dual spaces S and S d . Analogous to the affine case, we can define primal and dual bases for the homogeneous case.

Direct approach by polynomial fitting

The first idea is merely to solve, in a given polynomial basis, the polynomial system, induced by the equation

F h (x) - i k i (x) d = 0,
with respect to the coefficients of the linear forms k i . We call this the direct approach. In the tensor framework, even if the rank is supposed to be known, attempts to solve this problem have not entailed efficient algorithms (cf. section 1). In the polynomial framework, it is easy to see that we end up with an over-determined polynomial system of n+d d equations in r(n + 1) unknowns. This description of the problem is not optimal, since it introduces r! redundant solutions corresponding to permutations of the linear forms. Another drawback is that polynomials involved are of degree d in the coefficients k i,n , which are too high from the computational point of view. In fact, our approach does not involve the solution of polynomial systems of degree higher than 2.

Different views using duality

We consider the following Veronese map of degree d

ν : S1 → S d k(x) → k(x) d .
which sends a linear (homogeneous) polynomial to its d-th power. Recall that the (monomial) basis of S 1 is the set {x 0 , x 1 , . . . , x n }, while the basis of

S d is the set {x a } |a|≤d = {x a0 0 x a1 1 • • • x an n } a1+•••+an=d , viz.
the set of all the monomials in x 0 , x 1 , . . . , x n of total degree d. The cardinality of the basis is

(x) = k 0 x 0 + • • • + k n x n corresponds to k(x) d = i0+•••+in=d d io,...,in k i0 0 • • • k in n x i0 0 • • • x in n . In terms of vectors, k = [k 0 , . . . , k n ] ⊤ corresponds to the vector [. . . , d i0,...,in k i0 0 • • • k in n , . . . ] ⊤ . Another Veronese map, also of degree d, is δ : ( C n+1 → b S d z → ½z, which sends a point z = (z 0 , z 1 , . . . , z n ) to [. . . , z i0 0 • • • z in n , . . . ].
Recall that the linear functionals {d a } |a|=d consist a basis of S d . It holds that S 1 ∼ = ( C n+1 . The map τ :

τ : S d → S d [ . . . , z i0 0 . . . z in n , . . . ] ⊤ → [ . . . , d i0,...,in z i0 0 . . . z in n , . . . ]
⊤ is an isomorphism. The map remains an isomorphism even if we restrict it to the images of the maps ν and δ, that is ν(S 1 ) and δ( ( C n ), respectively. The inverse of τ is the map τ -1 : S d → S d . Consider a polynomial in S d , that is the d-th power of a linear form, say k(x) d . If we apply the map τ -1 to this polynomial, then we have that τ -1 (k(x) d ) = ½ k ; that is the linear form that gives the evaluation of a polynomial (homogeneous of degree d) over the point k = [k 0 , k 1 , . . . , k n ] ⊤ . Let us now revisit the problem of decomposition. Initially we are given a polynomial F h (x) ∈ S d . The decomposition F h (x) = r i=1 λ i k i (x) d corresponds to a secant variety in ν(S 1 ) ⊂ S d . Using the properties of the isomorphism τ and its inverse we can gain another view of the problem. If we apply τ -1 to F h we compute its dual, that is Φ = τ -1 (F h ). The decomposition of the latter, i.e. Φ =

r i=1 λ i ½ ki , is a linear combination of elements in δ( ( C n ) ⊂ R d . Overall, it holds that τ -1 (F h (x)) = τ -1 r i=1 λ i k i (x) d = r i=1 λ i τ -1 (k i (x) d ) = r i=1 λ i ½ ki = Φ. Moreover, τ (Φ) = τ ( r i=1 λ i ½ Li ) = r i=1 λ i τ (½ ki ) = r i=1 λ i k i (x) d = F h (x).
All the previous views of the decomposition problem are equivalent. The results obtained for any of them could be translated for the other.

Quotient algebra and duality

The idea of the algorithm is to exploit the properties of Φ ∈ R, that we assume that is known up to degree d. More precisely, we consider the symmetric bilinear form À Φ : (p, q) → Φ(p q), the matrix of which in the mono- mial basis is (Φ(x a+b )) a,b∈IN n . Let I Φ be the kernel of

À Φ . Proposition 3.2 If Φ = r i=1 λ i ½ ki with λ i = 0 and k i ∈ ( C n , then p ∈ I Φ iff p(k i ) = 0 for i = 1, . . . , r.
In other words, the common roots of all the polynomials in I Φ define the linear terms in the tensor decomposition of F . In order to compute the zeros of I Φ , we may use a well-known theorem (see e.g. [START_REF] Cox | Using Algebraic Geometry[END_REF][START_REF] Cox | Ideals, Varieties, and Algorithms[END_REF][START_REF] Elkadi | Introduction á la résolution des systémes polynomiaux[END_REF]), which we apply to the zero-dimensional ideal I Φ : Theorem 3.3 The eigenvalues of the matrices Å a and Å ⊤ a , of the linear operators that correspond to the multiplication by a in R modulo I Φ , and its transposed, are {a(k 1 ), . . . , a(k r )}. The common eigenvectors of the matrices (Å ⊤ xi ) 1≤i≤n are (up to a scalar) ½ ki , i = 1, . . . , r.

If we denote by À E

Φ the restriction of À Φ to a vector space E of dimension r on which À Φ is invertible, we have the relation À E a⋆Φ = Å ⊤ a À E Φ , where À a⋆Φ : (p, q) → Φ(a p q). Thus the solution of the generalized eigenvalue problem

(À E a⋆Λ -ζ À E Λ )v = Ç yields the eigenvec- tor À E Λ v of Å ⊤
a , which are by Th. 3.3, the evaluations ½ ki . From these eigenvectors, we deduce the linear fac- tors in the tensor decomposition. The coefficients λ i (i = 1, . . . , r) can then be computed by solving a linear system of size r.

ALGORITHM

The algorithm that we will present for decomposing a symmetric tensor as sum of rank 1 symmetric tensors generalizes the algorithm of Sylvester [START_REF] Sylvester | Sur une extension d'un théorème de Clebsch relatif aux courbes du quatrième degré[END_REF], devised for dimension 2 tensors, see also [START_REF] Comas | On the rank of a binary form[END_REF].

Overview

Algorithm 1: Symmetric tensor decomposition Input: A homogeneous polynomial f (x 0 , x 1 , . . . , x n ) of degree d.

Output: A decomposition of f as f = r i=1 λ i k i (x) d with r minimal. -Compute the coefficients of f * : c α = a α d α -1 , for |α| ≤ d, α = (α 1 , . . . , α n ); -r := 1; -Repeat 1. Compute a set B of monomials of degree ≤ d connected to 1 with |B| = r; 2. Find parameters h s.t. det(À B Φ ) = 0 and the operators Å i = À B xiΦ (À B Φ ) -1 commute. 3.
If there is no solution, restart the loop with r := r + 1. 4. Else compute the n × r eigenvalues ζ i,j and the eigenvectors v j , s.t. Å i v j = ζ i,j v j , i = 1, . . . , n, j = 1, . . . , r. until the eigenvalues are simple.

-Solve the linear system in (λ j ) j=1,...,k : Φ = r j=1 λ j 1 vj where v j ∈ K n are the eigenvectors found in step 4.

Let us briefly comment on the computation process. A basis connected to 1, is a basis containing 1 where each element different from 1, is the product of a variable by another element of the basis. Consider the homogeneous polynomial f (x) in (1) that we want to decompose. We may assume without loss of generality, that for at least one variable, say x 0 , all its coefficients in the decomposition are non zero, i.e. k i,0 = 0, for 1 ≤ i ≤ r. We dehomogenize f with respect to this variable and we denote this polynomial by f a := f (1, x 1 , . . . , x n ). We want to decompose the polynomial f a (x) ∈ R d as a sum of powers of linear forms, i.e.

f (x) = r i=1 λ i (1 + k i,1 x 1 + • • • + k i,n x n ) d = r i=1 λ i k i (x) d .
Equivalently, we want to decompose its corresponding dual element f * ∈ R d as a linear combination of evaluations over the distinct points

k i := (k i,1 , • • • , k i,n ): f * = r i=1 λ i ½ ki (
we refer the reader to the end of Sec- tion 3.2).

Assume that we know the value of r. If we know the value of Φ on polynomials of degree high enough, it allows us to compute the tables of multiplication modulo the kernel of À Φ . By Theorem 3.3, if we solve the generalized eigenvector problem (À x1⋆Φζ À Φ )v = Ç, then we can recover the points of evaluation k i . By solving a linear system, we will then deduce the value of λ i , . . . , λ r . For certain (big) values of r it can happen that not all the elements of the corresponding matrices are known. In this case, we use the property that the matrices of multiplication commute, and we form a system, the solutions of which are these unknown elements. We refer the reader to [START_REF] Comon | Symmetric tensor decomposition[END_REF][START_REF] Brachat | Symmetric tensor decomposition[END_REF] for details, and we present an example to illustrate the algorithm.

Example

(1) Convert the symmetric tensor to the corresponding homogeneous polynomial.

Assume that we are given a tensor of dimension 3 and order 5, and that the corresponding homogeneous polynomial is f = -1549440 x0x1x2 (2) Compute the matrix of the linear form.

We form a n+d-1 d × n+d-1 d matrix, the rows and the columns of which correspond to the evaluation of the dual of the polynomial over all the monomial {x a } |a|≤d , using the map a j0 j1 ... jn → a j0 j1 ... jn j0! j1! ... jn! d!

, where a j0 j1 ... jn is the coefficient of the monomial x j0 1 • • • x jn n . Part of the corresponding matrix follows.

       1 x 1 x 2 x 2 1 x 1 x 2 x 2 2 x 3
      
The whole matrix is 21 × 21. For reasons of space we present only the first 7 columns. Notice that we do not know the elements in some positions of the matrix (in the 7th column). In general we do not know the elements that correspond to monomials with (total) degree higher than 5.

(3) Extract a principal minor of full rank.

We should re-arrange the rows and the columns of the matrix so that the first principal minor is of full rank, r. We call this minor À Φ . In order to do that we try to put the matrix in row echelon form, using elementary row and column operations. In our example the 4 × 4 principal minor is of full rank, so there is no need for re-arranging the matrix. The matrix À Φ is

À Φ =    38 -24 36

  

Notice that the columns of the matrix correspond to the monomials {1, x 1 , x 2 , x 2 1 }.

(4) We compute the "shifted" matrix À x1⋆Φ .

If the columns of À Φ correspond to set of some mono- mials, say {x a }, then the columns of À x1⋆Φ correspond to the set of monomials {x 1 x a }.

, x 2 1 , x 1 x 2 , x 3 1 }, i.e.
the monomials of À Φ , {1, x 1 , x 2 , x 2 1 }, multiplied by x 1 . We assume for the moment that all the elements of the matrices À x1⋆Φ and À Φ are known. If this is not the case, then we can compute the unknown entries of the matrix, using necessary and sufficient conditions of the quotient algebra; it holds that Å xi Å xj -Å xj Å xi = Ç.

We refer the reader to [START_REF] Comon | Symmetric tensor decomposition[END_REF] for details.

(5) We solve the equation

(À x1⋆Φ -ζ À Φ )v = 0.
We solve the generalized eigenvalue/eigenvector problem using one of the well-known techniques [START_REF] Golub | Matrix Computations[END_REF].

We multiply the (generalized) eigenvectors by À Φ and we normalize the resulting vectors so that the first element is 1, and we read the solutions from the coordinates of the (normalized) eigenvectors, according to Th.

  

The coordinates of the eigenvectors correspond to the elements {1, x 1 , x 2 , x 2 1 }. Thus, we can recover the coefficients of x 1 and x 2 in the decomposition from coordinates of the eigenvectors. Thus, The polynomial admits a decomposition f = λ 1 (x 0 -12x 1 -3x 2 ) 5 + λ 2 (x 0 -2x 1 + 3x 2 ) 5 + λ 3 (x 0 + 2x 1 + 3x 2 ) 5 + λ 4 (x 0 + 12x 1 -13x 2 ) 5 .

It remains to compute λ i 's. We can do this easily by solving an over-constrained linear system, which we know that always has a solution, since the decomposition exists. Doing that, we deduce that λ 1 = 5, λ 2 = 15, λ 3 = 15 and λ 4 = 3. We obtain the following minimum decomposition of the polynomial as a sum of powers of linear forms: f = 5(x 0 -12x 1 -3x 2 ) 5 + 15(x 0 -2x 1 + 3x 2 ) 5 + 15(x 0 + 2x 1 + 3x 2 ) 5 + 3(x 0 + 12x 1 -13x 2 ) 5 that is the corresponding tensor is of rank 4.

CONCLUSIONS AND FUTURE WORK

We proposed an algorithm for symmetric tensor decomposition, extending the algorithm of Sylvester to dimensions higher than 2. The algorithm decomposes symmetric tensors when the rank is sub-generic and when the decomposition is unique. In order for the algorithm to work for any rank, we should be able to extend the quotient matrix defined in Sec. 3.3. We will report on this in the near future. We are currently working on an efficient C++ implementation of the algorithm.

  +an≤d represents the elements of the monomial basis of the vector space R d . It contains n+d-1 d

- 1 .

 1 Under the action of ν, a linear polynomial k

  3 + 2417040 x0x1 2 x2 2 + 166320 x0 2 x1x2 2 -829440 x0x1 3 x2 -5760 x0 3 x1x2 -222480 x0 2 x1 2 x2 + 38 x0 5 -497664 x1 5 -1107804 x2 5 -120 x0 4 x1 + 180 x0 4 x2 + 12720 x0 3 x1 2 + 8220 x0 3 x2 2 -34560 x0 2 x1 3 -59160 x0 2 x2 3 +831840 x0x1 4 +442590 x0x2 4 -5591520 x1 4 x2 + 7983360 x1 3 x2 2 -9653040 x1 2 x2 3 + 5116680 x1x2 4 .

  In our example

			-24	1272	-288	-3456	
	À x1⋆Φ	 	1272 -3456 -288 -7416	-7416 5544	166368 -41472	  ,
			-3456 166368 -41472 -497664
	the columns of which correspond to the monomi-
	als {x 1

We thank the anonymous referees for their comments that helped us improve the paper.

0 This work has been partially supported by contract ANR-06-BLAN-0074 "Decotes"


André L F De Almeida

Xavier Luciani
email: luciani@i3s.unice.fr.

Pierre Comon
email: pcomon@i3s.unice.fr.

BLIND IDENTIFICATION OF UNDERDETERMINED MIXTURES BASED ON THE HEXACOVARIANCE AND HIGHER-ORDER CYCLOSTATIONARITY

Keywords: Blind identification, underdetermined mixtures, cyclostationarity, hexacovariance. 1

In this work, we consider the problem of blind identification of underdetermined mixtures in a cyclostationary context relying on sixthorder statistics. We propose to exploit the cyclostationarity at higher orders by taking into account the knowledge of source cyclic frequencies in the sample estimator of the observation hexacovariance. Two blind identification algorithms based on the proposed estimator are considered and their performances are tested by means of computer simulations. Our simulation results show that significant improvements can be obtained when both second and fourth-order cyclo-stationarities are exploited.

INTRODUCTION

Blind identification and blind source separation methods have been successfully applied in multidisciplinary contexts, including radiocommunications, sonar, radar, biomedical signal processing and data analysis, just to mention a few. A widespread class of these methods rely on independent component analysis by means of higher-order statistics [START_REF] Comon | Independent Component Analysis, a new concept ?[END_REF]. In this context, a problem that have attracted a particular interest is that of blind identification of underdetermined mixtures (i.e. when we have more sources than sensors). Several solutions have been proposed in the literature to solve this problem (see, e.g. [START_REF] Cardoso | Super-symmetric decomposition of the fourth-order cumulant tensor. Blind identification of more sources than sensors[END_REF][START_REF] Comon | Blind identification and source separation in 2x3 underdetermined mixtures[END_REF][START_REF] Albera | Blind identification of overcomplete mixtures of sources (BIOME)[END_REF][START_REF] Comon | Blind identification of under-determined mixtures based on the characteristic function[END_REF][START_REF] De Lathauwer | Fourthorder cumulant-based blind identification of underdetermined mixtures[END_REF][START_REF] De Lathauwer | Blind identification of underdetermined mixtures by simultaneous matrix diagonalization[END_REF] and references therein). The solution proposed in these works may resort to second, fourth or sixth-order statistics of the output data.

In several applications such as radiocommunications and passive listening, the sources may be nonstationary, and they are often (quasi)-cyclostationary. This property appears as soon as the observations are oversampled and/or when the different sources have different bandwidth. The work [START_REF] Ferreol | On the behavior of current second and higher order blind source separation methods for cyclostationary sources[END_REF] addressed the behavior of secondand four-order blind source separation algorithms in a cyclostationary context. The authors proposed to exploit the cyclostationary property of the sources by adding a correction term to the standard sample estimator of the quadricovariance which takes into account the known cyclic frequencies of the received sources. The results presented in [START_REF] Ferreol | On the behavior of current second and higher order blind source separation methods for cyclostationary sources[END_REF] show that the performance of fourth-order statistics based blind source separation algorithms can be considerably improved when the proposed estimator is used. The cyclostationarity property has no yet been considered for cumulant estimators of orders higher than four. The works [START_REF] Albera | Sixth order blind identification of under-determined mixtures -BIRTH[END_REF][START_REF] Albera | Blind identification of underdetermined mixtures based on the hexacovariance[END_REF] exploit sixth-order statistics in the blind identification problem. Following these works, a family of blind identification algorithms called BIOME (Blind Identification of Overcomplete MixturEs) was proposed in [START_REF] Albera | Blind identification of overcomplete mixtures of sources (BIOME)[END_REF]. Although powerful, these algorithms do not take into account the cyclostationary nature of the sources since they rely on the standard sample estimate of the hexacovariance.

In this work, we propose to exploit the higher-order cyclostationarity in the blind identification problem. Motivated by the results of [START_REF] Ferreol | On the behavior of current second and higher order blind source separation methods for cyclostationary sources[END_REF][START_REF] Ferreol | Second-order blind separation of first-and second-order cyclostationary sources[END_REF] on one hand, and of [START_REF] Albera | Blind identification of overcomplete mixtures of sources (BIOME)[END_REF] on the other hand, we propose to take into account the known cyclic frequencies of the sources in the calculation of the empirical estimator of the hexacovariance. Two blind identification algorithms based on the proposed estimator are tested. The first one is the 6-BIOME algorithm of [START_REF] Albera | Blind identification of overcomplete mixtures of sources (BIOME)[END_REF], also referred to as "BIRTH" in the 6th order case, while the second one is a direct minimization of a tensor model fitting error by an iterative algorithm.

PROBLEM DEFINITION AND ASSUMPTIONS

Consider a noisy mixture of P statistically independent narrowband sources received by an array of M sensors. The vector y(k) containing the discrete-time version of the complex envelopes of the received signal at the sensor outputs can be modeled according to the following classical linear model:

y(k) = P p=1 sp(k)hp + n(k) . = Hs(k) + n(k), (1)
where

H = [h 1 , . . . , h P ] ∈ C M ×P , s(k) = [s 1 (
k), . . . , s P (k)] T ∈ C P and n(k) ∈ C M are the mixing matrix, source and noise random vectors, respectively. It is assumed that for any fixed time index k, s(k) and n(k) are statistically independent. We are interested in the case where the received source signals are cyclostationary and have a nonzero carrier residue. These properties are generally verified in interception or passive listening applications. In this scenario, the input-output model (1) may be too restrictive so that we adopt following observation model:

y(k) = P p=1 s p (k)e -(2π∆f p kT s +φ p ) h p + n(k) = Hs(k) + n(k), (2)
where s(k) = [s1(k), . . . , sP (k)] T is the new source signal vector, with s p (k) . = s p (k)e -(2π∆f p kT s +φ p ) , while ∆f p and φ p are, respectively, the carrier residue and phase of the p-th source and T s is the sampling period. The sources are assumed to have the same symbol period T , while the observations are sampled at the Nyquist rate, i.e. T s ≥ T /2.

Before to proceed, the following hypotheses are assumed:

H1. The sources s 1 (k), . . . , s P (k) are non-Gaussian, cyclostationary, cycloergodic, and mutually uncorrelated up to order 6;

H2. The noise vector n(k) is stationary and ergodic, following a complex-valued Gaussian distribution;

H3. The sixth order marginal source cumulants are not null and have all the same sign;

H4. The mixing matrix H does not contain collinear columns.

The problem is to identify the mixing matrix H (up to trivial column permutation and scaling) and, possibly, the source vector s(k) from the only knowledge of the observation vector y(k) by means of an estimation of its associated sixth-order statistics. Recall that we are interested in the so-called underdetermined case, which means that we have P > M .

After the identification of the mixing matrix, we consider a Maximum-A-Posteriori (MAP) criterion to estimate the source vector s(k) by means of an exhaustive search over the known finite alphabet of the sources. This methodology will be used for evaluating the Bit-Error-Rate (BER) of the proposed blind identification algorithms exploiting higher-order cyclostationarity.

HIGHER-ORDER STATISTICS IN THE PRESENCE OF CYCLOSTATIONARITY

As discussed in [START_REF] Ferreol | On the behavior of current second and higher order blind source separation methods for cyclostationary sources[END_REF][START_REF] Ferreol | Second-order blind separation of first-and second-order cyclostationary sources[END_REF], in the cyclostationary context the covariance function of s p (k) admits a Fourier series expansion over the set Γp = {αp} of cyclic frequencies, where the coefficients of this expansion correspond to the cyclic covariances of the p-th source. Therefore, the covariance matrix of y(k) contains cyclic frequencies of all the sources belonging to the set {Γ 1 , . . . , Γ P }. In the general case, the N -th order cyclic covariance of y(k) associated with a given cyclic frequency α is defined by:

c i p+1 ,...,i N i 1 ,...,i p ,y (α) = c i p+1 ,...,i N i 1 ,...,i p ,y (k)e -2παkT s d ( 3
)
where c i p+1 ,...,i N i 1 ,...,ip,y (k) are the N -th order output cumulants and

f (k) d . = lim K→∞ (1/K) K k=1
f (k) corresponds to the discretetime temporal mean operation of f (k) over an infinite number of samples. We assume that the set {Γ1, . . . , ΓP } of cyclic frequencies of the sources is known. It is worth noting that a set of these cyclic frequencies depends on the nonzero carrier residues of the sources ∆f 1 , . . . , ∆f P . The main motivation for considering the cyclostationarity property in the blind identification and source separation problems comes from the fact that, when the observations are (quasi)-cyclostationary, the time-averaged cumulants c i p+1 ,...,i N i 1 ,...,i p ,y (i.e. those calculated at the zero cyclic frequency only) generate, as K becomes infinite, an "apparent" (biased) estimation of the cumulants instead of the true cumulants. This is a consequence of the time dependence of the statistics of the data for (quasi)-cyclostationary sources [START_REF] Ferreol | On the behavior of current second and higher order blind source separation methods for cyclostationary sources[END_REF]. This estimation bias can, however, be compensated by adding a correction term to the sample estimator that takes into account all the nonzero cyclic frequencies present in the observed data. We propose to exploit the cyclostationarity property by means of sixth-order statistics (hexacovariance). It is worth mentioning that the hexacovariance has been considered in [START_REF] Albera | Blind identification of underdetermined mixtures based on the hexacovariance[END_REF][START_REF] Albera | Blind identification of overcomplete mixtures of sources (BIOME)[END_REF] for blind identification of the underdetermined mixtures. However, these works have not exploited the cyclostationarity property of the sources. Here, we rely on the results of [START_REF] Ferreol | On the behavior of current second and higher order blind source separation methods for cyclostationary sources[END_REF][START_REF] Ferreol | Second-order blind separation of first-and second-order cyclostationary sources[END_REF] (which have considered the quadricovariance) to address the performance of two blind identification algorithms based on the hexacovariance. These algorithms are presented in Section 4. Our goal is to evaluate the potential improvements obtained when higher-order cyclic moments are exploited. Proposed hexacovariance estimator: Let c i 4 ,i 5 ,i 6 i 1 ,i 2 ,i 3 ,y be an element of the hexacovariance of the observations y(k). The cyclic estimator c i 4 ,i 5 ,i 6 i 1 ,i 2 ,i 3 ,y of the hexacovariance based on K snapshots of the data is given by:

c i 4 ,i 5 ,i 6 i 1 ,i 2 ,i 3 ,y . = r i 4 ,i 5 ,i 6 i 1 ,i 2 ,i 3 (0) -[3] α,β α+β=0 r i 4 i 1 ,i 2 ,i 3 (α) r i 6 i 5 (β) -[9] α,β α+β=0 r i 4 ,i 5 i 1 ,i 2 (α) r i 6 i 3 (β) -[3] α,β α+β=0 r i 1 ,i 2 (α) r i 4 ,i 5 ,i 6 i 3 (β) + 2[9] α,β,γ α+β+γ=0 r i 1 ,i 2 (α) r i 4 i 3 (β) r i 5 ,i 6 (γ) + 2[6] α,β,γ α+β+γ=0 r i 4 i 1 (α) r i 5 i 2 (β) r i 6 i 3 (γ) (4)
where

r i 4 ,i 5 ,i 6 i 1 ,i 2 ,i 3 (0) . = 1 K K k=1 y i 1 (k)y i 2 (k)y i 3 (k)y * i 4 (k)y * i 5 (k)y * i 6 (k)
is the estimated sixth-order moment at the zero frequency and

r i p+1 ,...,i p+q i 1 ,...,i p (α) . = 1 K K k=1 yi 1 (k) • • • yi p (k)y * i p+1 (k) • • • y * i p+q (k)e -2παkT s .
is the estimated (p + q)-th order moment evaluated at the cyclic frequency α, where p + q ≤ 6. In the compact expression (4), [n] denotes the McCullagh bracket notation representing the existence of n monomials of the same order that arises by permuting separately either superscripts or subscripts [12] [13]. Finally, note that the values taken by α, β, and γ satisfying α + β + γ = 0, are those of the known cyclic frequencies. Remark: The proposed estimator of the hexacovariance is only "approximately unbiased" since we have ignored the bias introduced by the estimated sixth-order moment. This approximation has two motivations. First, the complexity associated with the calculation of the unbiased sixth-order moment estimator is prohibitive. Second, the performances obtained with the proposed approximation are satisfactory, as it will be shown in Section 5.

BLIND IDENTIFICATION ALGORITHMS

In this section, we exploit the higher-order cyclostationarity of the sources by considering two blind identification algorithms capable of identifying underdetermined mixtures. The first one is the sixthorder version of the BIOME family of algorithms and is called 6-BIOME [START_REF] Albera | Blind identification of overcomplete mixtures of sources (BIOME)[END_REF]. The second algorithm is based on the same modeling, but the 6th order cumulants are stored in a 3rd order tensor instead of a matrix. Before presenting these algorithms, we introduce some notation and properties associated with the matrix representations of the hexacovariance tensor. Thanks to the multilinearity proprety of cumulants, the hexacovariance c i 4 ,i 5 ,i 6 i 1 ,i 2 ,i 3 ,y of the observations y(k) = Hs(k) is a sixth-order rank-P tensor C ∈ C M ×M ×M ×M ×M ×M which can be written as

C y = P p=1 κ p (h p ⊗ ⊗ ⊗ h p ⊗ ⊗ ⊗ h p ⊗ ⊗ ⊗ h * p ⊗ ⊗ ⊗ h * p ⊗ ⊗ ⊗ h * p ), (5)
where ⊗ ⊗ ⊗ denotes the outer product between vectors and κp is the marginal 6th order cumulant of the p-th source. The latter model is sometimes referred to as "PARAFAC".

Symmetric matrix factorization

The overall information contained in the hexacovariance tensor Cy defined in ( 5) can be organized in a symmetric matrix C y ∈ C M 3 ×M 3 defined as follows [START_REF] Albera | Blind identification of overcomplete mixtures of sources (BIOME)[END_REF]:

Cy = P p=1 κp hp ⊗ hp ⊗ h * p hp ⊗ hp ⊗ h * p = H 3 ∆s(H 3 ) H , ( 6
)
where ⊗ and denote, respectively the Kronecker and Khatri-Rao products, H 3 = H H H * ∈ C M 3 ×P and ∆ s ∈ C P ×P is a diagonal matrix containing the marginal source cumulants κ 1 , . . . , κ P along the main diagonal. The 6-BIOME algorithm briefly, which is presented in Section 4.1, relies on model [START_REF] De Lathauwer | Fourthorder cumulant-based blind identification of underdetermined mixtures[END_REF].

Non-symmetric matrix factorizations

We can organize the information contained in the hexacovariance tensor in alternative (non-symmetric) matrix forms. Let us consider the following one:

C y = H H H H * H * H (3)T , ∈ C M 5 ×M (7)
or, alternatively,

C y = H (1) H (2) H (3)T , ( 8
)
where H (1) = H H H, H (2) = H * H * and H (3) = H * ∆ s . The factorization ( 8) is an equivalent third-order PARAFAC model of dimensions M 3 ×M2 ×M representing the hexacovariance tensor. In Section 4.2, this factorization is exploited for blind identification by means of the Levenberg-Marquardt (LM) algorithm.

Sixth-order BIOME (6-BIOME)

The 6-BIOME algorithm is reminiscent of the BIRTH (Blind Identification using Redundancies in the daTa Hexacovariance matrix) algorithm proposed in [START_REF] Albera | Sixth order blind identification of under-determined mixtures -BIRTH[END_REF] and later improved in [START_REF] Albera | Blind identification of underdetermined mixtures based on the hexacovariance[END_REF]. It exploits the multilinear algebraic structure of the hexacovariance by solving a joint approximate diagonalization problem based on the symmetrically unfolded matrix factorization (6) of the estimated hexacovariance (4). Following the idea of [START_REF] Albera | Blind identification of overcomplete mixtures of sources (BIOME)[END_REF][START_REF] Albera | Sixth order blind identification of under-determined mixtures -BIRTH[END_REF], we can write the square-root of (6) as:

(C y ) 1/2 = H 3 ∆ s V (9)
where V is a unitary matrix. The 6-BIOME algorithm is summarized as follows:

6-BIOME ALGORITHM 1. Estimate the hexacovariance tensor using all the cyclic moments (4) and form Cy;

2. Compute the square-root (Cy) 1/2 ∈ C M 3 ×P from the Eigen-Value Decomposition (EVD) of Cy;

3. Slice (Cy) 1/2 into M matrix blocks Γm ∈ C M 2 ×P ; 4. Find V by solving a simultaneous diagonalization problem from M (M -1)/2 Hermitian matrices Θ m,n . = Γ † m Γ n ; 5. Calculate H 3 = (C y ) 1/2 V H ; 6. Arrange each of the P columns of H 3 in vector b m ∈ C M 3 ; 7. Transform each vector bm ∈ C M 3 in a set of M matrices B ∈ C M ×M
and calculate the dominant eigenvector h m of each of these matrices;

8. Form H, the columns of which are the vectors h m .

Identification with Levenberg-Marquardt (PARAFAC-LM)

The second algorithm is based on the 3rd order tensor representation (7) of decomposition [START_REF] Comon | Blind identification of under-determined mixtures based on the characteristic function[END_REF]. The proposed approach consists of iteratively fitting this 3rd order storage of the hexacovariance tensor using the Levenberg-Marquardt (LM) method 2 .We consider the minimization of the following quadratic cost function:

f (p) = 1 2 e(p) 2 F = 1 2 e H (p)e(p), (10)
where e(p) = vec C y -( H (1) H (2) ) H (3)T ∈ C M 6 ×1 is the residue and p is the parameter vector defined as:

p =   p H (1) p H (2) p H (3)   =   vec( H (1)T ) vec( H (2)T ) vec( H (3)T )   ∈ C 3M P ×1 , ( 11
)
and the LM update is given as follows:

p(i + 1) = p(i) -J H (i)J(i) + λ(i)I -1 g(i), (12)
where J(i) denotes the Jacobian matrix:

J(i) = J H (1) (i) J H (2) (i) J H (3) (i) ∈ C M 6 ×3M P , ( 13
)
and g(i) denotes the gradient vector:

g(i) =   g H (1) (i) g H (2) (i) g H (3) (i)   ∈ C 3M P ×1 ( 14
)
After convergence, an estimate of the mixture H up to column permutation and scaling is obtained from the estimated parameter vector p H (3) . 1. Estimation error for 6-BIOME and PARAFAC-LM using classical and cyclostationary hexacovariance estimators.

Estimation error

SIMULATION RESULTS

We evaluate the performance of 6-BIOME and PARAFAC-LM algorithms when using the cyclostationarity-based hexacovariance estimator proposed in Section 3. The results were obtained from 100 Monte Carlo runs. For each run, the noisy mixture of cyclostationary sources is generated from a simulator of radiocommunication signals developed by A.Chevreuil [START_REF] Chevreuil | Channel simulator[END_REF], which allows the control of the transmission parameters as well as the parameters defining the radio channel. We have considered a reference carrier of 10MHz and a common symbol period of T = 0.4 ms. The sources are modulated using Binary Phase Shift Keying (BPSK). The pulse shaping filter is a square root raised cosine with roll-off factor 0.3 and the sampling rate of the observed data at the receiver is T s = T /2. A uniform linear array of sensors separated by half a wavelength is considered. At each run, the angles of arrival of the sources are randomly drawn between 0 • et 80 • according to a uniform distribution. The Signalto-Noise Ratio (SNR) is fixed at 10 dB in all simulations. At each run, the performance is evaluated from the following normalized error measure: e(H,

H) = H -H 2 F H 2 F .
First, we are interested in evaluating the estimation error assuming 3 sources and 2 sensors. Table 1 shows the median value of the estimation error obtained with 6-BIOME and PARAFAC-LM algorithms using both classical ("class") and cyclostationary ("cyclo") hexacovariance estimators for different data block lengths. It can be seen that the exploitation of higher-order cyclostationarity improves the performance of both algorithms. Note also that the proposed PARAFAC-LM algorithm offers better results than 6-BIOME.

Table 2 shows the Bit-Error-Rate (BER) performance of 6-BIOME and PARAFAC-LM algorithms averaged over 100 Monte Carlo runs.

For each run, the BER is calculated a posteriori using a Maximum-a-Posteriori (MAP) sequence estimator based on the estimated mixing matrix using 5000 observations. As a reference for comparison, the performance of the perfect MAP estimator, which assumes perfect knowledge of the mixing matrix, is shown. For both algorithms, an improved BER performance is observed when higher-order cyclostationarity is exploited. Such a BER improvement is more pronounced in the case (P, M ) = (4, 3) where the performance of the classical hexacovariance estimator is limited due to the higher number of parameters to be estimated.

CONCLUSION

We have addressed the problem of blind identification of underdetermined mixtures in a cyclostationary context by exploiting sixthorder statistics. A "corrected" hexacovariance estimator has been presented which takes into account the second-and fourth-order cyclic moments. Using the proposed estimator, we have assessed the performance of 6-BIOME and PARAFAC-LM algorithms relying on different matrix factorizations of the estimated hexacovariance. Our results show that both algorithms benefit from the exploitation of higher-order cyclostationarity, thus offering an improved identification of the mixing matrix. We have also observed an improved performance of PARAFAC-LM over the 6-BIOME algorithm.

Table 2 .

 2 BER performance of 6-BIOME and PARAFAC-LM obtained with a MAP estimator. SNR=10dB.

	BER
	(P, M )	(3, 2)	(4, 3)
	6-BIOME class.	0.0653 0.1103
	6-BIOME cyclo.	0.0503 0.0771
	PARAFAC-LM class.	0.0315 0.0914
	PARAFAC-LM cyclo.	0.0248 0.0676
	Perfect MAP (H known)	0.0137 0.0064

This work has been supported in part by Amesys contract DP021371 CSE002 "aIntercom".

Note that the LM algorithm has been used in a different context to fit

3rd order tensors in[START_REF] Tomasi | A comparison of algorithms for fitting the parafac model[END_REF].


Laurent Galluccio
email: gallucci@i3s.unice.fr

Olivier Michel
email: olivier.michel@gipsa-lab.inpg.fr

Pierre Comon
email: pcomon@i3s.unice.fr

Arbres de recouvrement minimaux duaux et application à la segmentation non supervisée

Cet article propose de nouvelles approches de segmentation non supervisée. Nous proposons d'exploiter les propriétés d'une nouvelle mesure de distance reposant sur la construction d'arbres de recouvrement minimaux duaux : le « Dual Rooted Prim Tree »(DRooPi), pour construire une paire de classes candidates. Ces classes correspondent à une partition de l'ensemble des sommets connectés par un Droopi. La partition est obtenue par coupure simple du plus grand segment de l'arbre de recouvrement minimal partiel qu'est Droopi, dont les deux racines sont définies par la paire de sommets choisie au hasard. Des fonctions de consensus sont ensuite calculées sur l'ensemble des classes associées à chaque paire de sommets. Pour obtenir la partition finale, un algorithme de classification spectrale est appliqué avec comme mesure de distance les fonctions de consensus définies. Des résultats sont obtenus sur divers données synthétiques et réelles.

Introduction

Une tâche de segmentation non supervisée consiste à définir une partition d'un champ de données multidimensionnel en un ensemble de classes pour lesquelles aucune inofrmation a priori n'est disponible (nombre, forme, tailles respectives...). Nous ne discuterons pas dans ce papier du choix de la métrique utilisée pour définir une mesure de ressemblance entre deux « vecteurs » du champ de données analysé.

Parmi les multiples approches développées pour résoudre ce problème, nous citerons par exemple les méthodes de classification hiérarchique, les algorithmes de partitionnement . . . [START_REF] Jain | Data clustering : a review[END_REF] Les résultats obtenus par différentes méthodes de segmentation, voire par une seule méthode mais avec des paramètres d'initialisation différents, sont en général très divers. Par conséquent, des approches ont été développées afin de combiner un ensemble de partitions « candidates », visant à produire une segmentation finale unique. Certaines de ces approches utilisent des fonctions de consensus basées par exemple sur des matrices de co-association [START_REF] Fred | Combining multiple clusterings using evidence accumulation[END_REF], des votes [START_REF] Ayad | Cumulative voting consensus method for partitions with a variable number of clusters[END_REF], des représentations en hypergraphe [START_REF] Strehl | Cluster ensemble -a knowledge reuse framework for combining multiple partitions[END_REF], des mélanges de modèles [START_REF] Topchy | Clustering ensembles : Models of consensus and weak partitions[END_REF].

Dans cette communication, nous proposons d'utiliser le concept de combinaison de méthodes et une nouvelle méthode de segmentation non supervisée basée sur des arbres à recouvrement minimal (minimal spanning tree : MST) à racines duales. Une caractéristique attractive de la segmentation basée sur les MSTs est sa capacité à prendre en compte une information géométrique aussi bien locale que globale sur la distribution des données à partitionner. Dans [START_REF] Griskschat | Dual rooted-diffusions for clustering and classification on manifolds[END_REF], Griskschat et al. ont montré que les arbres de diffusion à racines duales pour la classification permettent de définir une matrice de similarité qui prend en compte une information de voisinage alors que la distance Euclidienne échoue à capturer cette information. En fait, un arbre à racines duales permet d'extraire des classes à l'intérieur de l'ensemble de données initiales sans définir une partition des données, car tous les sommets n'ont pas besoin d'être connectés. De plus, il présente une variabilité plutôt importante par rapport à l'ensemble d'initialisation des paramètres. Dans ce sens, il peut être considéré comme une méthode de partitionnement faible. Dans notre méthode, la partition finale est obtenue en appliquant un algorithme de classification spectrale [START_REF] Ng | On spectral clustering : Analysis and an algorithm[END_REF] sur des matrices de co-association qui ont montré être des matrices d'affinité admissibles. Toutes ces approches seront abordées dans la Section 2. Dans la section 3, des résultats expérimentaux seront présentés sur des données synthétiques et réelles. Soient X un ensemble de N points ∈ R L et {C 1i , . . . , C Ki } = P i une partition de X en K classes, telle qu'elle est obtenue par A i , où i ∈ [1, . . . M ] ; A i désigne ici un algorithme ou un jeu de paramètres d'initialisation d'un algorithme de segmentation particulier. Soit P = (P 1 , . . . , P M ) l'ensemble des segmentations. Soit G = (V, E) un graphe non dirigé où V = (v 1 , . . . , v N ) est l'ensemble de N sommets et E = (e 1 , . . . , e N -1 ) correspond à l'ensemble de segments. Le poids d'un segment mesure la dissimilarité ou la séparation entre deux sommets. Dans cet article, seule la distance Euclidienne sera considérée bien que l'on peut souligner que d'autres métriques peuvent être utilisées [START_REF] Michel | Unsupervised clustering with mst : Application to asteroïd data[END_REF]. Très souvent en pratique, la méthode utilisée pour construire P est l'algorithme K-means, peu coûteux et facilement implantable. L'inconvénient est que K-means tend à générer des classes convexes, ce qui dans le cas général peut être une limitation. Nous proposons ici de construire l'ensemble P en réalisant une coupure unique (« single cut ») sur le MST obtenu par la réunion de deux sous-MSTs (sMSTs), initialisé en deux sommets distincts et dont la construction est stoppée dès lors que les deux sMSTs possèdent un sommet commun. Plusieurs algorithmes permettent de construire un MST [START_REF] Marchette | Random graphs for statistical pattern recognition[END_REF]. Dans ce papier, la méthode proposée est basée sur les propriétés de construction de l'algorithme de Prim [START_REF] Prim | Shortest connection networks and some generalizations[END_REF]. Le MST est totalement connecté, acyclique (aucune boucle), unique : indépendant du sommet initial (si il n'existe pas d'égalité dans la matrice de distance). L'algorithme de Prim connecte au graphe partiellement connecté à l'itération i le sommet non encore connecté le plus proche du graphe (au sens d'une métrique arbitraire). Les sous-arbres obtenus après N -1 itérations de l'algorithme de Prim sont des MSTs définis sur l'ensemble des N sommets connectés. Le poids d'un MST, défini comme la somme des poids des segments connectant les sommets, est minimal.

Arbre de recouvrement dual

Dans [START_REF] Griskschat | Dual rooted-diffusions for clustering and classification on manifolds[END_REF], Griskschat et al. ont proposé une nouvelle mesure de distance entre deux points (ou sommets) d'un champ de données multidimensionnel, basée sur le temps exprimé en nombre d'itérations de l'algorithme de construction de Prim, nécessaire pour observer la « collision »des deux sMSTs. A chaque itération, les deux sMSTs sont incrémentés1 . Les itérations sont réalisées jusqu'à ce que les deux arbres aient un sommet en commun ; le graphe final n'est pas donc forcément totalement connecté. Tous les points de X ne sont donc pas forcément inclus dans l'un ou l'autre des deux sMSTs quand les itérations sont stoppées. Nous proposons ici une variante de cet algorithme, qui consiste à n'incrémenter à chaque itération qu'un seul des deux sMSTs suivant un critère de coût minimal. Pour une itération donnée, un seul sommet sera connecté à l'un des deux sous-arbres (et non les deux) qui devra être celui de coût minimal. On notera l le coût ou le poids du segment utilisé pour cette connexion. L'intérêt de cette approche réside dans la propriété suivante : permettre une meilleure prise en compte du voisinage de chaque racine dans la définition d'une nouvelle distance. Ce processus de construction s'arrête également quand les deux sous-arbres se rencontrent. N iter désignera le nombre d'itérations effectuées jusqu'à ce que l'algorithme de construction soit stoppé. L'arbre final obtenu en considérant l'ensemble des sommets et segments appartenant à l'un ou l'autre des sMSTs associés à chacune des deux racines sera désigné dans la suite par le nom de « Dual Rooted Prim (Droopi) Tree ». Le graphe Droopi enraciné en x i ∈ X et x j ∈ X sera noté DR(x i , x j ). Différentes mesures d(x i , x j ) de proximité peuvent être proposées à partir du DR(x i , x j ) :

d iter = N iter , d leng = Niter iter=1 l iter , d max = max iter∈[1,Niter] l iter .
Le graphe Droopi satisfait quelques propriétés intéressantes, l'une d'elle sera beaucoup utilisée dans le reste du papier : pour un couple de sommets {v 1 , v 2 } servant comme racines des deux sous-arbres, le dernier segment construit, de poids noté l max -connecte les deux sous-arbres ensemble, -est toujours le plus long (avec un poids maximal) parmi l'ensemble des segments des deux sous-arbres, -définit une distance d max Par conséquent, chaque couple de sommets racines mène à la construction d'un graphe Droopi qui connecte un sousensemble de sommets, qui peut être facilement partitionné en deux classes en appliquant une coupure simple pour enlever le dernier segment. Pour chaque couple de noeuds racines considéré, les membres d'une classe reçoivent le même label et les labels sont enregistrés.

Pour chaque réalisation de l'algorithme de partitionnement de graphe avec différentes initialisations,

P i = {C 1i , C 2i }, i ∈ [1, M ]
représente la partition de l'ensemble de sommets connectés par un couple de sommets {v 1i , v 2i } (où i indique ce choix particulier de racines).

Il est important de souligner que C 1i ∪ C 2i = V et par conséquent, l'algorithme produit deux classes labelisées et un ensemble de sommets non labelisés. Ceci mènera à quelques affinements dans la définition de mesure de co-association, comme décrite dans la sous-section suivante.

Mesures de co-association

Suivant les travaux de Fred et Jain [START_REF] Fred | Combining multiple clusterings using evidence accumulation[END_REF], nous proposons de combiner les résultats de classification (l'ensemble des P i ) afin de calculer une matrice de similarité par une approche d'accumulation de preuves. L'idée directrice en est que des points appartenant à la même classe (inconnue) seront plus souvent affectés à une classe commune C ki par les différentes méthodes de classification ou par les différentes réalisations considérées d'une segmentation partielle (tous les sommets de X ne sont pas pris en compte dans DR(x i , x j )). Par conséquent, prendre en compte la fréquence (mesure de co-association) avec laquelle les points sont labelisés dans une classe commune permet de mettre en évidence certains groupes dominants. Cette mesure a été souvent utilisée dans la littérature [START_REF] Strehl | Cluster ensemble -a knowledge reuse framework for combining multiple partitions[END_REF][START_REF] Fred | Combining multiple clusterings using evidence accumulation[END_REF]. Pour l'ensemble des M partitions réalisées,

co assoc(x i , x j ) = n(x i , x j ) M , (1)
où n(x i , x j ) correspond au nombre de fois qu'une paire de points x i et x j a été classée dans le même groupe parmi les M partitions réalisées. Tous les points ne sont pas labelisés par la réalisation de DR(x i , x j ) (Fig. 1) ; afin de prendre en compte cette information dans la mesure de similarité, nous proposons de définir la mesure de co-association définie co assoc2,

co assoc2(x i , x j ) = n(x i , x j ) m(x i , x j ) , (2)
où m(x i , x j ) est le nombre de partitions dans P où x i ET x j ont été labelisés (0 ≤ m(x i , x j ) ≤ M ).

Bien que ces mesures de co-association reposent sur des heuristiques facilement compréhensibles, elles peuvent être sensiblement modifiées pour mieux correspondre aux caractéristiques du problème considéré. Afin d'aborder le problème de sommets non labelisés existant dans l'approche Droopi, on peut choisir de considérer ces sommets comme étant des éléments d'une classe de rejet. Si nn(x i , x j ) correspond au nombre de fois où aussi bien x i que x j appartiennent à cette classe de rejet, une modification naturelle de la définition (1) pourrait être par conséquent co assoc mod = n+nn M . Comme cette mesure accorde une importance identique aux classes détectées et à la classe de rejet, nous proposons ici de considérer :

co assoc3(x i , x j ) = n(x i , x j ) m(x i , x j ) + nn(x i , x j ) M , (3)
Comme m ≤ M , cela revient à appliquer un plus fort coefficient aux classes détectées et donc un plus faible pour la classe de rejet. Bien que nous ne pouvons pas établir ceci à partir de considérations théoriques, nous observons expérimentalement qu'utiliser cette dernière mesure tend à classer les observations aberrantes ou « outliers »ensemble.

Dans cette étude, nous proposons d'exploiter ces mesures dans le cas particulier d'utilisation d'algorithmes de classification spectrale. Dans ce la matrice d'affinités A considérée est définie par

A(i, j) =∝ β exp + co assoc(x i , x j ) σ (i, j) ∈ [1, N ] 2 , où β et σ sont des constantes à déterminer.
Les résultats présentés sont obtenus par utilisation de l'algorithme de classification spectrale de Ng et al. [START_REF] Ng | On spectral clustering : Analysis and an algorithm[END_REF], qui exploite les propriétes des valeurs propres normalisées du Laplacien normalisé d'un graphe totalement connecté décrit par la matrice A(i, j) 2 . Cette méthode sera nommée par la suite « Evidence Accumulation Clustering by Droopi tree Cut (EAC-DC) ».

Etude expérimentale

Afin de valider et d'évaluer les performances des approches proposées, des tests sur des données synthétiques et des données réelles extraites de la banque de données de « UCI Machine Learning Repository » [START_REF] Asuncion | UCI machine learning repository[END_REF] ont été réalisés. Dans toutes nos expériences, pour chaque algorithme testé, un ensemble de M = 100 paramètres d'initialisation tirés aléatoirement a été considéré.

Nous comparons nos approches avec divers algorithmes basés sur des fonctions de consensus sur l'ensemble des classes :

-classification par accumulation de preuves basée sur une mesure de co-association avec un algorithme hiérarchique « average link »(EAC-AL) [START_REF] Fred | Combining multiple clusterings using evidence accumulation[END_REF], -représentation en graphe et en hypergraphe comme décrite dans [START_REF] Strehl | Cluster ensemble -a knowledge reuse framework for combining multiple partitions[END_REF] (CSPA, HGPA et MCLA), -vote cumulatif, voir [START_REF] Ayad | Cumulative voting consensus method for partitions with a variable number of clusters[END_REF] pour plus de détails et de définitions (URCV, RCV, ACV), -partition médiane (QMI) [START_REF] Topchy | Clustering ensembles : Models of consensus and weak partitions[END_REF]. Ces méthodes utilisent l'algorithme K-means pour former leur ensemble de classes. Il est à noter que, bien que nous considérons le nombre de classes a priori connu, les méthodes citées sont non supervisées.

Résultats sur des données synthétiques

Le tableau 1 reporte les résultats obtenus avec les différentes fonctions de consensus introduites précédemment (Fig. 2 sont représentés les résultats de nos approches sur ces données). La mesure de performance est définie par le rapport entre le nombre d'objets bien classés sur le nombre d'objets total (la référence est connue). Nos approches arrivent à détecter correctement les classes, malgré leurs formes, tandis que les autres méthodes échouent à trouver ces classes non-convexes.

Résultats sur des données réelles

Le premier ensemble de données considéré est « Iris »(150 points dans un espace caractéristique à 4 dimensions) [START_REF] Asuncion | UCI machine learning repository[END_REF]. Il contient trois classes de variétés de fleurs 2. Il faut noter ici qu'en toute généralité, dans le cas où M peut être faible, toutes les distances d(x i , x j ) peuvent n'être pas définies. Dans ce cas, on choisira arbitrairement de fixer A(i, j) = 0.

Figure 1 -

 1 Figure 1 -Arbres de recouvrement minimaux duaux construit sur un ensemble de données. Les symboles X indiquent les sommets racines. Le segment en pointillé correspond au dernier segment connecté.

Figure 2 -

 2 Figure 2 -Résultats obtenus avec nos approches proposées EAC-DC sur différents ensemble de données synthétiques quelle que soit la mesure de co-association utilisée. Les labels des classes sont indiqués par des symboles.

Table 1 -

 1 Résultats obtenus sur les données synthétiques une classe bien séparée des autres et deux classes proches. Le deuxième ensemble de données « breast cancer Wisconsin »a été analysé. Il est constitué de 683 vecteurs de 9 dimensions caractéristiques séparés en deux classes. Le troisième ensemble de données testé est « Wine », composé de 178 vecteurs de 13 dimensions incluant trois classes de cépages différentes. Les résultats obtenus par les différents algorithmes sont présentés dans le tableau 2.On constate que la méthode développée donne de meilleures résultats sur les données Iris et breast cancer. En ce qui concerne les données Wine, les performances de nos approches sont plutôt convaincantes.

	Méthodes	Demi-lunes Deux spirales Deux cercles
	EAC-DC (1, 2, 3)	1	1	1
	EAC-AL	0.7675	0.5550	0.5200
	QMI	0.7550	0.5050	0.5150
	URCV	0.7600	0.5000	0.5200
	RCV	0.7650	0.5300	0.5200
	ACV	0.5000	0.5000	0.5000
	CSPA	0.7300	0.5025	0.5550
	HGPA	0.5000	0.5000	0.5000
	MCLA	0.7275	0.5050	0.5400

avec

Table 2 -

 2 Résultats obtenus sur les données réelles

	Méthode	Iris	Breast Cancer Wisconsin	Wine
	EAC-DC 1	0.9067	0.9678	0.7022
	EAC-DC 2 0.9267	0.9605	0.7022
	EAC-DC 3	0.9200	0.9606	0.7022
	EAC-AL	0.8733	0.9429	0.6910
	CSPA	0.9200	0.8448	0.7135
	HGPA	0.9200	0.6501	0.7022
	MCLA	0.8933	0.9575	0.7247
	URCV	0.8800	0.9590	0.6742
	RCV	0.9000	0.9663	0.6854
	ACV	0.9067	0.9649	0.6966
	QMI	0.8800	0.9356	0.6011

L'incrémentation est ici à comprendre au sens où un nouveau sommet est connecté à l'arbre qui se construit.

ConclusionDans ce papier, nous avons présenté une nouvelle méthode de segmentation non supervisée, basée sur le concept d'accumulation de preuves et sur la combinaison de résultats de multiples segmentations. Nous avons proposé d'exploiter les propriétés d'une nouvelle construction de graphe nommée Droopi pour établir l'ensemble de classes. Cette construction permet de mieux capturer la géométrie intrinsèque globale et locale des données. De multiples réalisations de la segmentation par une coupure sur le MST sont effectués à partir d'un ensemble de paramètres d'initialisation choisis aléatoirement et permettent de construire un ensemble de classes grâce à de nouvelles mesures de co-association proposées. L'extension de la méthode proposée avec des arbres à K racines est en cours d'étude. Les performances de cette nouvelle méthode ont été évaluées sur un ensemble de données synthétiques et réelles, soulignant des caractéristiques des algorithmes proposés très intéressantes et prometteuses.


L Albera

A Kachenoura

A Karfoul

P Comon

L Senhadji

One decade of biomedical problems using ICA: a full comparative study

Keywords: Independent Component Analysis (ICA), biomedical applications, comparative study

This communication aims at giving some insights into the use of Independent Component Analysis (ICA) for solving biomedical problems. First the concept of ICA is reviewed and different classes of ICA methods are described. Next a survey on most encountered biomedical problems solved using ICA is detailed. Finally a comparative performance study of thirteen ICA algorithms is performed on biomedical signals.

I. WHAT IS ICA?

Given one realization of a N-dimensional random vector process whose m-th vector is defined by:

x

[m] = As[m] + n[m] ( 1
)
where s[m] is a P-dimensional source vector with statistically independent components, ν[m] is a noise vector which is independent from the source vector and A is the so-called mixing matrix. The ICA problem consists in finding, only from the data, a (N×P) matrix W , called separator, such that:

y[m] = W T x[m] ( 2
)
is an estimate of s[m] to within a diagonal matrix and a permutation. But how to use the statistical independence of the source components in order to restore them?

Entropy and mutual information. A natural way to measure statistical dependence consists in using the Mutual Information (MI) [START_REF] Haykin | Source separation: Models, concepts, algorithms and performance in Unsupervised Adaptive Filtering[END_REF], wich can be easily related to the (Shannon) Differential Entropy (DE) or its normalized version called negentropy. The INFOMAX method solves the ICA problem by maximizing an output DE using the natural gradient algorithm [START_REF] Haykin | Source separation: Models, concepts, algorithms and performance in Unsupervised Adaptive Filtering[END_REF]. The PICA [START_REF] Karvanen | Blind separation methods based on pearson system and its extensions[END_REF] algorithm uses the parametric Pearson model in the minimization of the MI. The FastICA technique extracts each source by maximizing an approximation of the negentropy by means of an approximate Newton iteration. Another way to measure statistical dependence, less natural but easier to compute consists in using cumulants.

Cumulants. The r-th order cumulants are the coefficients of the Taylor expansion of the second characteristic function about the origin. They enjoy very attractive properties such that i) if at least two components or groups of components are statistically independent, then all cumulants involving these components are null, ii) if a variable is Gaussian, then its HO marginal cumulant is null, iii) cumulants are symmetric arrays, say the value of their entries does not change by permutation of their indices and iv) cumulants satisfy the multi-linearity property. The SOBI [START_REF] Haykin | Source separation: Models, concepts, algorithms and performance in Unsupervised Adaptive Filtering[END_REF], TFBSS [START_REF] Fevotte | Two contributions to blind sources separation using time-frequency distributions[END_REF] and SOBIUM [START_REF] Lathauwer | Blind identification of underdetermined mixtures by simultaneous matrix diagonalization[END_REF] techniques jointly diagonalize a set of Second Order (SO) cumulant matrices built from the data while COM2 [START_REF] Haykin | Source separation: Models, concepts, algorithms and performance in Unsupervised Adaptive Filtering[END_REF], STOTD [START_REF] Lathauwer | Independent component analysis and (simultaneous) third-order tensor diagonalisation[END_REF] and JADE [START_REF] Haykin | Source separation: Models, concepts, algorithms and performance in Unsupervised Adaptive Filtering[END_REF] diagonalize in different ways a Fourth Order (FO) cumulant matrix computed from the data. The ALS-FUBI ACDC algorithm diagonalises a set of FO cumulant matrices calculated from the data, as FO-BIUM, but using a non-orthogonal diagonalization scheme. Eventually, the ELScaf [START_REF] Comon | Blind identification of underdetermined mixtures based on the characteristic function[END_REF] use pseudocumulants of the data defined as the derivatives of the second characteristic function in points different form the origin.

II. A BIBLIOGRAPHICAL SURVEY ON BIOMEDICAL APPLICATIONS

In the last decade, ICA-based methods have been widely used in the field of biomedical engineering, especially to analyze electrophysiological signals.

A. Functional brain imaging

Promising results have been reported in electroencephalographical signal processing using ICA techniques. They include Evoked Potentials (EP) enhancement, categorized brain signals detection, spindles detection and estimation, and artifacts re-duction. The analysis of two particular neurophysiological signals is investigated in this context. P-300 evoked potentials. P-300 is a positive ERP, which occurs over the parietal cortex with a latency of about 300 ms after rare or task-relevant stimuli. The P-300 can be obtained in all stimulus modalities (auditory as well as visual and somatosensory modalities) and can even be produced by the omission of a stimulus in a regular train of stimuli [START_REF] Niedermeyer | Electroencephalography: Basic Principles, Clinical Applications, and Related Fields[END_REF]. Due to the poor SNR as well as to the presence of artifacts (such as ocular, muscular and cardiac activities), the P-300 wave can be buried in the signal collection. Hence, the main objective when applying ICA to P-300-based BCI systems is two-fold: i) to denoise the EEG signal in order to enhance the SNR of the P-300 and ii) to separate ERP responses to target and non-target ones. The first point was considered by Bayliss and al. in [START_REF] Bayliss | Single trial p300 recognition in a virtual environnement[END_REF]. Authors described an experiment demonstrating the existence of a P-300 wave when facing red stoplights and the absence of this signal when facing yellow stoplights in a virtual driving environment. Bayliss et al. pointed out that most of artifacts were due to eye movements. They showed that an ICA technique was able to separate the background EEG signal and eye movements from the P-300 signal. In [START_REF] Xu | BCI competition 2003-data set iib: Enhancing p300 wave detection using ica-based subspace projections for BCI applications[END_REF], Xu et al. dealt with the second point and proposed to enhance the P-300 wave detection in the P-300 speller paradigm used to record the database IIb of BCI Competition 2003 [START_REF] Blankertz | The BCI competition 2003: progress and perspectives in detection and discrimination of EEG single trials[END_REF]. As classical methods for enhancing the detection of P-300 components are time consuming, authors proposed to use an ICA technique in the training phase. The key issue of this study is the selection of the meaningful Independent Components (ICs). Indeed, according to the prior physiological knowledge, authors proposed two additional post-processing steps, namely the temporal manipulation of ICs and the spatial manipulation of ICs (see [14, section II] for details). They showed that the proposed algorithm for P-300 detection based on ICA provided a perfect accuracy (100%) in the competition.

Mu rhythm and other activities from sensorimotor cortex. EEG contains a fairly wide frequency spectrum. Nevertheless, the relevant frequency range from the psychophysiological viewpoint lies between 0.1 Hz and 100 Hz [START_REF] Niedermeyer | Electroencephalography: Basic Principles, Clinical Applications, and Related Fields[END_REF]. For example, Beta rhythm is associated with active thinking and attention whereas the Mu rhythm is affected by movements or movement imagery. QIN et al. [START_REF] Qin | Motor imagery classification by means of source analysis for brain-computer interface applications[END_REF] present a pilot study aimed at classifying motor imagery, using ICA as a spatio-temporal filter. More precisely, the study was focused on the Mu rhythm which decreases or desynchronizes with movement or movement imagery. Authors conduct an experi-ment where the subject was asked to imagine right or left-hand movement according to a timetable. EEG were recorded using 59 electrodes but only those located around sensorimotor cortex were used in the study. ICA was used to extract ICs related to the left and the right motor imagery task. The obtained results show that ICA plays an important role in extracting a useful feature that identifies the imagined hand movement. A promising classification rate (about 80 %) of left or right-hand movement imagery was obtained on human subject studies, based only on a single trial and without any training procedure.

B. Electrocardiogram signal analysis

The electrocardiogram (ECG) reflects the electrical activity of the heart which is usually recorded with surface electrodes placed on the chest, arms and limbs. Although the ECG recording techniques are very effective, the distortions caused by noises and interfering physiological signals are still very important. Let us consider in more details the two following classical applications.

Noninvasive fetal ECG extraction. The FECG contains a lot of useful information which can provide valuable clinical indication about the fetal well-being and then allow for an early diagnosis of fetal cardiac abnormalities and other pathologies. One of the challenging tasks is to reliably detect and enhance FECG in a non-invasive fashion by using several skin electrodes located on a pregnant woman's body. In fact, the obtained ECG also contains the Mother ECG (MECG), and other interfering signals including the respiration, the Elec-troMyoGram (EMG), the electrodes movements, etc. De Lathauwer et al. [START_REF] De Lathauwer | Fetal electrocardiogram extraction by blind source subspace separation[END_REF] show that the ECG data recorded from mother's skin electrodes can be modeled as a linear static mixture of independent sources. Indeed, The transfer function between the bioelectric sources to body surface electrodes is assumed to be linear, resistive and the hight propagation velocity of the electrical signal in the human tissues validates the instantaneous assumption of the model. Besides, the authors state in accordance with that the MECG-subspace is characterized by a three-dimensional vector signal, whereas the dimension of FECG-subspace is subject to changes during the pregnancy period. Applying the COM2 ICA method to eight mother surface lead ECG, the authors show that the full three-dimension MECG subspace is well reconstructed whereas the FECG subspace is extracted in two different components.

Atrial activity extraction for atrial fibrillation. Atrial Fibrillation (AF) is associated with in-creased mortality and hospitalization for most people, and understanding of pathological mechanisms underlying AF using non-invasive diagnosis tools such as surface ECG is crucial to improve the patients treatment strategies. However, due to the low Signal to Noise Ratio (SNR) of the Atrial Activity (AA) on the surface 12-lead ECG, the analysis of AF episode remains difficult. The AF problem was tackled by means of ICA methods in [START_REF] Rieta | Atrial activity extraction for atrial fibrillation analysis using blind source separation[END_REF]. The authors justify, on the one hand, three basic considerations about the AA and the Ventricular Activity (VA) and, on the other hand, the fashion which both activities are recorded from the surface electrodes: i) the independence between VA and AA, ii) the sub-Gaussian and super-Gaussian distributions of AA and VA, respectively and iii) the linear static mixing model followed by the data. The authors applied the FastICA algorithm on the 12-lead ECG of seven patients suffering from AF. The obtained results show that three separated sources have a more sub-Gaussian distribution and hence are candidates to be related to AA. Four sources are associated with Gaussian noise and artifact, whereas five extracted components which present a super-Gaussian distribution contain a VA. One important problem that arises when ICA is used in biomedical context is to automatically select and classify independent source of interest. Typically, in the above example the question is: how to choose the most informative AA source among the three extracted components with a sub-Gaussian distribution? The authors solve this problem by exploiting some prior information about spectral content of AA during AF episode. Indeed, the AA signal exhibits a narrowband spectrum with a main frequency of between 3.5 -9 Hz. Thus, they apply a spectral analysis over all the sources with sub-Gaussian distribution and choose the component which presents a major peak at frequency f = 6.31 Hz, as the AA source.

III. PERFORMANCE COMPARISON OF THIRTEEN ICA METHODS

Performance criterion.

In order to quantitatively compare several ICA methods, we used the Signal to Interference-plus-Noise Ratio (SINR) of each source after separation as performance criterion [START_REF] Haykin | Source separation: Models, concepts, algorithms and performance in Unsupervised Adaptive Filtering[END_REF]. More precisely, the SINR of the p-th source at the i-th output of the separator W = [w 1 , . . . , w P ] is defined by:

SINR p [w i ] = π p |w i T a p | 2 /(w i T R ν p w i ) (3)
where π p represents the power of the p-th source, w i the i-th column of the separator W and R ν p is the total noise covariance matrix for the p-th source, corresponding to the estimated data covariance matrix in the absence of the component p. On the basis of these definitions, the restitution quality of the p-th source at the output of the separator W is evaluated by computing the maximum of SINR p [w i ] with respect to i for 1 ≤ i ≤ P. This quantity is denoted by SINRM p The performance of a source separator W is defined by the line vector SINRM(W ) = (SINRM 1 [W ], . . . , SINRM P [W ]). In a given context, a separator W (1) is better than another one W (2) for retrieving the source p, provided that SINRM p [W (1) ] > SINRM p [W (2) ].

Data generation. We consider the case of the extraction of the Mu rhythm presented in section A. when seven surface electrodes, located around the sensorimotor cortex, are used to record EEG data. In such a context the surface observations can be considered as a noisy mixture of one source of interest, namely the Mu signal, and artifact sources such as an ocular. The intracerebral Mu wave located in the motor cortex is simulated using the parametric model of Jansen [START_REF] Jansen | Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns?[END_REF] whose parameters are selected to derive a Mu-like activity. The ocular signal is issued from our polysommnographic database. As far as the additive noise is considered, it is modeled as the sum of instrumental and physiological noises. Therefore, a Gaussian vector process is used to simulate the instrumental noise while a brain volume conduction of 800 independent EEG sources, generated using the Jansen model, is simulated to produce a surface background EEG activity. Finally the mixing matrix is defined as the concatenation of two columns modeling the head volume conduction of the Mu rythm and the ocular activity toward the surface.

Computer results. The SINRM of each source at the output of SOBI, SOBI rob (a SOBI modified version based on a robust whitening procedure), SOBIUM, TFBSS, PICA, INFOMAX, STOTD, COM2, JADE, FastICA (FastICA def and FastICA sym , the deflation and symmetric versions, respectively), ALS-FUBI ACDC and ELScaf was computed as a function of the number of samples (with a sampling rate of 256 Hz) and for a SNR of 5dB. Considering the ocular activity's extraction (i.e. SINRM 1 in the left column), we note a higher performance of SOBI rob , ALS-FUBI ACDC and SOBIUM (i.e. -0.5dB with respect to the optimal separator (SMF)). Regarding the other algorithms, their performance is still reasonable (i.e. -2dB from SMF) except for INFOMAX and PICA which separators seem to be biased in this context. As far as the Mu wave's extraction is concerned (i.e. SINRM 2 in the right column), we observe the very good behavior of the SOBI rob method with respect the other ones. This is due to the fact that SOBI rob is not affected by the spatially correlated noise. The ALScaf algorithm fails in the extraction of the Mu wave whereas ALS-FUBI ACDC requires a large amount of samples to converge. The quasi-Gaussian distribution of the Mu wave may explain such a behavior. In the case of the ALScaf method, the use of more grid points or the exploitation of several HO derivatives should improve its performance. The STOTD, SOBI, COM2, JADE, TFBSS, SOBIUM, INFOMAX, FastICA sym and FastICA def globally give acceptable results (i.e. -1.5dB from the SMF). Note that the PICA algorithm provides poor separation of the Mu source.

IV. CONCLUSION

The concept of ICA is briefly reviewed in this paper allowing to cover a wide range of techniques. In addition, ICA for biomedical engineering is tackled by giving a survey on its widespread applications. A comparative performance analysis on electrophysiological data reproducing the real scalp EEG recordings is conducted. The obtained results show that the performance of the ICA algorithms depend on the electrophysiological nature of the sources that is extracted. Hence, the collection and the exploitation of statistical and physiological information on the sources of interest, such as temporal color or distribution (Gaussian or not), can help to choose the more appropriate ICA method. In our study context, the SOBI rob method gives the best results probably due to the source temporal nonwhiteness and the noise temporal whiteness.

Fig. 1 :

 1 Fig. 1: Criterion SINRM as a function of data samples for N = 7 surface electrodes and P = 2 sources and with a SNR of 5dB.

ACKNOWLEDGEMENTS

This work is supported in part by the French Government under two ANR Contracts, namely DECOTES (06-BLAN-0074) and mv-EMD (BLAN07-0314-02).