Chloé Dimeglio

Jean-Michel Loubes

Elie Maza

MANIFOLD EMBEDDING FOR CURVE REGISTRATION

Keywords: AMS 2000 subject classifications: Primary 62G05; secondary 62M99 Manifold learning, Intrinsic statistics, Structural statistics, Graph-based methods, Curve alignment, Curve registration, Warping Model, Functional data

We focus on the problem of finding a good representative of a sample of random curves warped from a common pattern f . We first prove that such a problem can be moved onto a manifold framework. Then, we propose an estimation of the common pattern f based on an approximated geodesic distance on a suitable manifold. We then compare the proposed method to more classical methods.

1. Introduction. The outcome of a statistical process is often a sample of curves {f i , i = 1, . . . , m} showing an unknown common structural pattern, f , which characterizes the behaviour of the observations. Examples are numerous, among others, growth curves analysis in biology and medicine, quantitative analysis of microarrays in molecular biology and genetics, speech signals recognition in engineering, study of expenditure and income curves in economics. . . . Hence, among the last decades, there has been a growing interest to develop statistical methodologies which enables to recover from the observation functions a single "mean curve" that conveys all the information of the data.

A major difficulty comes from the fact that there are both amplitude variation (in the y-axis) or phase variation (in the x-axis) which prevent any direction extraction of the mean, median, correlations or any other statistical indices for a standard multivariate procedure such as principal component analysis, and canonical correlations analysis, see [START_REF] Kneip | Statistical tools to analyze data representing a sample of curves[END_REF] or [START_REF] Ramsay | Functional data analysis[END_REF] and references therein. Indeed, the classical cross-sectional mean does not provide a consistent estimate of the function of interest f since it fails to capture the structural characteristics in the sample of curves as quoted in [START_REF] Ramsay | Curve registration[END_REF]. Hence, curve registration methods (also called curve alignment, structural averaging, or time warping) have been proposed in the statistical literature. We refer to, just to name a few, [START_REF] Sakoe | Dynamic programming algorithm optimization for spoken word recognition[END_REF] in Spoken Word Recognition domain, [START_REF] Kneip | Statistical tools to analyze data representing a sample of curves[END_REF] for Landmark Registration, [START_REF] Silverman | Incorporating parametric effects into functional principal components analysis[END_REF] for a functional principal component analysis, [START_REF] Wang | Alignment of curves by dynamic time warping[END_REF] for Dynamic Time Warping, [START_REF] Ramsay | Curve registration[END_REF] for Continuous Monotone Registration, Rønn [2001] for shifted curves, [START_REF] Liu | Functional convex averaging and synchronization for time-warped random curves[END_REF] for functional convex averaging, [START_REF] Gervini | Nonparametric maximum likelihood estimation of the structural mean of a sample of curves[END_REF] for maximum likelihood estimation, [START_REF] Gamboa | Semi-parametric estimation of shifts[END_REF] for shifts estimation, James [2007] for alignment by moments, and [START_REF] Dupuy | Non parametric estimation of the structural expectation of a stochastic increasing function[END_REF] for Structural Expectation estimation.

This issue is closely related to the problem of finding the mean of observations lying in a space with an unknown, non necessarily euclidean, underlying geometry. The problem is thus twofold.

First, the mere definition of the mean should be carefully studied. Indeed, let E = {X 1 , . . . , X n } be a sample of i.i.d random variables of law X ∈ M where M is a submanifold of R p . If we denote by d the Euclidean distance on R p , then the classical sample mean, or Fréchet sample mean, defined by (1) µ = arg min

µ∈R p n i=1 d 2 (X i , µ)
is not always a good representative of the given sample E, and, obviously, of the underlying population. Using the geometry of the manifold, it seems natural to replace Criterion (1) by

µ I = arg min µ∈M n i=1 δ 2 (X i , µ)
where δ is the geodesic distance on manifold M, giving rise to the intrinsic mean, whose existence and properties are studied, for instance, in [START_REF] Bhattacharya | Large sample theory of intrinsic and extrinsic sample means on manifolds[END_REF]. When dealing with functional data, we assume that the functions f i can be modeled as variables with values on a manifold, and curve registration amounts to considering an intrinsic statistic that reflects the behaviour of the data. In the following we will consider, for α > 0,

(2) In particular, for α = 1, we will deal with µ 1 I , the intrinsic sample median. Second, previous construction relies on the choice of the embedding which may not be unique, then the manifold itself and its underlying geodesic distance. Actually we only have at hand a sample of random variables which are sought to be a discretization of an unobserved manifold. Over the last decade, some new technics have been developed to find and compute the natural embedding of data onto a manifold and to estimate the corresponding geodesic distance, see for instance de Silva and [START_REF] De | Unsupervised learning of curved manifolds[END_REF] for a review of global (Isomap type) and local (LLE type) procedures, while applications have been widely developed, see for instance [START_REF] Pennec | Intrinsic statistics on riemannian manifolds: basic tools for geometric measurements[END_REF].

µ α I = arg min µ∈M n i=1 δ α (X i , µ
In the following, we will consider an approximation, achieved with a graph theory approach inspired by works on manifold learning and dimension reduction [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF]]. We will first show that curve registration for parametric transformations can be solved using a manifold geodesic approximation procedure. Then, we will highlight that this enables to recover a mean pattern which conveys the information of a group of curves. This pattern is used for curve classification for simulated data and real data which consists in predicting a particular landscape using the reflectance of the vegetation.

This article falls into the following parts. Section 2 is devoted to the construction of the approximated geodesic distance. In Section 3, we describe the manifold framework point of view for curve registration. We then explain how to estimate a representative of a sample of warped curves. The performance of this estimator is then studied in Section 4 using simulated data, and in Section 5 with a real data set. Concluding remarks are given in Section 6. Proofs are gathered in Section 7.

2. A graph construction for topology estimation. Let X be a random variable with values in an unknown connected and geodesically complete Riemannian manifold M ⊂ R p . We observe an i.i.d sample E = {X i ∈ M, i = 1, . . . , n} with distribution X. Set d the Euclidean distance on R p and δ the induced geodesic distance on M. Our aim is to estimate intrinsic statistics defined by Equation (2). Since the manifold M is unknown, the main issue is to estimate the geodesic distance between two points on the manifold, that is δ (X i , X j ).

Let γ ij be the geodesic path connecting two points X i and X j , that is the minimum length path on M between points X i and X j . Denoting by L (γ) the length of a given path γ on M, we have that δ (X i , X j ) = L (γ ij ).

In the Isomap algorithm, [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF] propose to learn manifold topology from a graph connecting k-nearest neighbors for a given integer k. In the same way, our purpose is to approximate the geodesic distance δ with a suitable graph connecting nearest neighbors. Our approximation is carried out in three steps. Thereafter, we denote g ij a path connecting two points X i and X j on a given graph, and L (g ij ) the length of such a path.

Step 1. Consider K = (E, E) the complete Euclidean graph associated to sample E, that is the graph made with all the points of the sample E as vertices, and with edges

E = {{X i , X j } , i = 1, . . . , n -1, j = i + 1, . . . , n} .
For an Euclidean graph, the edge weights are the edge lengths, that is, the Euclidean distances between each pair of points.

Step 2. Let T = (E, E T ) be the Euclidean Minimum Spanning Tree (EMST) associated to K, that is, the spanning tree that minimizes

{X i ,X j }∈E T d (X i , X j ) .
The underlying idea in this construction is that, if two points X i and X j are relatively close, then we have that δ (X i , X j ) ≈ d (X i , X j ). This may not be true if the manifold is very twisted and if too few points are observed, and may induce bad approximations, hence the algorithm will produce a good approximation for relatively regular manifolds. It also generally requires a large number of sampling points on the manifold in order to guarantee the quality of this approximation. This drawback is well known when dealing with graph based approximation of the geodesic distance. Then, the graph T is a connected graph spanning K which mimics the manifold M. Furthermore, an approximation of the geodesic distance δ (X i , X j ) is provided by the sum of all the euclidean distance of the edges of the shortest path on T connecting X i to X j , namely δ (X i , X j ) = min

g ij ∈T L (g ij ) .
However, this approximation is too sensitive to perturbations of the data, and hence, very unstable. To cope with this problem, we propose to add more edges between the data to add extra paths in the data sample and thus to increase stability of the estimator. The idea is that paths which are close to the ones selected in the construction of the EMST could provide alternate ways of connecting the edges. Close should be here understood as lying in balls around the observed points. Hence, these new paths between the data are admissible and should be added to the edges of the graph. This imsart-aoas ver. 2009/02/27 file: ManiWarp.hyper9470.tex date: May 27, 2011 provides redundant information but also stabilizes the constructed distance, and may also provide an answer to the the main defect of the algorithm that considers that two points very close with respect to the Euclidean distance are also close with respect to the geodesic distance.

Step 3. Let B (X i , ǫ i ) ⊂ R p the open ball of center X i and radius ǫ i defined by

ǫ i = max {X i ,X j }∈E T d (X i , X j ) . Let graph K ′ = (E, E ′ ) defined by {X i , X j } ∈ E ′ ⇐⇒ X i X j ⊂ n i=1 B (X i , ǫ i )
where

X i X j = {X ∈ R p , ∃λ ∈ [0, 1], X = λX j + (1 -λ)X i } .
Then, K ′ is the graph which gives rise to our estimator of the distance δ :

(

3) δ (X i , X j ) = min g ij ∈K ′ L (g ij ) .
Hence, δ is the distance associated with K ′ , that is, for each pair of points X i and X j , we have δ (X i , X j ) = L (γ ij ) where γij is the minimum length path between X i and X j associated to K ′ . We note that, the 3-steps procedure described above contains widespread graph-based methods to achieve our purpose. In this article, our graph-based calculations, such as MST estimation or shortest path calculus, were carried out with the R Language [R Development Core Team, 2010] with the igraph package for network analysis [START_REF] Csardi | The igraph software package for complex network research[END_REF].

An example of this 3-steps procedure and its behaviour when the number of observations increases are displayed respectively in Figure 1 and Figure 2. In Figure 1, points X 1 i , X 2 i i are simulated as follows :

(4)

X 1 i = 2i -n -1 n -1 + ǫ 1 i and X 2 i = 2 2i -n -1 n -1 2 + ǫ 2 i
where ǫ 1 i and ǫ 2 i are normaly distributed with mean 0 and variance 0.01. In Figure 2, we give some results of graph K ′ for n ∈ {10, 30, 100, 300}. We can see in such a figure that graph K ′ tends to be close to the manifold t, t 2 ∈ R 2 , t ∈ R . The main difference between our algorithm and the Isomap algorithm lies in the treatment of points which are far from the others. Indeed, the first

Data

Step 1

Step 2

Step 3 step of the original Isomap algorithm consists in constructing the k-nearest neighbor graph or the ǫ-nearest neighbor graph for a given integer k or a real ǫ > 0. Hence, points which are not connected to the biggest graph, since they are too distant, are not used for the construction of the estimated distance. Such a step is not present in our algorithm since in the applications we consider a distant point is not always an outlier. Hence, we do not exclude any points, and rather, for the construction of the EMST, all points of the data set are connected. Moreover, the Isomap algorithm requires the choice of parameters which are closely related to the local curvature of the manifold (see, for instance, [START_REF] Balasubramanian | The isomap algorithm and topological stability[END_REF]). This involves a heavy computing phase which is crucial for the quality of the construction, while, in our version we tend to give an automatic selection of parameters. We will show in Section 3 that both procedures used for curve registration behave in a similar way and over performs other standard feature extraction methods. In the following section, we present a new application of manifold learning to the curve alignment problem.

3. Application to curve alignment. Consider a function f : R → R, which will be the pattern to be recovered, observed in a translation effect framework. Let A be a real valued random variable with unknown distribution on an interval (b, c) ⊂ R. The observation model is defined by (5)

X j i = f (t j -A i ) , i ∈ {1, .
. . , n}, j ∈ {1, . . . , m}, where (A i ) i are i.i.d random variables drawn with distribution A which model the unknown translation parameters, while (t j ) j ∈ R m stand for the measurement points.

This situation usually happens when individuals experience similar events, which are explained by a common pattern f , and when the starting times of the events are not synchronized. Such a model has been studied, for instance, in [START_REF] Silverman | Incorporating parametric effects into functional principal components analysis[END_REF][START_REF] Birgitte | Nonparametric maximum likelihood estimation for shifted curves[END_REF]. This issue has also received a specific attention in a semi-parametric framework in [START_REF] Gamboa | Semi-parametric estimation of shifts[END_REF] or [START_REF] Castillo | Estimation of the distribution of random shifts deformation[END_REF]. In these works, among others, shift parameters are estimated, which enables to align the observations and thus to get rid of the translation issue. Model (5) falls also under the generic warping model proposed in [START_REF] Maza | Estimation de l'espérance structurelle d'une fonction aléatoire[END_REF] and in [START_REF] Dupuy | Non parametric estimation of the structural expectation of a stochastic increasing function[END_REF] which purpose is to estimate the underlying structure of the curves. For this, the authors define the structural median f SM of the data. In the case of translation effects, it corresponds to ( 6)

f SM = f (• -med(A))
with med(A) the median of A. Hence, a natural estimator of the structural median f SM , related to Model (5), would be

(7) f SM = f t 1 -med(A) , f t 2 -med(A) , . . . , f t m -med(A)
with med(A) the median of sample (A i ) i . However, we first note that the translation parameters (A i ) i are not observed, and, as a consequence, that the median med(A) can not directly be calculated. Then, the function f is also unknown, so, estimating med(A) is not enough to calculate f SM . Our purpose is to show that our manifold point of view provides a direct estimate of f SM without the prior estimation of med(A).

In order to use the manifold embedding approach, define

X : R → R m a → X(a) = (f (t 1 -a) , f (t 2 -a) , . . . , f (t m -a))
and set C = {X(a) ∈ R m , a ∈ R} .

As soon as f ′ = 0, the map X : a → X(a) provides a natural parametrization of C which can thus be seen as a submanifold of R m of dimension 1. The corresponding geodesic distance is given by

δ (X(a 1 ), X(a 2 )) = a 2 a 1 X ′ (a) da .
The observation model ( 5) can be seen as a discretization of the manifold C for different values (A i ) i . Finding the median of all the shifted curves can hence be done by understanding the geometry of space C, and thus approximating the geodesic distance between the curves.

The following theorem states that the structural median f SM defined by Equation ( 7) is equivalent to the median with respect to the geodesic distance on C, that is

µ 1 I = arg min µ∈C n i=1 δ (X i , µ) ,
which provides a geometrical interpretation of the structural median.

imsart-aoas ver. 2009/02/27 file: ManiWarp.hyper9470.tex date: May 27, 2011

Theorem 1. Under the assumption that f ′ = 0, we get that

µ 1 I = f SM .
Previous theorem can be extended to the more complex case of parametric deformations of the type

X : R 3 → R m (a, b, c) → X(a, b, c) = (af (t 1 -b) + c, . . . , af (t m -b) + c)
as soon as a = 0 and f ′ = 0. Such a model has been described in [START_REF] Vimond | Efficient estimation for a subclass of shape invariant models[END_REF] and in [START_REF] Bigot | Semiparametric estimation of shifts on compact lie groups for image registration[END_REF]. In this case, the submanifold is obviously of dimension 3.

In an even more general framework, when the observations can be modeled by a set of curves warped one from another by an unobservable deformation process, this estimate enables to recover the main pattern. It relies on the assumption that all the data belong to a manifold whose geodesic distance can be well approximated by the graph structure of the modified minimal spanning tree described in Section 2.

Finally, we propose the following estimator of the structural median δ (X i , µ) , using the geodesic distance δ, estimated by the algorithm described in Section 2.

The numerical properties of this estimator is studied using simulations in Section 4, and for real data sets in Section 5.

Simulations.

We consider the target function f : R → R defined by f (t) = t sin(t). We simulate deformations of this function on j = 1, . . . , m = 100 equally distributed points t j of the interval [-10, 10], according to the following model :

(9) Y i (t j ) = A i f (B i t j -C i ) , i =, . . . , n, j = 1, . . . , m,
where (A i ) i and (C i ) i are i.i.d uniform random variables on [-10, 10] while (B i ) i is an i.i.d sample of a uniform distribution on [-1, 1]. We finally obtain a data set of n = 100 curves where each differs from the initial function f by a translation and an amplitude deformation. The data is displayed on the left graph of Figure 3.

We then consider four estimators of the function f . The first one, which minimizes the approximated geodesic distance, defined by Equation (8), will be referred to as the structural median estimator. The second one is obtained by the Curve Alignment by Moments procedure (CAM) developped by James [2007]. The third one is the template obtained with the Isomap strategy, with the "isomap" function of the R package vegan [START_REF] Jari Oksanen | Community Ecology Package[END_REF]. The last one is the mere mean of the data.

We recall here that the CAM procedure consists on extracting the mean pattern by synchronization of the moments of the simulated curves. For this, James [2007] introduces the feature function concept for a given function g, defined as I g (t) :

I g (t) ≥ 0 and I g (t)dt = 1
and the moments

µ (1) g = tI g (t)dt and µ (k) g = t -µ (1) g k I g (t)dt, k ≥ 2.
Then, the CAM procedure align the curves by warping their moments, for instance, the amplitude of the peaks, at the location they occur, the variance around these peaks, and so on. This method relies on the choice of a proper feature function, for instance I (l)

g (t) = |g (l) (t)|/ |g (l) (s)
|ds for a given l ≥ 0, on an approximation of the functions by splines, and the selection of the number of moments to be synchronized. Hence, it highly depends on the choice of these tuning parameters. We have chosen the optimal value of the parameters over a grid.

These four estimators are shown on Figure 3. With the CAM or the mere mean procedure, the average curve does not reflect the structure of the initial curves, or the amplitude of their variations. On the contrary, the structural median extracted by Manifold Warping has the characteristics of the closest target curve, but is also its best approximation together with the pattern obtained with the Isomap strategy. Note here that our version of the algorithm for curves provide the same kind of template and is parameter free while parameters governing the dimension of the manifold embedding must be chosen for the Isomap procedure. Nevertheless, both procedures are competitive and lead to similar performance. and 5, we present two data sets corresponding to reflectance patterns of two landscapes in the same region with the same period. However, the reflectance depends on the vegetation whose growth depends on the weather condition and the behavior in soil. It is therefore relevant to consider that these profiles are deformations in translation and/or amplitude of a single representative function of the reflectance behaviour of each landscape in this region at this time.

Our aim is to build a classification procedure. For this, we will use a labeled set of curves and extract from each group of similar landscape a representative profile. Then, we will allocate a new curve to the group whose representative curve will be the closest. That is the reason why it is important to obtain a pattern which captures the structure of the curves. We will use three different ways to get a representative group of curves, the mean curve, the CAM method and our method, referred to as the Manifold Warping. We will compare their classification performance together with a usual classification procedure : the classical k-nearest neighbours.

In Figure 4, we observe that the CAM average oversmoothes the peaks of activity at times 12 and 22 to make them almost nonexistent. This is a clear defect since, according to the experts of landscape remote sensing, these peaks of activity are representative of the nature of landscape. Indeed, these peaks convey essential informations which determines, among other things, the type of landscape. On the other hand, these changes are very well rendered by the pattern obtained by Manifold Warping. The same conclusions can be drawn in Figure 5 for an other landscape. In this application domain, extracting a curve by Manifold Warping is best able to report data as reflecting their structure and thus to obtain a better representative. Now, we try to identify "unknown" landscapes by comparing each curve to the mean pattern of each group. The allocation rule is built using the Euclidean distance. Note that here we have sought to classify the landscapes, not using the whole curve which would correspond to a whole year of observation but using only a part of the curves, namely all the observations before t = 30. To benchmark our procedure, we compare our performance to the method of the k-nearest neighbors classification. Finally, we obtain the confusion matrices displayed in Tables 1 and2. We get a much better discrimination of landscapes with the method consisting in estimating a representative by Manifold Warping than by the CAM method or by classical mean.

Conclusion.

By using an Isomap inspired strategy, we have extracted from a pattern of curves, a curve which, playing the role of the mean, serves as a pattern conveying the information of the data. In some cases, in particular when the structure of the deformations entails that the curve can be embedded into a manifold regular enough, we have shown that this corresponds to finding the structural expectation of the data, devel- oped in [START_REF] Dupuy | Non parametric estimation of the structural expectation of a stochastic increasing function[END_REF], which improves the performance of other mean extraction methods. This enables to derive a classification strategy that assigns a curve to the group, whose representative curve is the closest, with respect to the chosen distance. Of course, the performance of this allocation rule deeply relies on the good performance of the pattern extraction.

One of the major drawbacks of this methodology are that first a high number of data are required in order to guarantee a good approximation of the geodesic distance at the core of this work. Actually, note that the number of observations, i.e the sampling rate of the manifold highly depends on the regularity of the manifold such that the assumption that the euclidean path between two observations follow approximatively the geodesic path. Hence, the data set should be carefully chosen for the manifold to be smooth enough. We point out that an enhancement could come from a prior registration procedure first applied to the curve and then the manifold warping procedure applied to the registered data.

The second drawback which may also be viewed as an advantage, is the following : the extracted pattern is a curve that belong to the observations. One the one hand, it may contains noise if the data are noisy observations, but on the other hand it thus guarantees that the pattern shares the mean properties and specifies of the observations. A solution when the noise must be removed is either to directly smooth the resulting pattern or to consider the neigbourhood of the extracted pattern with respect to the approximated geodesic distance and then use a kernel estimator with these observations to obtain a regularized mean curve.

Nevertheless, we promote this procedure when a large amount of data are available and when the sets of similar curves share a common behaviour which fully characterizes the observations, coming from an economic, physical or biological model for instance. This methods has been applied with success to a large amount of cases. Numerical packages for R or Matlab are available on request. where D is the following distance on ]b, c[ :

D (A i , α) = α A i
X ′ (a) da .

In the following, let A (i) i the ordered parameters. That is

A (1) < A (2) < • • • < A (n) .
Then, for a given α ∈]b, c[ such that A (j) < α < A (j+1) , we get that

C(α) = jD α, A (j) + j-1 i=1 iD A (i) , A (i+1)
+ (n -j)D α, A (j+1) +

n-1 i=j+1

(n -i)D A (i) , A (i+1) .

For the sake of simplicity, let n = 2q + 1. It follows that med(A) = A (q+1) . Moreover, let α = A (j) with j < q + 1. By symmetry, the case j > q + 1 will hold. Then, we rewrite C (α) as

C (α) = j-1 i=1 iD A (i) , A (i+1) + n-1 i=j (n -i)D A (i) , A (i+1)
and, by introducing A (q+1) , we get that

C(α) = j-1 i=1 iD A (i) , A (i+1) + q i=j iD A (i) , A (i+1) + q i=j (n -2i)D A (i) , A (i+1) + n-1 i=q+1 (n -i)D A (i) , A (i+1) .
Finally, we notice that

C(α) = C A (q+1) + q i=j (n -2i)D A (i) , A (

Fig 1 .

 1 Fig 1. Construction of a subgraph K ′ from Simulation (4) with n = 30 points. On the top left, a simulated data set. On the top right, the associated complete Euclidean graph K (Step 1). On the bottom left, the EMST associated with the complete graph K (Step 2). On the bottom right, the associated open balls and the corresponding subgraph K ′ (Step 3).

5.Fig 3 .

 3 Fig 3. On the left, a simulated data set of warped curves from Model (9) and an estimation of f with the mere mean (white line). On the right, the target function f (red dotted line), an estimation of the structural median by Manifold Warping (green solid line), an estimation obtained by Isomap (blue dot-dashed line), and an estimation obtained with the CAM procedure (black dashed line).

Fig 4 .

 4 Fig 4. On the left, the first landscape data. On the right, the CAM representative estimation (black dashed line) and the Manifod Warping estimation (green solid line).

Fig 5 .

 5 Fig 5. On the left, the second landscape data. On the right, the CAM representative estimation (black dashed line) and the Manifod Warping estimation (green solid line).

  7. Appendix. imsart-aoas ver. 2009/02/27 file: ManiWarp.hyper9470.tex date: May 27, 2011Proof of Theorem 1. Take µ = X(α) with α ∈]b, c[, we can write

  i+1) > C A (q+1) . imsart-aoas ver. 2009/02/27 file: ManiWarp.hyper9470.tex date: May 27, 2011And the result follows sinceµ 1 I = arg min X(α)∈C C(α) = X A (q+1) = X med(A) = f SM .

  Classical mean and k-nearest neighbors confusion matrices.

	Pixel	Mean classification	k-nn classification
	Reference Landscape1 Landscape2	Landscape1 Landscape2
	Landscape1	12	9	15	6
	Landscape2	0	20	2	18
			Table 2

imsart-aoas ver. 2009/02/27 file: ManiWarp.hyper9470.tex date: May 27, 2011


André Burnol

Jennifer Hellal

Claire Bény

Jacques-André Bourdoiseau

Catherine Crouzet

Catherine Joulian

Fabienne Battaglia-Brunet

Comparison of arsenic and mercury release linked to iron and sulfur biogeochemical processes

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Comparison of arsenic and mercury release linked to iron and sulfur biogeochemical processes André Burnol (1), Jennifer Harris-Hellal (1), Jacques-André Bourdoiseau (2), Catherine Crouzet (1), Claire Beny (1), and Fabienne Battaglia-Brunet (1)

(1) BRGM, EPI, Orléans, France (a.burnol@brgm.fr), (2) Université de La Rochelle, LEMMA, La Rochelle, France Understanding and predicting the mobility of redox-sensitive metals, such as arsenic (As) and mercury (Hg), in aquifers is important to reduce their adverse effects. Both these metals have a strong affinity for iron (hydr)oxides that control their sequestration and their bioavailability in soils and aquifers, in strong relation to microbial activity and physico-chemical conditions.

The aim of this work was to study the mineralogical transformations, induced by bacterial consortia collected in polluted environments, leading to arsenic or mercury release and speciation in controlled column systems under anoxic conditions.

The experimental setup consists of columns filled in the lower half with sand and in the upper half with a mixture of sand and iron oxides, initially doped with either arsenite As(III) or mercury Hg(II). The system is inoculated with the bacterial consortia of polluted aquifers and a nutritive medium is then circulated through the system in which biological and physico-chemical parameters were followed.

The effect of the activity of sulfate-reducing bacteria on arsenic mobilization was demonstrated at the interface between the Fe-rich zone and the zone without Fe: the biotic reduction of solid Fe(III) or the abiotic reduction by HS-are both possible. The increase of the energy gain for both reducers (iron and sulfate reducers) is also suggested and this phenomenon could explain the elevated As concentrations seen in certain aquifers (e.g. Bengal delta aquifers).

Other column experiments are presently being carried out to elucidate mercury redox transformations in similar conditions: the reduction of Hg(II) to volatile Hg0 and methylation to HgCH3-are both likely to take place in such anoxic conditions and could increase the mobility of mercury.


Ahmad Karfoul

Student member, IEEE Laurent Albera

Member, IEEE Gwénaël Birot

BLIND UNDERDETERMINED MIXTURE IDENTIFICATION BY JOINT CANONICAL DECOMPOSITION OF HO CUMULANTS

Keywords: Blind Underdetermined Mixture Identification (BUMI), underdetermined mixture, canonical decomposition, PARAFAC, INDSCAL, ICA, BSS

A new family of cumulant-based algorithms is proposed in order to blindly identify potentially underdetermined mixtures of statistically independent sources. These algorithms perform a joint CANonical Decomposition (CAND) of several higher order cumulants through a CAND of a 3-way array with special symmetries. These techniques are studied in terms of identifiability, performance and numerical complexity. From a signal processing viewpoint, the proposed methods are shown i) to have a better estimation resolution and ii) to be able to process more sources than the other classical cumulant-based techniques. Secondly, from a numerical analysis viewpoint, we deal with the identifiability and the convergence speed of several 3-way array decomposition procedures such as the symmetric ACDC scheme. We also show how to accelerate the iterative CAND algorithms by using differently the symmetries of the considered 3-way array. Next, from a multi-linear algebra viewpoint the paper aims at giving some insights on the uniqueness of a joint CAND of several hermitian multi-way arrays compared to the CAND of only one array. In addition, this allows us to extend the concept of Virtual Array (VA) to the case of a combination of several VA's.

I. INTRODUCTION

B LIND Underdetermined Mixture Identification (BUMI) is addressed in this paper. More particularly, we focus on noisy static mixtures of statistical independent sources. In such a context, the BUMI problem is defined by: Problem 1: Given a vector random process {x[k]} k∈N , find a (N × P ) mixing matrix A = [a 1 , • • • , a P ] potentially underdetermined such that:

∀ k ∈ N, x[k] = As[k] + ν[k] (1)
where a p is the p-th column vector of A, s

[k] = [s 1 [k], • • • , s P [k]
] T is a random vector whose P components are statistically independent and where ν[k] is a (N × 1) Gaussian noise vector, independent of s[k].

The resolution of problem 1 is of a great research interest given its importance in a great variety of signal processing Manuscript received October 10, 2008; accepted August 16, 2009. First published September 09, 2009; current version published January 13, 2010. The associate editor coordinating the review of this manuscript and approving it for publication was Dr. Alper Tunga Erdogan. This work is supported by the French Government, in one part by the ANR DECOTES Contract (06-BLAN-0074), in second part by the mv-EMD Contract(BLAN07-0314-02) and in third part by the Brittany region Contract ACOMB-SOLARE.

The authors are with the INSERM U642, Rennes F-35000, France, and the Université de Rennes 1, LTSI, Rennes F-35000, France. applications, such as blind source separation [START_REF] Ans | Adaptative neural architectures: Detection of primitives[END_REF], [START_REF] Kachenoura | ICA: a potential tool for BCI systems[END_REF] or estimation of direction of arrivals [START_REF] Ferreol | Fourth order blind identification of underdetermined mixtures of sources (FOBIUM)[END_REF]. Then, many solutions based on the concept of independent component analysis [START_REF] Comon | Independent Component Analysis, a new concept ?[END_REF], [START_REF] Jutten | Blind separation of sources, part I: An adaptive algorithm based on neuromimetic architecture[END_REF] were proposed in order to solve problem 1 when A is an overdetermined mixture, say when the number N of sensors is greater or equal than the number P of sources. For instance some methods iteratively maximize the opposite of the mutual information [START_REF] Bell | An information-maximization approach to blind separation and blind deconvolution[END_REF], [START_REF] Karvanen | Blind separation methods based on pearson system and its extensions[END_REF], [START_REF] Pham | Fast algorithm for mutual information based independent component analysis[END_REF] or the negentropy [20, chapter 8]. Other algorithms resort to cumulants, which allow for a measure of statistical independence [START_REF] Pesquet | Cumulant-based independence measures for linear mixtures[END_REF], [START_REF] Thirion | New criteria for blind signal separation[END_REF] less natural but easier to compute. In such a way, semi-algebraic methods were defined in order to benefit from the structure of the data cumulants as a function of the mixture. Particularly, some methods [START_REF] Belouchrani | A blind source separation technique using secondorder statistics[END_REF], [START_REF] Fety | Méthodes de traitement d'antenne adpatées aux radiocommunications[END_REF], [START_REF] Tong | AMUSE: a new blind identification algorithm[END_REF] exploit Second Order (SO) cumulants only, while others use either both SO and Fourth Order (FO) cumulants [START_REF] Cardoso | Blind beamforming for nongaussian signals[END_REF], [START_REF] Comon | Independent Component Analysis, a new concept ?[END_REF], [START_REF] Lathauwer | Independent component analysis and (simultaneous) third-order tensor diagonalisation[END_REF], [START_REF] Moreau | A generalization of joint-diagonalization criteria for source separation[END_REF], [START_REF] Pesquet | Cumulant-based independence measures for linear mixtures[END_REF], [START_REF] Thirion | New criteria for blind signal separation[END_REF] or FO cumulants only [START_REF] Albera | ICAR, a tool for blind source separation using fourth order statistics only[END_REF]. Contrary to the methods based on HO cumulants only [START_REF] Albera | ICAR, a tool for blind source separation using fourth order statistics only[END_REF], the methods using SO cumulants [START_REF] Belouchrani | A blind source separation technique using secondorder statistics[END_REF], [START_REF] Cardoso | Blind beamforming for nongaussian signals[END_REF], [START_REF] Comon | Independent Component Analysis, a new concept ?[END_REF], [START_REF] Fety | Méthodes de traitement d'antenne adpatées aux radiocommunications[END_REF], [START_REF] Lathauwer | Independent component analysis and (simultaneous) third-order tensor diagonalisation[END_REF], [START_REF] Moreau | A generalization of joint-diagonalization criteria for source separation[END_REF], [START_REF] Pesquet | Cumulant-based independence measures for linear mixtures[END_REF], [START_REF] Thirion | New criteria for blind signal separation[END_REF], [START_REF] Tong | AMUSE: a new blind identification algorithm[END_REF] attempt to identify an orthogonalized version of the mixture. In practice, the latter converge faster as a function of the number of samples provided that the noise covariance is known, while the formers are asymptotically insensitive to the presence of a Gaussian noise with unknown covariance.

Nevertheless, in practical fields, the overdetermined mixture assumption does not always hold. For example, in radiocommunications when the probability of receiving more sources than sensors increases with the reception bandwidth, then it can be necessary to solve problem 1 in the underdetermined case, for which P > N . Such a context may also occur in biomedical engineering, especially when the number of electrodes would be strongly reduced in order to propose an ambulatory recording system [START_REF] Poree | Blind source separation for ambulatory sleep recording[END_REF]. Identifiability results on underdetermined source separation were introduced a few years ago [START_REF] Eriksson | Complex random vectors and ica models: identifiability, uniqueness and separability[END_REF], [START_REF] Taleb | On underdetermined source separation[END_REF]. They showed that underdetermined mixtures could be identified up to a trivial matrix, that is a diagonal matrix and a permutation, provided that all sources are non-Gaussian. Many methods [START_REF] Albera | Blind Identification of Overcomplete Mixtures of sources (BIOME)[END_REF], [START_REF]Blind identification and source separation in 2x3 underdetermined mixtures[END_REF], [START_REF] Estevao | Blind channel identification algorithms based on the parafac decomposition of cumulant tensors: the single and multiuser cases[END_REF], [START_REF] Ferreol | Fourth order blind identification of underdetermined mixtures of sources (FOBIUM)[END_REF], [START_REF] Karfoul | Canonical decomposition of even higher order cumulant arrays for blind underdetermined mixture identification[END_REF], [START_REF] Lathauwer | Blind identification of underdetermined mixtures by simultaneous matrix diagonalization[END_REF], [START_REF] Lathauwer | Fourthorder cumulant-based blind identification of underdetermined mixtures[END_REF] have been developed since then. They use implicitly [START_REF]Blind identification and source separation in 2x3 underdetermined mixtures[END_REF], [START_REF] Ferreol | Fourth order blind identification of underdetermined mixtures of sources (FOBIUM)[END_REF], [START_REF] Lathauwer | Blind identification of underdetermined mixtures by simultaneous matrix diagonalization[END_REF], [START_REF] Lathauwer | Fourthorder cumulant-based blind identification of underdetermined mixtures[END_REF] or explicitly [START_REF] Estevao | Blind channel identification algorithms based on the parafac decomposition of cumulant tensors: the single and multiuser cases[END_REF], [START_REF] Karfoul | Canonical decomposition of even higher order cumulant arrays for blind underdetermined mixture identification[END_REF] a CANonical Decomposition (CAND) [START_REF] Lathauwer | A link between the canonical decomposition in multilinear algebra and simultaneous matrix diagonalization[END_REF] of either SO [START_REF] Lathauwer | Blind identification of underdetermined mixtures by simultaneous matrix diagonalization[END_REF], FO [START_REF]Blind identification and source separation in 2x3 underdetermined mixtures[END_REF], [START_REF] Estevao | Blind channel identification algorithms based on the parafac decomposition of cumulant tensors: the single and multiuser cases[END_REF], [START_REF] Ferreol | Fourth order blind identification of underdetermined mixtures of sources (FOBIUM)[END_REF], [START_REF] Lathauwer | Fourthorder cumulant-based blind identification of underdetermined mixtures[END_REF] or Higher Order (HO) [START_REF] Albera | Blind Identification of Overcomplete Mixtures of sources (BIOME)[END_REF], [START_REF] Karfoul | Canonical decomposition of even higher order cumulant arrays for blind underdetermined mixture identification[END_REF] cumulants of the data. They decompose either one cumulant array [START_REF] Albera | Blind Identification of Overcomplete Mixtures of sources (BIOME)[END_REF], [START_REF] Estevao | Blind channel identification algorithms based on the parafac decomposition of cumulant tensors: the single and multiuser cases[END_REF], [START_REF] Karfoul | Canonical decomposition of even higher order cumulant arrays for blind underdetermined mixture identification[END_REF], [START_REF] Lathauwer | Fourthorder cumulant-based blind identification of underdetermined mixtures[END_REF] or several ones at a given statistical order [START_REF]Blind identification and source separation in 2x3 underdetermined mixtures[END_REF], [START_REF] Lathauwer | Blind identification of underdetermined mixtures by simultaneous matrix diagonalization[END_REF]. In fact, CAND is an extension of the Singular Value Decomposition (SVD) to the case of multi-way arrays without requiring any orthogonality constraint. Other algorithms [START_REF] Comon | Blind identification of underdetermined mixtures based on the characteristic function[END_REF], [START_REF] Taleb | An algorithm for the blind identification of N independent signal with 2 sensors[END_REF] take advantage of pseudo-cumulants, which are given by the (weighted) derivatives of the second characteristic function at some points different from the origin. However, all previous methods have some limits. Indeed, under some computational cost constraints and for N ≥ 3 they do not succeed in exploiting all the information contained in the statistical tools they used, except the FOOBI2 method [START_REF] Lathauwer | Fourthorder cumulant-based blind identification of underdetermined mixtures[END_REF]. But the latter cannot process more than N 4 (N 4 ≤ N 2 ) sources, where N 4 is the number of sensors, which virtually appear when a Fourth Order (FO) cumulant is used in practical contexts such as radiocommunications [START_REF] Chevalier | On the virtual array concept for the fourth-order direction finding problem[END_REF].

In order to overcome these limitations and due to the attractive properties of cumulants which make them easy to handle, we propose in this paper a new family of BUMI cumulant-based algorithms. These algorithms perform a joint CAND of several HO cumulants through a CAND of a special 3-way array with special symmetries. These techniques are studied in terms of identifiability, performance and numerical complexity. From a signal processing viewpoint, the proposed methods are shown i) to have a better estimation resolution and ii) to be able to process more sources than the other classical cumulant-based techniques. Secondly, from a numerical analysis viewpoint, we deal with the identifiability and the convergence speed of several 3-way array decomposition procedures such as the symmetric ACDC scheme [START_REF] Yeredor | Non-orthogonal joint diagonalization in the leastsquares sense with application in blind source separation[END_REF]. We also show how to accelerate the iterative CAND algorithms by using differently the symmetries of the considered 3-way array. Next, from a multi-linear algebra viewpoint the paper aims at giving some insights on the uniqueness of a joint CAND of several hermitian multi-way arrays compared to the CAND of only one array. In addition, this allows us to extend the concept of Virtual Array (VA) to the case of a combination of several VA's.

Thus, cumulants are presented in section II. In addition, some basic definitions in multi-linear algebra and a tensor reformulation of problem 1 are given in section III. In Section IV, the new family of methods is presented. Next, detailed identifiability and numerical complexity studies are proposed in sections V and VI, respectively. Finally, simulation results and conclusion are pointed out in sections VII and VIII, respectively.

II. STATISTICAL TOOLS

Cumulants have many attractive properties which make them easy to handle [START_REF] Kachenoura | ICA: a potential tool for BCI systems[END_REF], [START_REF] Mccullagh | Tensor Methods in Statistics[END_REF]. In particular, they can be expressed as a function of moments of smaller order, which is shown by the Leonov and Shiryaev formula [START_REF] Mccullagh | Tensor Methods in Statistics[END_REF], [START_REF] Spooner | The cumulant theory of cyclostationarity time-series -Part. II: Development and applications[END_REF]. The latter formula is well-known in the case of non-delayed cumulants [START_REF] Mccullagh | Tensor Methods in Statistics[END_REF], [START_REF] Spooner | The cumulant theory of cyclostationarity time-series -Part. II: Development and applications[END_REF]. But since in the following we will use the time dependence of each source through cumulants, we propose to rewrite this formula hereafter in the case of delayed cumulants. To this end, let's give some notations and definitions.

Definition 1: A partition of a set Ω(q) = {1, 2, • • • , q} of integers is a set of M (1 ≤ M ≤ q) nonempty
disjunctive subsets whose union is equal to Ω(q). Now, let Ω (j) m (M, q) be the m-th (1 ≤ m ≤ M ) subset of the j-th (1 ≤ j ≤ J(M, q)) partition of size M of Ω(q). This notation is illustrated in the following example:

Example 1: The different partitions of set {1, 2, 3} are:

Ω (1) 1 (1, 3) = {{1, 2, 3}} {Ω (1) 1 (2, 3), Ω (1) 2 (2, 3)} = {{1, 2}, {3}} {Ω (2) 1 (2, 3), Ω (2) 2 (2, 3)} = {{1, 3}, {2}} {Ω (3) 1 (2, 3), Ω (3) 2 (2, 3)} = {{2, 3}, {1}} {Ω (1) 1 (3, 3), Ω (1)
2 [START_REF] Ans | Adaptative neural architectures: Detection of primitives[END_REF][START_REF] Ans | Adaptative neural architectures: Detection of primitives[END_REF], Ω Thus, we have the following proposition: Proposition 1: Let {x[k]} k∈N be an N -dimensional random process. Its q-th order delayed cumulants are given by:

Cum{x n1 [k], x n2 [k-τ 1 ] ξ(1) , • • •, x nq [k-τ q-1 ] ξ(q-1) } = q M =1 (-1) (M -1) (M -1)! J(M,q) j=1 M m=1 E[ ω∈Ω j m (M,q) x nω [k-τ ω-1 ] ξ(ω-1) ] (2)
where E[.] is the mathematical expectation,

x n [k -τ ] is the n-th component of vector x[k -τ ], τ is a set {τ 1 , • • • , τ q-1 }
of q-1 delays, ξ(ω) equals ±1 with the conventions, ξ(0) = 1, τ 0 = 0, x 1 = x and x -1 = x * denoting by * the complex conjugate operator.

It is noteworthy that the estimation of cumulants using equation ( 2) is appropriate when the processed sources are stationary. But, in practice, sources are rarely so. In such a case, time averaged cumulants have to be used instead of (2), such as

C nq+1,••• ,n2q n1,••• ,nq,x [τ ] =< Cum{x n1 [k], x n2 [k-τ 1 ] ξ(1) , • • •, x nq [k-τ q-1 ] ξ(q-1) } > (3)
where < . > is the discrete-time temporal mean operation [START_REF] Albera | Blind Identification of Overcomplete Mixtures of sources (BIOME)[END_REF]. Such a statistics enjoy the multilinearity property and they can be used as a measure of independence. Thus for sake of convenience, time averaged cumulant are considered throughout this paper.

III. TOWARD A MULTI-WAY ARRAY FORMULATION

This section aims at showing how problem 1 can be solved using a joint CAND of several HO cumulant arrays through a multilinear algebra framework. Recall that multilinear algebra is defined as the algebra of q-way arrays (q > 2), that is, the arrays whose elements are addressed by more than two indices. It is noteworthy that multilinear algebra can be related to tensor algebra through a multilinear map, called Segre map [START_REF] Comon | Symmetric tensors and symmetric tensor rank[END_REF]. Now some definitions are necessary in order to understand how the BUMI problem can be reformulated in a multilinear algebra framework. First let's begin by introducing a very useful product operator:

Definition 2: The outer product T = u (1) • • • • • u (q) of q (q ≥ 2) vectors u (i) ∈ C Ni (1 ≤ i ≤ q) is a rank-1 array
n1 • • • u (q)
nq .

Then the matrix rank notion can be extended to the case of q-way (q > 2) arrays [START_REF] Kruskal | Three-way arrays : Rank and uniqueness of trilinear decompositions, with application to arithmatic complexity and statistics[END_REF], [START_REF]Rank, decomposition, and uniqueness for 3-way and N -way arrays[END_REF], as following: Definition 3: The rank of an array T ∈ C N1ו••×Nq , denoted by rk(T ), is the minimal number of rank-1 arrays belonging to C N1ו••×Nq that yield T in a linear combination.

Despite the similarities between matrices and q-way (q > 2) arrays, the rank of the latter ones may exceed their dimensions. Now, let's see how the matrix SVD can be extended to the case of multi-way arrays.

Definition 4: The CAND of a q-way (q > 2) array T ∈ C N1ו••×Nq is the linear combination of P = rk(T ) rank-1 q-way arrays belonging to C N1ו••×Nq , that is:

T = P p=1 λ p u (1) p • u (2) p • • • • • u (q) p (4)
where for each integer p, λ p is a complex number and vectors u (i)

p are complex N i -dimensional vectors.
Nevertheless, the q matrices

U (i) = [u (i) 1 , • • • , u (i)
P ], given by (4) and well-called loading matrices of T , are not necessarily unitary contrary to the matrix SVD. Figure 1 shows an example of CAND in the case of a 3-way array. In addition, equation [START_REF] Bell | An information-maximization approach to blind separation and blind deconvolution[END_REF] shows that the q vectors u (i) p can be arbitrarily scaled as long as their product remains the same. Also the rank-1 terms can be permuted without modifying the sum. Then we call the CAND of a q-way array (q > 2) unique when it is only subject to these trivial indeterminacies. Important results on CAND uniqueness were obtained by Kruskal [START_REF] Kruskal | Three-way arrays : Rank and uniqueness of trilinear decompositions, with application to arithmatic complexity and statistics[END_REF]. He showed that the CAND of a q-way array T is essentially unique when the following sufficient condition is satisfied:

2P + q -1 ≤ q i=1 k(U (i) ) (5)
where k(T ) denotes the kruskal rank of T defined by the maximal number r such that any set of r columns of T is linearly independent [START_REF] Lathauwer | Blind identification of underdetermined mixtures by simultaneous matrix diagonalization[END_REF]. Kruskal's condition was firstly proposed for real 3-way arrays. Later it was extended to arrays with order higher than three [START_REF] Sidiropoulos | On the uniqueness of the multilinear decomposition of n-way arrays[END_REF] and for complex HO arrays [START_REF] Sidiropoulos | Parallel factor analysis in sensor array processing[END_REF]. The CAND uniqueness of specific 3-way arrays was also studied. For instance the case when one loading matrix has a full column rank [START_REF] Jiang | Kruskal's permutation lemma and the identification of CANDECOMP/PARAFAC and bilinear models with constant modulus constraints[END_REF] and the case where the HO array is tall in one direction [START_REF] Lathauwer | A link between the canonical decomposition in multilinear algebra and simultaneous matrix diagonalization[END_REF]. In practice, the computation of CAND may require to handle unfolding matrices or vectors of multi-way arrays. A useful way to unfold square even order multi-way arrays is described below.

Definition 5: Let T be a square 2q-way (q > 1) array of dimension N , and let be the integer part of q/2. Then the (i, j)-th component of the unfolding matrix mat 1 (T ) of size (N q × N q ) is given by:

(mat 1 (T )) i,j = T n1,••• ,n q-,n q-+1 ,••• ,nq,nq+1,••• ,n 2q-,n 2q-+1 ,••• ,n2q (6)
where i=(n

1 -1)N q-1 + • • • + (n q--1)N + (n 2q-+1 - 1)N -1 +• • •+(n 2q -1)N +n 2q and j=(n q+1 -1)N q-1 +• • •+ (n 2q--1)N + (n q-+1 -1)N -1 + • • • + (n q -1)N + n q .
For instance, let T be a square 6-way array of dimension N . Then = 1 and the (i, j)-th component of the unfolding matrix mat 1 (T ) is given by:

(mat 1 (T )) i,j = T n1,n2,n3,n4,n5,n6 (7)
where i = (n 1 -1)N 2 + (n 2 -1)N + n 6 and j = (n 4 -1)N 2 + (n 5 -1)N + n 3 .

Another way to transform higher order arrays into matrices is presented hereafter and is illustrated by figure 2.

Definition 6: Let T ∈ C N1ו••×Nq be a q-way (q > 2)
array. Then the (n i , m)-th component of the unfolding matrix mat

(i) 2 (T ) ∈ C Ni×Ni+1•••NqN1•••Ni-1 associated to the i-th mode (1 ≤ i ≤ q)
of T is given by:

(mat

(i) 2 (T )) ni,m = T n1,••• ,ni-1,ni,ni+1,••• ,nq (8)
where

m = (n i+1 -1)N i+2 • • • N q N 1 • • • N i-1 + (n i+2 - 1)N i+3 • • • N q N 1 • • • N i-1 +• • • +(n q -1)N 1 N 2 • • • N i-1 +(n 1 - 1)N 2 N 3 • • • N i-1 + (n 2 -1)N 3 N 4 • • • N i-1 +• • •+ n i-1 .
In addition, let's introduce the following q-way array-to-vector transformation:

Definition 7: Let T ∈ C N1ו•
•×Nq be a q-way (q > 1) array.

Then the components of the N 1 • • • N q -dimensional vector vec(T ) are defined by:

vec(T ) (n1-1)N2•••Nq+(n2-1)N3•••Nq+•••+nq = T n1,••• ,nq (9)
Conversely, let unvec be the inverse operator such as unvec (vec(T )) = T .

Hermitian q-way arrays are considered in this paper, hence the following definition: Definition 8: A complex q-way (q ≥ 2) array T is called hermitian when its CAND takes the following form:

T = P p=1 λ p u ξ(1) p • u ξ(2) p • • • • • u ξ(q) p (10)
where the P values λ p are real and where application ξ(i) is defined in proposition 1, that is, if the (P × P ) diagonal matrix Λ = diag{[λ 1 , • • • , λ P ]} is real and if all the loading matrices of T are equal within a conjugate transformation.

Inserting the observation model ( 1) into (2) and using i) the statistical independence of the P source components, ii) the statistical independence between the sources and the noise and iii) the noise Gaussian distribution, we get for any integer q ≥ 2:

C 2q, x [τ ] = P p=1 C 2q,sp [τ ] a •q p • a * p •q (11)
where

C 2q,sp [τ ] =< Cum{s p [k], s p [k -τ 1 ], . . . , s p [k - τ q-1 ], s p [k -τ q ] * , . . . , s p [k -τ 2q-1 ] * } > and where a •q
is the outer product of q vectors a. Formula (11) conveys the multilinearity property enjoyed by cumulant and moment arrays [START_REF] Kachenoura | ICA: a potential tool for BCI systems[END_REF]. Then, provided that the P sources have at least M non-zero 2q-th order delayed marginal cumulants, problem 1 can be reformulated as following:

Problem 2: Given M 2q-th (q ≥ 2) order cumulant arrays C 2q, x [τ (m) ] of {x[k]} k∈Z associated to a set {τ (m) } of (2q - 1)-dimensional time lags, find its joint CAND.
But, is it really possible to compute such a joint CAND, especially when P is greater than the maximum number of sources which can be processed by the classical cumulantbased techniques [START_REF] Albera | ICAR, a tool for blind source separation using fourth order statistics only[END_REF], [START_REF] Albera | Blind Identification of Overcomplete Mixtures of sources (BIOME)[END_REF], [START_REF] Belouchrani | A blind source separation technique using secondorder statistics[END_REF], [START_REF] Cardoso | Blind beamforming for nongaussian signals[END_REF], [START_REF] Comon | Independent Component Analysis, a new concept ?[END_REF], [START_REF]Blind identification and source separation in 2x3 underdetermined mixtures[END_REF], [START_REF] Comon | Blind identification of underdetermined mixtures based on the characteristic function[END_REF], [START_REF] Estevao | Blind channel identification algorithms based on the parafac decomposition of cumulant tensors: the single and multiuser cases[END_REF]- [START_REF] Fety | Méthodes de traitement d'antenne adpatées aux radiocommunications[END_REF], [START_REF] Karfoul | Canonical decomposition of even higher order cumulant arrays for blind underdetermined mixture identification[END_REF], [START_REF] Lathauwer | Blind identification of underdetermined mixtures by simultaneous matrix diagonalization[END_REF], [START_REF] Lathauwer | Fourthorder cumulant-based blind identification of underdetermined mixtures[END_REF], [START_REF] Lathauwer | Independent component analysis and (simultaneous) third-order tensor diagonalisation[END_REF], [START_REF] Moreau | A generalization of joint-diagonalization criteria for source separation[END_REF], [START_REF] Pesquet | Cumulant-based independence measures for linear mixtures[END_REF], [START_REF] Thirion | New criteria for blind signal separation[END_REF], [START_REF] Tong | AMUSE: a new blind identification algorithm[END_REF]? Yes, it suffices to choose the cumulant order (2q) greater or equal than four and six using the approaches proposed in sections IV-B and IV-A, respectively. Nevertheless, for the sake of simplicity and readability, we limited the presentation of our methods and their identifiability study to a fixed statistical order 2q equal to six and identical for all the proposed algorithms. An extension to all other order would be then straightforward from sections IV and V.

IV. TOWARD A 3-WAY ARRAY DECOMPOSITION

We present in this section a new family of methods, named SIBI (SIxth order Blind mixture Identification), in order to solve problem 2 when q = 3. This family is essentially composed of two classes of techniques, namely the SAD (Semi-Algebraic Decomposition) and ALS (Alternate Least Square) SIBI approaches. Both classes can be described through two steps.

The first step consists in identifying the P rank-1 arrays A (3) p given by:

∀1 ≤ p ≤ P, A (3) p = a p • a p • a * p ( 12
)
from a CAND of a special 3-way array H x of size (N 3 ×N 3 × M ) whose ( 1 , 2 , m)-th entry is given by:

H 1, 2,m,x = C n4,n5,n6 n1,n2,n3,x [τ (m) ] (13)
where

1 = n 6 + (n 2 -1)N + (n 1 -1)N 2 and 2 = n 3 + (n 5 - 1)N + (n 4 -1)N 2 .
On the other hand, H x can also be built by stacking up the M hexacovariance matrices

C 6,x [τ (m) ] = mat 1 (C 6, x [τ (m) ]) of size (N 3 × N 3
) in the third direction. From ( 11) and ( 13), we get:

Proposition 2: Given the P M -dimensional vectors b p = [C 6,sp [τ (1) ], • • •, C 6,sp [τ (M ) ]] T , H x has the following CAND:

H x = P p=1 a (3) p • a (3) p * • b p ( 14
)
where a

(3) p = a p ⊗ a p ⊗ a * p with ⊗ the Kronecker product operator and where

A (3) = [a (3) 1 , • • • , a (3)
P ], A (3) * and B = [b 1 , • • • , b P ] are the three loading matrices of H x .
Once the P N 3 -dimensional vectors a (3) p have been identified from the CAND of H x , each array A (3) p is found by computing A (3) p = unvec(a

(3) p ).
The second step consists in canonically decomposing each rank-1 array A (3) p in order to identify each column vectors a p of A. It is noteworthy that the way to implement this second step is the same for all SIBI methods and is based on the following proposition:

Proposition 3: Each (N × N 2 ) unfolding matrix A (3) p = mat (1) 2 (A (3)
p ) is the horizontal concatenation of N (N ×N ) matrices, which are of the form γ a p a p H with γ ∈ C.

In fact, γ is equal to the (n, p)-th component, A n,p , of A when the n-th horizontal matrix block of A (3) p is concerned. Consequently, the eigenvector jointly associated with the largest eigenvalue of each of the N matrix blocks of A (3) p allows to estimate a p . Such an eigenvector can be computed using a joint diagonalization procedure [START_REF]Jacobi angles for simultaneous diagonalization[END_REF], [START_REF] Chabriel | Algebraic joint zero-diagonalization and blind sources separation[END_REF], [START_REF] Yeredor | Non-orthogonal joint diagonalization in the leastsquares sense with application in blind source separation[END_REF], [START_REF]On using exact joint diagonalization for noniterative approximate joint diagonalization[END_REF].

A. The SAD-SIBI approach

A semi-algebraic solution is proposed hereafter in order to find decomposition [START_REF] Comon | Blind identification of underdetermined mixtures based on the characteristic function[END_REF]. It consists in handling the M frontal slices H x (:, :, m) = C 6,x [τ (m) ] instead of H x itself. Then according to [START_REF] Cardoso | Blind beamforming for nongaussian signals[END_REF] and [START_REF] Comon | Independent Component Analysis, a new concept ?[END_REF] we can easily show that: Proposition 4: Each hexacovariance matrix C 6,x [τ (m) ] has the following algebraic structure:

C 6,x [τ (m) ] = A (3) ζ 6,s [τ (m) ] A (3)H ( 15
)
where

ζ 6,s τ (m) =diag{[C 6,s1 [τ (m) ], • • • , C 6,s P [τ (m) ]]}.
Since A (3) is not generally orthogonal, equation [START_REF] Eriksson | Complex random vectors and ica models: identifiability, uniqueness and separability[END_REF] does not really correspond to the eigenvalue decomposition of C 6,x τ (m) , except when this latter has been prewhitened. Let's see how to perform such a prewhitening. Under assumptions 1 and 2 given by: Assumption 1: The (N 3 × P ) matrix A (3) is full column rank.

and:

Assumption 2: The (M ×P ) matrix B is full column rank.

respectively, the following proposition holds:

Proposition 5: There is at least one rank-P positive semidefinite matrix, C 6,x , which is a linear combination of the M matrices C 6,x [τ (m) ].

Its proof is given in appendix A. In some practical contexts, for example when the P non-delayed source marginal Sixth Order (SixO) cumulants have the same sign , such as in radiocommunications, a simple linear combination is given by C 6,x = C 6,x [0] where 0 = {0, 0, 0, 0, 0}. Otherwise, in the general case, we have to resort to numerical methods such that those presented in [START_REF] Tong | A finite-step global convergence algorithm for the parameter estimation of multichannel MA processes[END_REF], [START_REF] Zaidi | Positive definite combination of symmetric matrices[END_REF] to compute C 6,x [τ ]. In addition, it is noteworthy that, matrix C 6,x may be non-hermitian and then should be replaced by (C 6,x +C 6,x H )/2 in order to ensure the existence of its eigen-value decomposition. For the sake of simplicity, we will assume in the sequel that matrix C 6,x is hermitian. Then proposition 5 ensures that C 6,x has a square root, which can be computed as following:

C 1/2 6,x = E s Λ 1/2 s ( 16
)
where E s is the (N 3 × P ) unitary matrix whose columns are the eigenvectors associated with the P largest eigenvalues of C 6,x sorted in the diagonal matrix Λ s . Since two square roots of a positive semidefinite matrix are equal up to a unitary matrix, we have:

C 1/2 6,x = A (3) ζ6,s 1/2 U • H ( 17
)
where ζ 6,s , ζ 1/2 6,s , U • are the linear combination described in proposition 5 of the M matrices (ζ 6,s [τ (m) ]+ζ 6,s [τ (m) ] H )/2, a square root of the positive definite diagonal matrix ζ 6,s and a unitary matrix. Consequently, the inverse matrix, C -1/2 6,x , of C 1/2 6,x allows for a prewhitening of each matrix C 6,x [τ (m) ],

say each matrix C 6,x [τ (m) ] = C -1/2 6,x C 6,x [τ (m) ](C -1/2 6,x ) H
is unitarily similar to a diagonal matrix indexed by m as described by the following equation:

C 6,x [τ (m) ] = U • D[τ (m) ]U • H (18)
where D τ (m) = ζ 6,s τ (m) ζ-1 6,s is diagonal. Thus, a joint orthogonal diagonalization with the well-known semialgebraic JAD algorithm [START_REF]Jacobi angles for simultaneous diagonalization[END_REF] allows for an identification of U • up to a unitary trivial matrix. Nevertheless, U • is identifiable provided that, according to the results of [START_REF] Belouchrani | A blind source separation technique using secondorder statistics[END_REF], the M (2q -1)dimensional delays τ (m) are chosen such that: Assumption 3: For each couple (p, q) of sources, there exists at least one set τ (m) of delays for which:

C 6,sp [τ (m) ]/ ζ6,sp,p = C 6,sq [τ (m) ]/ ζ6,sq,q (19)
It is noteworthy that assumption 3 is true under assumption 2. Indeed, assumption 2 implies that the (M ×P ) matrix B ζ6,s is full column rank and so that all columns of B ζ6,s are different. Hence equation [START_REF] Givens | Computation of plain unitary rotations transforming a general matrix to triangular form[END_REF]. Thus, under assumption 2, the JAD algorithm diagonalizes the M matrices C 6,x [τ (m) ] by minimizing the following criterion:

Ψ 1 (U ) = M m=1 P p =q=1 |(U H C 6,x [τ (m) ]U ) p,q | 2 (20)
with respect to (w.r.t.) to a unitary matrix U expressed as the product of Givens rotations [START_REF] Givens | Computation of plain unitary rotations transforming a general matrix to triangular form[END_REF]. Once matrix U • has been identified, matrix A (3) can be computed up to a trivial matrix multiplying matrix C

1/2 6,x (17)
by the estimate of U • . This approach will be referred to as SAD-SIBI JAD in the sequel.

B. The ALS-SIBI approaches

Another way to find decomposition [START_REF] Comon | Blind identification of underdetermined mixtures based on the characteristic function[END_REF] consists in iteratively and alternatively looking for the three loading matrices A (3) , A (3) * and B of H x in the least square sense. More precisely, from proposition 2 and definition 6, we have: Proposition 6: Given the three unfolding matrices H (1) x = mat [START_REF] Albera | ICAR, a tool for blind source separation using fourth order statistics only[END_REF] 2 (H x ), H (2) x = mat

2 (H x ) and H (3) x = mat

2 (H x ) of H x , we get H (1) x =A (3) (A (3) * B) T , H (2) x =A (3) * (B A (3) ) T and H (3) x = B(A (3) A (3) * ) T where is the Khatri-Rao product operator (columns-wise kronecker product) [START_REF] Smilde | Multi-way Analysis[END_REF].

So, the following criterion can be minimized w.r.t. each of the three matrices T (1) , T (2) and T (3) in order to compute the loading matrices A (3) , A (3) * and B of H x up to a trivial matrix:

Ψ (1)
3 (T (1) ,T (2) ,T (3) ) = H (1) x -T (1) (T (2) T (3) ) T 2 F (21) using the Alternate Least Square (ALS) procedure [START_REF] Smilde | Multi-way Analysis[END_REF] where . F denotes the Frobenius norm, giving rise to the ALS-SIBI UNS method. In fact, at each iteration it, we compute the three estimates T (i) it (1 ≤ i ≤ 3) one at a time by fixing both the others as described below:

T (1)
it = H (1) x ((T

(2) it-1 T (3) it-1 ) ) T T (2)
it = H (2) x ((T

(3) it-1 T (1) it ) ) T T (3)
it = H (3) x ((T

(1) it T (2) it ) ) T (22)
where T denotes the pseudo-inverse of T . Consequently, matrices T

(1)

it , T (2)
it and T

(3)

it should converge to matrices

A (3) , A (3)
* and B, respectively, provided that the CAND ( 14) is unique, which is ensured under assumptions 2 and the following assumption:

Assumption 4: The (N 6 × P ) matrix A (3) A (3) * is full column rank.

It is noteworthy that one is entitled to expect some convergence cycles in which all factors will evolve in the same direction [START_REF] Rajih | Enhanced line search: A novel method to accelerate PARAFAC[END_REF] and consequently a slow convergence of the ALS-SIBI UNS method. An improved version of the ALS scheme, named Line Search (LS) ALS, was proposed (see [START_REF] Rajih | Enhanced line search: A novel method to accelerate PARAFAC[END_REF] and the references therein) in order to exit from these convergence cycles. This procedure is based on the relaxation factor concept, which consists in extrapolating the estimated loading matrices many iterations ahead. Then, we decided to apply this LS procedure to the ALS-SIBI UNS method, giving rise to the LSALS-SIBI UNS algorithm. The extrapolation of the i-th estimated loading matrix of H x is given by:

T (i) new = T (i) it-2 + β (i) (T (i) it-1 -T (i) it-2 ) (23)
where T (i) new is the matrix that will be used in the it-th iteration instead of

T (i) it-1 , the term (T (i) it-1 -T (i)
it-2 ) defines the direction of the cycle and β (i) is the relaxation factor associated to the i-th loading matrix. In practice, the triplet (β (1) , β (2) , β (3) ) of relaxation factors is chosen such that β (1) = beta (2) = β (3) = it 1/n where n is a fixed integer [START_REF] Rajih | Enhanced line search: A novel method to accelerate PARAFAC[END_REF]. Recent works [START_REF] Rajih | Enhanced line search: A novel method to accelerate PARAFAC[END_REF] showed that the convergence speed of the LSALS procedure could be also increased, especially when loading matrices are strongly correlated. This induced an enhanced version of the LSALS algorithm, named Enhanced Line Search (ELS)ALS, in the real [START_REF] Rajih | Enhanced line search: A novel method to accelerate PARAFAC[END_REF] and the complex [START_REF] Lathauwer | An enhanced line search scheme for complex-valued tensor decompositions. Application in DS-CDMA[END_REF] fields. More precisely, the ELS procedure consists, at each iteration it of the LSALS method, in finding the optimal triplet (β (1) , β (2) , β (3) ). So, when the ELS scheme is applied to LSALS-SIBI UNS giving rise to the ELSALS-SIBI UNS method, the following criterion:

Ψ (2)
3 (β) = H (1) x -G

(it) (G (1)
(it) G

(it) ) T 2 F (3)
is minimized at each iteration it w.r.t. β (1) , β (2) and β (3) where

G (i) (it) = T (i) it-2 +β (i) (T (i) it-1 -T (i)
it-2 ). The ELS scheme allows to attain the final solution of a given cycle in only one step. Note that the ELS procedure was implemented in the case of three equal relaxation factors only, which is suboptimal but less time-consuming.

However, by applying the previous ALS-like procedures to H x , it amounts to assume a total independence between the loading matrices, which is clearly wrong in our case. Thus an alternative way of the ALS minimization would consists in using the symmetries that exist between some loading matrices of the considered multi-way array. In our case, the first and the second loading matrices of H x are equal up to a complex conjugate as shown in [START_REF] Comon | Blind identification of underdetermined mixtures based on the characteristic function[END_REF]. The ACDC-like procedures [START_REF] Yeredor | Non-orthogonal joint diagonalization in the leastsquares sense with application in blind source separation[END_REF], [39, Appendix 7.2] exploit such symmetries. In fact the ACDC algorithm [START_REF] Yeredor | Non-orthogonal joint diagonalization in the leastsquares sense with application in blind source separation[END_REF] minimizes the following weighted least square criterion w.r.t. V and the M diagonal

(P ×P ) matrices Λ[τ (m) ]: Ψ 2 (V,Λ) = M m=1 w m C 6,x [τ (m) ]-VΛ[τ (m) ]V H 2 F (25)
where w m ∈ R + * . More particularly, the last minimization problem is devised into two steps. The first one, named Alternating Columns (AC), looks for the minimizer of ( 25) w.r.t. a selected column vector v p of V keeping the p-1 other column vectors and the M matrices Λ[τ (m) ] fixed. The second step, called Diagonal Centers (DC), minimizes equation ( 25) w.r.t. the diagonal matrices Λ[τ (m) ] keeping V fixed. So, the ACDC algorithm alternates, until convergence, between the P consecutive AC sweeps, one by column vector of V , and the DC step. Note that, an efficient way to initialize the ACDC algorithm was recently proposed [START_REF]On using exact joint diagonalization for noniterative approximate joint diagonalization[END_REF]. In the following, we will refer to as the ALS-SIBI ACDC method when the ACDC procedure is used. But, are the ACDC-like procedures [START_REF] Yeredor | Non-orthogonal joint diagonalization in the leastsquares sense with application in blind source separation[END_REF], [START_REF] Rajih | Enhanced line search: A novel method to accelerate PARAFAC[END_REF]Appendix 7.2] the only way to decompose H x in a least square sense using its symmetries? No, there is another way such as the solution that we propose hereafter. Furthermore, we combined this new 3-way array CAND method with the ELS scheme in order to accelerate its convergence, giving rise to the ELSALS SYM procedure. In fact, the latter is mainly based on the following proposition: Proposition 7: Let H x be the 3-way array of size (N 3 × N 3 × P ) whose third (P × N 6 ) unfolding matrix is given by H

x = B H (3) x . Then the p-th (N 3 ×N 3 ) frontal slice of H x , H x (:, :, p), is a rank-1 matrix equal to a

(3) p a (3) p H .
The proof is straightforward. Indeed, from proposition 2 we get H

(3) x = I P (A (3) A (3) * ) T , which implies that the pth frontal slice of H x is equal to A (3) diag{I P (p, :)}A (3)H where I P (p, :) is the p-th row of the (P ×P ) identity matrix I P . Consequently, the eigenvector associated with the largest eigenvalue of H x (:, :, p) allows for an identification of the p-th column vector of A (3) up to a scalar factor. Now, let's see how the ELSALS SYM scheme uses this result. First the estimate, T

it-1 , of the third loading matrix, B, of H x is computed at iteration it -1 from equation [START_REF] Jutten | Blind separation of sources, part I: An adaptive algorithm based on neuromimetic architecture[END_REF]. Next, using proposition 7, we compute the (P × N 6 ) matrix

T (3) it-1 = T (3) it-1 H (3)
x and build the 3-way array T it-1 of size (N 3 × N 3 × P ) whose third (P × N 6 ) unfolding matrix is equal to

T (3) it-1 .
So, the eigenvector associated with the largest eigenvalue of the p-th frontal slice of T it-1 gives the p-th column vector of

T (1) it-1 . Thus, T (2) it-1 is computed as the conjugate of T (1) it-1 .
Eventually, the ELS procedure is used in order to accelerate the convergence of the proposed algorithm. Note that all these steps are repeated until convergence.

The convergence speed of the five proposed ALS-SIBI approaches, namely ALS-SIBI UNS , LSALS-SIBI UNS , ELSALS-SIBI UNS , ALS-SIBI ACDC and ELSALS-SIBI SYM , is studied now. More precisely, we display the median value of Ψ (1) 3 [START_REF] Jiang | Kruskal's permutation lemma and the identification of CANDECOMP/PARAFAC and bilinear models with constant modulus constraints[END_REF] over the independent trials as a function of the number of iterations at the output of the five ALS-SIBI methods. In fact, the latter methods are used to canonically decompose the 3-way array H x built from equation ( 14) with M = 16

5-dimensional delays τ

(m) = [τ (m) 1 , τ (m) 2 , τ (m) 3 , τ (m) 3
, 0] and the estimated SixO delayed marginal cumulants of P = 3 angularly close QPSK sources. The latter sources are assumed to impinge on a Uniform Linear Array (ULA) of N = 2 sensors. In addition, all the considered approaches are initialized using the truncated HOSVD procedure [START_REF] De Lathauwer | A multilinear singular value decomposition[END_REF]. Figure 3 shows, on the one hand, the good contribution of the ELS procedure to the ALS-SIBI UNS method and on the other hand the interest for the ELSALS-SIBI UNS algorithm in using symmetries of H x in the way we propose.

V. IDENTIFIABILITY

An identifiability study of the SIBI family is proposed hereafter under assumptions 1-3. In particular, the M time lag vectors τ (m) are chosen such that a well-conditioned matrix B is obtained. Furthermore, we have to select at least an M = N 2q time lag vectors. Note that, in practice, the choose of these time lag vectors relies on some prior information on sources. Our identifiability study is fourfold. First from a BUMI viewpoint, it gives the maximum number, P max , of sources which can be processed by the SIBI family. Secondly, from a numerical analysis viewpoint, it tackles the identifiability of the SAD, ACDC and ELSALS SYM procedures when a 3-way array with special symmetries is considered. Next, from a multilinear algebra viewpoint it gives some insights on the uniqueness of a joint CAND of several hermitian multiway arrays. Finally, from the VA theory viewpoint, we extend the VA concept [START_REF] Chevalier | On the virtual array concept for higher order array processing[END_REF] to the case of a combination of several VA's. Recall that the VA theory relies on the existence of N 2q sensors which virtually appear when using the 2q-th order statistical matrix C 2q,x [0] [START_REF] Chevalier | On the virtual array concept for higher order array processing[END_REF]. Regardless the application area, N 2q is equal to the maximal rank, rk max , of the (N q × P ) matrix A (q) :

N 2q = rk max (A (q) ) = rk max (A ( +1) A * ) (26)
where , A q denote the integer part of q/2 and the Khatri-Rao product of q matrices A, respectively. Some values of N 2q in radiocommunications context are given in [START_REF] Chevalier | On the virtual array concept for higher order array processing[END_REF] for several kinds of antennas. For instance, when a ULA of N sensors is used, N 2q is given by N 2q = q(N -1) + 1.

Let's now consider the SAD-SIBI JAD method. P max is the maximum value of P for which assumptions 1 and 2 are true. Then we have:

P max = rk max (A (3) ) = N 6 (27)
Regarding the five iterative ALS-SIBI approaches, their identifiability study consists in looking for the rank maximal of the 3-way array H x (14) they decompose. This should be achieved provided that the CAND of H x is unique. Finding this rank is usually a matter of trial and error [START_REF] Lathauwer | A link between the canonical decomposition in multilinear algebra and simultaneous matrix diagonalization[END_REF]. Nevertheless, under assumption 2 and according to the results of [START_REF] Jiang | Kruskal's permutation lemma and the identification of CANDECOMP/PARAFAC and bilinear models with constant modulus constraints[END_REF] we have:

P max = rk max (H (3) x ) = rk max (B(A (3) A (3) * ) T ) (28)
where (M ×N 6 ) H (3) x denotes the unfolding matrix associated to the third direction of H x [START_REF] Comon | Blind identification of underdetermined mixtures based on the characteristic function[END_REF]. For M ≥ N 6 , the maximum rank of H (3) x is generically equal to N 6 . Recall that a property is well-called generic when it holds everywhere except for a set of Lebesgue measure 0. However, in some applications such as in radiocommunication contexts, matrix H (3) x may belong to this set. In this case, the previous property does not hold generally. Intrinsically, we have: rk max (H (3) x ) ≤ rk max (A (6)

) = N 12 (29)
where A (6) = A (3) A (3) * . But, since according to [START_REF] Chevalier | On the virtual array concept for higher order array processing[END_REF] N 12 is strictly lower than N 6 , the number maximal of sources which can be processed using the ALS-SIBI approaches is then given by: P max = rk max (A (6) ) = N 12 [START_REF] Lathauwer | Fourthorder cumulant-based blind identification of underdetermined mixtures[END_REF] So from a BUMI viewpoint and according to equations ( 27) and [START_REF] Lathauwer | Fourthorder cumulant-based blind identification of underdetermined mixtures[END_REF] the SIBI methods can process more than N 4 sources. Note that N 4 is the maximum number of sources the classical BUMI cumulant-based methods [START_REF] Albera | Blind Identification of Overcomplete Mixtures of sources (BIOME)[END_REF], [START_REF]Blind identification and source separation in 2x3 underdetermined mixtures[END_REF], [START_REF] Estevao | Blind channel identification algorithms based on the parafac decomposition of cumulant tensors: the single and multiuser cases[END_REF], [START_REF] Ferreol | Fourth order blind identification of underdetermined mixtures of sources (FOBIUM)[END_REF], [START_REF] Karfoul | Canonical decomposition of even higher order cumulant arrays for blind underdetermined mixture identification[END_REF], [START_REF] Lathauwer | Blind identification of underdetermined mixtures by simultaneous matrix diagonalization[END_REF], [START_REF] Lathauwer | Fourthorder cumulant-based blind identification of underdetermined mixtures[END_REF] can process. Table I gives the number P max as a function of N for a ULA and for several cumulantbased methods. From a numerical analysis viewpoint, we show that: i) the SAD scheme is limited to 3-way arrays whose rank is lower than their smallest dimension. ii) Both the ACDC and ELSALS SYM procedures succeed generically in performing the CAND of a rank-N 2 (N×N×M ) array if M ≥ N 2 . Consequently, a more relaxed uniqueness condition can be obtained when a joint CAND of several hermitian multi-way arrays is performed. Since according to [START_REF] Chevalier | On the virtual array concept for higher order array processing[END_REF] N 6 = rk max (C 6, x [0]), equation [START_REF] Lathauwer | Fourthorder cumulant-based blind identification of underdetermined mixtures[END_REF] shows that a maximal rank of N 12 (> N 6 ) can be achieved when C 6, x [0] and at least N 12 -1 other arrays C 6, x [τ ] with the same loading matrices are simultaneously considered. This result extends the virtual array theory [START_REF] Chevalier | On the virtual array concept for higher order array processing[END_REF] to the case where several 2q-th order cumulant array are simultaneously used. Hence, a twofold increase of the virtual array aperture.

VI. NUMERICAL COMPLEXITY

The numerical complexity of the SIBI family is addressed hereafter. It is computed and compared to the one of classical BUMI methods, namely SOBI [START_REF] Belouchrani | A blind source separation technique using secondorder statistics[END_REF], JADE [START_REF] Cardoso | Blind beamforming for nongaussian signals[END_REF], FastICA [20, chapter 8], COM2 [START_REF] Comon | Independent Component Analysis, a new concept ?[END_REF], FOOBI1 [START_REF] Lathauwer | Fourthorder cumulant-based blind identification of underdetermined mixtures[END_REF], FOOBI2 [START_REF] Lathauwer | Fourthorder cumulant-based blind identification of underdetermined mixtures[END_REF], FOBIUM [START_REF] Ferreol | Fourth order blind identification of underdetermined mixtures of sources (FOBIUM)[END_REF] and 6-BIOME [START_REF] Albera | Blind Identification of Overcomplete Mixtures of sources (BIOME)[END_REF]. More particularly, it is computed in terms of number of floating point (complex) operations (flops) required to identify the (N×P ) mixing matrix A from K data snapshots. Note that a complex flop is defined as the sum of a complex multiplication and a complex addition. In practice, only multiplications are counted which does not affect the magnitude of the numerical complexity. Table II shows the computational complexity of the aforementioned methods where f 4 (N ) = O[N 4 /8] and f 6 (N ) = O[N 6 /72] denote the number of free entries of C 4, x and C 6, x , respectively. In addition, D and I denote the desired resolution in the unit disk for the COM2 approach and the number of sweeps used in the joint diagonalization algorithms, respectively. J, J 1 , J 2 and J 3 are the maximal numbers of iterations used by FastICA, ELSALS-SIBI UNS , ELSALS-SIBI SYM and ALS-SIBI ACDC , respectively. J 4 is the maximal number of iterations used for the estimation of the relaxation factor in ELSALS-SIBI UNS and ELSALS-SIBI SYM . Eventually, M represents the number of delay lags used in SOBI, FOBIUM and the SIBI techniques, M 1 = max{N 6 , P }, m 1 = min{N 6 , P }, M 2 = max{M, P } and m 2 = min{M, P }.

Figure 4 shows the minimum numerical complexity of the SIBI family and the aforementioned classical BUMI methods, as a function of the number P of sources. Indeed, for each value of P the minimum number N min of sensors of a ULA is computed such that the identifiability condition of each method is still valid. Hence, N min = P for SOBI, JADE, FastICA, COM2 and N min ≤ P ≤ q(N min -1)+1 for the other algorithms allowing for an identification of underdetermined mixtures. Note that the parameter q is directly related to the identifiability condition of each method. For example, q = 2 for FOOBI1, FOOBI2, FOBIUM and 6-BIOME, whereas q = 3 for the SAD-SIBI JAD approach and q = 6 for both the ALS- SIBI ACDC and the ALS-SIBI techniques. M , I and D are chosen to be equal to 9, 1+floor{P 1/2 } and 0.05, respectively. As far as the iterative approaches are concerned, J and J 1 are chosen such that a good identification of the mixing matrix is guaranteed, hence J = 1000, J 1 = 800, J 2 = 400, J 3 = 600. On the other hand, we put J 4 = 20, which ensures a good estimation of the relaxation factor. Note that the latter three parameters are highly related to many factors such as SNR, direction of arrivals, number of sources, number of samples, etc. As depicted in figure 4, the ELSALS-SIBI approaches belong to the group of the more expensive methods such as FastICA whereas the SAD-SIBI technique is less expensive as 6-BIOME, FOBIUM, COM2 FOOBI1.

VII. SIMULATION RESULTS

Five performance evaluation studies of the SIBI family are presented hereafter, and more precisely the performance of SAD-SIBI JAD , ALS-SIBI ACDC and ELSALS-SIBI SYM for both underdetermined and overdetermined mixtures of sources. The performance analysis was achieved using the criterion D(A, A) = (α 1 , α 2 , • • • , α P ), which measures a pseudo-distance between mixture A and its estimate A where α p = min 1≤i≤P d(a p , a i ) with d the pseudo-distance between vectors defined by [START_REF] Ferreol | Fourth order blind identification of underdetermined mixtures of sources (FOBIUM)[END_REF]:

d(u, v) = 1 -u H v 2 / u 2 v 2 (31)
It is noteworthy that only criterion α 1 is displayed in figures 5,7, 8 and 9 since the P -1 other criteria α p have shown identical results. The P sources considered in this section are Quadrature Phase Shift Keying (QPSK) linearly modulated with pulse shape filter corresponding to a 1/2-Nyquist filter with a roll-off of 0.3. The latters are assumed to impinge on a ULA of N sensors and are chosen to be well angularly separated.

Both first studies consist in evaluating the performance of the SIBI family in a strongly underdetermined mixture which of both considered methods. It appears that both methods are effectively able to process up to seven sources from a ULA of two sensors. In addition, in comparison with figure 5, ELSALS-SIBI SYM seems to be less sensitive to a strong increase of the number P of sources even it also needs more samples to give a constant performance as P increases.

The third study consists in comparing, in a more sympathetic underdetermined mixture case, the SIBI family and cumulant-based methods such as FOBIUM [START_REF] Ferreol | Fourth order blind identification of underdetermined mixtures of sources (FOBIUM)[END_REF], 6-BIOME [START_REF] Albera | Blind Identification of Overcomplete Mixtures of sources (BIOME)[END_REF], as a function of the number of samples. For this purpose N = 2 sensors and P = 3 sources are considered with residual carriers given by f c1 T e = 0, f c2 T e = 0.35, f c3 T e = 0.7 and an SNR of 15 dB. While FOBIUM and 6-BIOME appear to be efficient in figure 7, especially when the number of samples increases, the SIBI family shows a higher performance even for a small number of samples. Now we evaluate the influence of SNR on the performance of the SAD-SIBI JAD , ALS-SIBI ACDC , ELSALS-SIBI SYM , FOBIUM and 6-BIOME using the same application context as the one used in figure 7. Under the assumptions, figure 8 shows that the five methods have a similar behavior for small values of SNR. Nevertheless, the SIBI techniques give better results as soon as SNR increases beyond 0 dB. Eventually, figure 9 shows the variations of α 1 as a function of the spatial correlation factor of a Gaussian noise at the output of the three considered SIBI approaches and seven classical cumulant-based methods such as COM2 [START_REF] Comon | Independent Component Analysis, a new concept ?[END_REF], Fas-tICA [20, chapter 8], FOOBI1, SOBI [START_REF] Belouchrani | A blind source separation technique using secondorder statistics[END_REF], JADE [START_REF] Cardoso | Blind beamforming for nongaussian signals[END_REF], 6-BIOME and FOBIUM. For this purpose, we considered 1000 samples, N = 3 sensors and P = 2 sources with residual carriers such that f c1 T e = 0, f c2 T e = 0.65 and a SNR of 5 dB for all sources. Note that the Gaussian noise model used in this simulation is the sum of an internal noise ν in [k] and an external noise ν out [k] of covariance matrices C in 2,ν [k, 0] and C out 2,ν [k, 0], respectively, such that:

C in 2,ν [0] r,q = σ 2 2 δ(r-q) C out 2,ν [0] r,q = σ 2 2 ρ |r-q| (32)
where σ 2 , ρ, C 2,ν [0] r,q = C in 2,ν [0] r,q + C out 2,ν [0] r,q are the total noise variance per sensor, the noise spatial correlation factor and the (r, q)-th component of the total noise covariance matrix, respectively. The computer results show on the one hand how the performance of the SIBI family is stable, even for a strong noise spatial correlation. In one hand, the classical methods seem to be less robust, especially COM2, FastICA, JADE and SOBI which require a prior spatial whitening based on SO cumulants.

VIII. CONCLUSION

We proposed a new class of BUMI methods based on the joint CAND of a set of HO cumulant arrays. This class of algorithms is described through SixO cumulants, giving rise to two SixO subfamilies of techniques, namely the SAD-SIBI and the ALS-SIBI methods. More particularly, these two families perform a CAND of a 3-way array built from a set of SixO cumulant arrays but in a different way. Indeed, the former approach is algebraic while the latter is iterative. All approaches were theoretically studied especially in terms of identifiability and computational complexity. In particular, the identifiability analysis showed that the SIBI algorithms can process more sources than the classical cumulant-based methods. In addition, computer results showed the good behavior of the SIBI family against classical methods, for both overdetermined and underdetermined mixtures of sources. Besides, from a numerical analysis viewpoint, we deal with the identifiability and the convergence speed of several 3-way array decomposition procedures such as the symmetric ACDC scheme. We also show how to accelerate the iterative CAND of a 3-way array, which has special symmetries. Next, from a multi-linear algebra viewpoint we give some insights on the uniqueness of a joint CAND of several hermitian multi-way arrays compared to the CAND of only one array. In addition, this allows us to extend the concept of Virtual Array (VA) to the case of a combination of several VA's. Indeed, the use of several VA's with the same number of virtual sensors amounts to use a VA with a higher aperture.

APPENDIX A PROOF OF PROPOSITION 5

The full column rank of B implies that there is at least one M -dimensional vector β with values in the complex field such that:

B T β = [1, 1, • • • , 1] T ∈ C P (33)
According to the definition of matrix B given in proposition 2, (33) is equivalent to the following equation: where β m is the m-th component of the M -dimensional column vector β and I P is the (P ×P ) identity matrix. Then, under assumption 1, the following matrix:

  3)} = {{1}, {2}, {3}}. with J(1, 3) = J(3, 3) = 1 and J(2, 3) = 3 different partitions.

Fig. 1 .

 1 Fig. 1. CAND of a rank-P 3-way array

Fig. 2 .

 2 Fig. 2. Unfolding matrix of a 3-way array with respect to the first mode, given by concatenating the vertical slices of the original 3-way array.

Fig. 3 .

 3 Fig. 3. Loss function Ψ (1) 3 at the output of the ALS-SIBI methods for P = 3 poorly angularly separated sources impinging on a ULA of N = 2 sensors.

N

Fig. 4 .

 4 Fig. 4. Minimal comlexity of the SIBI family and other BUMI methods as a function of the number of sources, for a ULA and 1000 data snapshots.

Fig. 7 .

 7 Fig. 7. α 1 for a ULA of 2 sensors and 3 QPSK's of SNR= 15 dB.

Fig. 8 .

 8 Fig. 8. α 1 for a ULA of 2 sensors and 3 QPSK's of 1000 samples.

Fig. 9 .

 9 Fig. 9. α 1 as a function of the spatial correlation factor of a Gaussian noise for a ULA of 3 sensors and 2 QPSK's of 1000 samples.

β

  m ζ 6,s [τ (m) ] = I P[START_REF] Moreau | A generalization of joint-diagonalization criteria for source separation[END_REF]

β

  m C 6,x [τ (m) ] = A (3) ( M m=1 β m ζ 6,s [τ (m) ])A (3)H (35)is a rank-P positive semidefinite matrix equal to A (3) A (3)H . Hence the result.

Computational complexity

N : number of sensors, P : number of sources, M : number of time lags, L: number of data samples, I: number of sweeps, J, J1, J2, J3, J4: number of iterations in FastICA, in the ELSALS-SIBIUNS, ELSALS-SIBISYM in the estimation of the relaxation factor in ELSALS-SIBI and in the ALS-SIBIACDC, respectively COM2 min{LN 2 /2 + 4N 3 /3 + P N L, 2LN 2 } + min{12If 4 (P )P 2 + 2IP cannot be processed by any classical cumulant method based on at most order 6, say P ≥ 4 sources from a ULA of only N = 2 sensors. The sources have the same Signal-to-Noise Ratio (SNR) equal to 15 dB. On the one hand, figure 5 shows the variation of α 1 at the output of the three SIBI methods as a function of the number of samples for P = 4 sources with carrier residuals such that f c1 T e = 0, f c2 T e = 0.25, f c3 T e = 0.5, f c4 T e = 0.75. In addition, the performance was compared for two numbers M of sets τ (m) of delays, say M = 9 and M = 27. Figure 5 shows a good behavior of the SIBI family. Besides, the semi-algebraic SAD-SIBI JAD method seems to be insensitive to the increase in the number M of used sets of delays, contrary to the iterative ALS-SIBI ACDC and ELSALS- SIBI SYM algorithms. As shown in section V, more statistical information may ensure the uniqueness of CAND and then a good quality of factor's identification. Since better results are obtained for M = 27 sets of delays, the latter were used in the following studies. On the other hand, the case of P = 7 sources with residual carriers given by f c1 T e = 0, f c2 T e = 0.15, f c3 T e = 0.3, f c4 T e = 0.45, f c5 T e = 0.6, f c6 T e = 0.75, and f c7 T e = 0.9, is considered in figure 6. Note that, in such a practical context, only ALS-SIBI ACDC and ELSALS-SIBI SYM are able to process as many sources, as mentioned in table I. More precisely, figure 6 shows the maximum α p , named α max , for each data sample at the output of ALS-SIBI ACDC and ELSALS-SIBI SYM , say the worst performance over all sources