Bahram Ehsandoust
email: bahram.ehsandoust@gipsa-lab

Bertrand Rivet
email: bertrand.rivet@phelma

Christian Jutten
email: christian.jutten@gipsa-lab.grenoble-inp.fr

Massoud Babaie-Zadeh

Nonlinear Blind Source Separation for Sparse Sources

Keywords: Blind Source Separation, Independent Component Analysis, Sparse Signals, Manifold Learning

Blind Source Separation (BSS) is the problem of separating signals which are mixed through an unknown function from a number of observations, without any information about the mixing model. Although it has been mathematically proven that the separation can be done when the mixture is linear, there is not any proof for the separability of nonlinearly mixed signals.

Our contribution in this paper is performing nonlinear BSS for sparse sources. It is shown in this case, sources are separable even if the problem is under-determined (the number of observations is less that the number of source signals). However in the most general case (when the nonlinear mixing model can be of any kind and there is no side-information about that), an unknown nonlinear transformation of each source is reconstructed. It is shown why the problem reconstructing the exact sources is severely ill-posed and impossible to do without any other information.

I. INTRODUCTION

In signal processing applications, it is usual to have a number of signals measured by some sensors, while each of them might be a mixture of a number of source signals. BSS is the problem of separating the original sources from the observed mixture, without any information about either the source signals or the mixing model. The basic model of the mixing model is depicted in Fig. 1, and can be expressed as

x(t) = f (s(t)) (1)
where f : R N → R M is the unknown mixing function, N and M are the number of source and observation signals respectively and s(t) = [s 1 (t), s 2 (t), . . . , s N (t)] T and x(t) = [x 1 (t), x 2 (t), . . . , x M (t)] T (T denotes the transpose of the vector) are the source and observation vectors respectively.

The BSS problem can also be written as find g s.t. y(t) = g(x(t)) = g(f (s(t))) = ŝ(t)

where g : R M → R N is the separating function which is to be estimated and y(t) = [y 1 (t), y 2 (t), . . . , y N (t)] T is the vector of reconstructed sources (the output of the BSS algorithm).

It should be noted that the goal of BSS is to "separate" the sources and not to "reconstruct" them. Depending on the mixing model, there are some ambiguities in reconstructing the sources that can only be resolved using some prior information about the source signals, which is out of the scope of BSS Fig. 1. BSS problem basic model (for example in case of linear mixtures, separated sources would be equal to the original ones up to a scaling factor and reordering).

BSS problem has been studied since the early works in 80s [START_REF] Hérault | Space or time adaptive signal processing by neural networks models[END_REF], and noticeable results have been achieved for some cases. The general idea of performing BSS is to design a separating function (g) such that each entry of the output signals y(t) only depends on a single original source. For example, it has been shown [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF] that if source signals are mutually independent, the mixing function (f ) is linear and the number of the sources is equal to the number of the observations (x(t)), one can perform the separation by making the outputs mutually independent as well as the sources (which is called ICA1 ). In this case, the "independence" is the property of the sources, which is tried to be recovered for the separation. There are some other works exploiting other properties to perform the separation (a comprehensive survey can be found in [START_REF] Comon | Handbook of Blind Source Separation: Independent component analysis and applications[END_REF]), such as non-negativity [START_REF] Cichocki | New algorithms for nonnegative matrix factorization in applications to blind source separation[END_REF], non-Gaussianity [START_REF] Comon | Independent component analysis, a new concept?[END_REF], nonstationarity [START_REF] Parra | Convolutive blind separation of non-stationary sources[END_REF] and non-whiteness [START_REF] Buchner | Blind source separation for convolutive mixtures exploiting nongaussianity, nonwhiteness, and nonstationarity[END_REF] (they are named as the "non-properties" [START_REF] Mei | Blind source separation based on cumulants with time and frequency non-properties[END_REF]).

Although there are many works and proposed algorithms on BSS for linear mixtures, nonlinear BSS has not been well developed ( [START_REF] Deville | An overview of blind source separation methods for linear-quadratic and post-nonlinear mixtures[END_REF] provides a short survey on nonlinear BSS). It is shown through counter examples that ICA does not separate the sources which are non-linearly mixed [START_REF] Hosseini | On the separability of nonlinear mixtures of temporally correlated sources[END_REF], [START_REF] Babaie-Zadeh | On blind source separation in convolutive and nonlinear mixtures[END_REF]. Therefore, except a few general approaches which try to investigate the separability of any kind of nonlinear mixture models [START_REF] Ehsandoust | Blind source separation in nonlinear mixture for colored sources using signal derivatives[END_REF], [START_REF] Levin | Performing nonlinear blind source separation with signal invariants[END_REF], studies on nonlinear BSS are focused on specific mixture model and source signals, out of which we can name Post Nonlinear (PNL) [START_REF] Taleb | Source separation in post nonlinear mixtures[END_REF] and Bi-Linear (or Linear Quadratic) mixtures [START_REF] Merrikh-Bayat | Linear-quadratic blind source separating structure for removing show-through in scanned documents[END_REF], Convolutive Post Nonlinear mixtures [START_REF] Babaie-Zadeh | On blind source separation in convolutive and nonlinear mixtures[END_REF] and conformal mappings [START_REF] Hyvärinen | Nonlinear independent component analysis: Existence and uniqueness results[END_REF].

In this paper, we consider the problem where the mixing model is nonlinear and the sources are sparse (the case of sparse sources has already been studied for linear mixtures in many works like [START_REF] Bofill | Underdetermined blind source separation using sparse representations[END_REF]). So the mixing model that is studied in this work is not limited to a specific kind and can be any function but invertible. We show that in this case, the sources can be separated up to a nonlinear function.

The rest of the paper is organized as follows. In the next section, the main idea is introduced and the approach and the algorithm for performing the separation are proposed. Then some simulation results are shown in the section III. A comprehensive discussion on the performance of the proposed approach and how to develop it for the future works is finally made in section IV.

II. THE PROPOSED ALGORITHM

Sparse signals in the time domain are the signals which most of the time take the value zero. In other words, the probability of their value being equal to zero is high. Therefore, if all of the sources are sparse, it will be even more rare that they take non-zero values at the same time. This is the main feature that is exploited in this approach to perform the separation. The proposed algorithm comprises two main steps of 1) Clustering the Observations and 2) Unfolding the Manifold, which are described as follows.

A. Step 1 -Clustering the Observations

Let us assume a 2 × 2 mixing system (2 source signals and 2 observations) as an example and

x 1 (t) = e s1(t) -e s2(t) (3)
x 2 (t) = e -s1(t) + e -s2(t) (4)

as the nonlinear mixing function (x 1 (t) and x 2 (t) are then centered). For investigating what would happen in the case of sparse sources from a geometric point of view, the scatter plots of the observations for the two cases whether the sources are sparse or not, are plotted in Fig. 2. As it is already mentioned, when the sources are sparse, the samples of the source vector are mainly concentrated around the axes because it is quite rare that both of the sources take a non-zero value at the same time. So in this case, the scatter plot of the observations (Fig. 2) contains two manifolds (each of which is the result of the transformation of one of the axes in the source space). As a consequence, if we consider the points where two sources are active (i.e. nonzero) at the same time as outliers, we can apparently distinguish the two manifolds in the scatter plot of the observation vector. This phenomenon can be also expressed mathematically as follows. let us rewrite (1) as

x 1 (t) x 2 (t) = f ( s 1 (t) s 2 (t) ) = f 1 (s 1 (t), s 2 (t)) f 2 (s 1 (t), s 2 (t)) (5)
where f 1 (•) and f 2 (•) are the components of the mixing matrix. For the times when only s 1 is active we have

x 1 (t) = f 1 (s 1 (t), 0) = h 11 (s 1 (t)) x 2 (t) = f 2 (s 1 (t), 0) = h 12 (s 1 (t)) (6)
Fig. 2. Comparing scatter plots of the source and observation vectors whether the sources are sparse or not and when only s 2 is active we have

x 1 (t) = f 1 (0, s 2 (t)) = h 21 (s 2 (t)) x 2 (t) = f 2 (0, s 2 (t)) = h 22 (s 2 (t)) (7)
where

h 11 (•) and h 21 (•) are R → R marginal functions of f 1 (•, •), and h 12 (•) and h 22 (•) are R → R marginal functions of f 2 (•, •).
Since when s 1 (t) is active, both x 1 (t) and x 2 (t) are functions of only one source (s 1 (t)), they would lie on a 1dimensional manifold in two dimensional space, and similarly for the samples when only s 2 (t) is active. In other words, when only s 1 (t) is active and at least one of h 11 and h 12 is invertible (we assume h 11 is invertible), we would have

x 2 (t) = h 12 (s 1 (t)) = h 12 (h -1 11 (x 1 (t))) = q 1 (x 1 (t)) (8
) and similarly when only s 2 (t) is active and assuming h 21 is invertible), we would have

x 2 (t) = h 22 (s 2 (t)) = h 22 (h -1 21 (x 1 (t))) = q 2 (x 1 (t)
). (9) Thus x 2 (t) = q 1 (x 1 (t)) and x 2 (t) = q 2 (x 1 (t)) would be the two manifolds that can be clustered in the observation space via some clustering (manifold learning) algorithms. The samples which lie on the manifold of q 1 , are assigned to s 1 (t) and the ones that lie on the manifold of q 2 , are assigned to s 2 (t).

It should be noted that although the above explanation was based on a simple 2 × 2 mixing system, the argument is still valid for higher dimensions as well. Because for any number of sources and observations, for sparse sources, only one of the sources is active in most of the times. Therefore, each axis in the N-dimensional source space transforms to a 1-dimensional manifold in the M-dimensional observation space.

B. Step 2 -Unfolding the Manifold

As mentioned before, each manifold in the observation space is a nonlinear transformation of one of the axes in the source space. The problem of estimating the inverse of this transformation is severely ill-posed. This can also be shown though mathematical expressions. From ( 6) and [START_REF] Mei | Blind source separation based on cumulants with time and frequency non-properties[END_REF] we have

s 1 (t) =      h -1 11 (x 1 (t)) = h -1 12 (x 2 (t)) x(t) ∈ q 1 0 x(t) ∈ q 2 Unknown o.w. (outlier) (10)
s 2 (t) =      0 x(t) ∈ q 1 h -1 21 (x 1 (t)) = h -1 22 (x 2 (t)) x(t) ∈ q 2 Unknown o.w. (outlier) (11)
while we have information about none of h ij (•) functions. The "unknown" values mentioned in [START_REF] Babaie-Zadeh | On blind source separation in convolutive and nonlinear mixtures[END_REF] and ( 11) are because, generally, we do not have enough information to separate the outliers (when more than one source are simultaneously active) and they actually cause the error in separation. More discussion on this can be found in the following.

In other words, it let u 1 (•) and u 2 (•) be two arbitrary R → R invertible functions. So, if instead of s 1 (t) and s 2 (t), u 1 (s 1 (t)) and u 2 (s 2 (t)) are mixed through the function f (•), ( 6) and ( 7) would become

x 1 (t) = h 11 (u 1 (s 1 (t))) = ĥ11 (s 1 (t)) x 2 (t) = h 12 (u 1 (s 1 (t))) = ĥ12 (s 1 (t)) (12)
(when only s 1 is active) and

x 1 (t) = h 21 (u 2 (s 2 (t))) = ĥ21 (s 2 (t)) x 2 (t) = h 22 (u 2 (s 2 (t))) = ĥ22 (s 2 (t)) (13)
(when only s 2 is active) where ĥij = h ij • u i for i, j ∈ {1, 2} are some other nonlinear functions of the original sources. Consequently, considering [START_REF] Deville | An overview of blind source separation methods for linear-quadratic and post-nonlinear mixtures[END_REF] and ( 9) we can say

x 2 (t) = ĥ12 (s 1 (t)) = ĥ12 • ĥ-1 11 (x 1 (t)) = h 12 • u 1 • u -1 1 • h -1 11 (x 1 (t)) = h 12 • h -1 11 (x 1 (t)) = q 1 (x 1 (t)) (14)
(when only s 1 is active) and

x 2 (t) = ĥ22 (s 2 (t)) = ĥ22 • ĥ-1 21 (x 1 (t)) = h 22 • u 2 • u -1 2 • h -1 21 (x 1 (t)) = h 22 • h -1 21 (x 1 (t)) = q 2 (x 1 (t)) (15)
(when only s 2 is active), which means that the q 1 and q 2 manifolds in the observation space remain the same. That is, the scatter plot of the observation vector does not contain any information about the probable nonlinear transformation of the sources. Therefore, the ambiguity of a nonlinear function in reconstructing the sources cannot be resolved. It is worth noting again that although the above description is made based on a simple 2 × 2 example, all arguments are still valid for the general case. In this case, each of the N sparse sources (i.e. each of the N axes in the source space) transforms to a 1-dimensional manifold in the M-dimensional observation space (which are called as q 1 , q 2 , . . . and q N ).

It should also be emphasized that the proposed algorithm works even for under-determined cases (where the number of observations is less than the number of source signals). This outstanding capability comes from the fact that the separation of sparse sources is based on separating N 1-dimensional manifolds, and even in a 2-dimensional space we can have any number of different 1-dimensional manifolds. Therefore, the minimum required number of the observation signals, despite of the number of sources, is always 2.

As a result, each of the observations can be taken as the separated sources. So the separated sources in the general form can be reconstructed as

ŷi (t) =      x j (t) x(t) ∈ q i 0 x(t) ∈ q k (k = i) Unknown o.w. (outlier) (16)
where ŷ(t) = [ŷ 1 (t), ŷ2 (t), . . . , ŷN (t)] T is the vector of the separated sources (here we use (•) because it is not the final results; after proposing a method for solving the problem of outliers, the final output is formulated in the following), i is the index of the separated source which is between 1 and N , k is any index between 1 and N , and 1 ≤ j ≤ M is the index of one of the observations (which is chosen arbitrarily).

As mentioned earlier, values of the separated signals for outliers (when more than one source is active), is "unknown", but considering some additional information about either the sources or the mixing model, leads to an estimation for them. As explained before, the proposed approach is based on the clustering of the manifolds in the observation space (each of which corresponds to one of the sources). Thus, the observation vector lies on a manifold when only one of the sources is active. When two or more sources are active, the observation sample is assumed to be outlier (in the manifold learning process) and is not assigned to any of the sources. This is why the sources can not be estimated in these cases and their values are said to be "unknown" in [START_REF] Babaie-Zadeh | On blind source separation in convolutive and nonlinear mixtures[END_REF], [START_REF] Ehsandoust | Blind source separation in nonlinear mixture for colored sources using signal derivatives[END_REF] and [START_REF] Bofill | Underdetermined blind source separation using sparse representations[END_REF]. Different manners can be suggested to estimate the value of each source for outliers depending on probable prior information about either the sources or the mixing function. The manners can be classified in the two following categories.

The first category of methods for estimating the sources in case of outlier observations, is based on using the sources themselves. In this case, we face the problem of estimating some missing values of a signal while we know it for most of the time. Designing the best restoration technique, one has to employ some known features of the signal (e.g. its bandwidth) if there is any. So the separated sources in this approach would be as

y i (t) =      x j (t) x(t) ∈ q i 0 x(t) ∈ q k (k = i) φ(ŷ i (t)) o.w. (outlier) (17)
where φ(•) is the chosen restoration function and ŷi (t) comes from [START_REF] Bofill | Underdetermined blind source separation using sparse representations[END_REF].

The other class of methods for doing the estimation is based on the mixing function and the separated manifolds. Using this approach, the value of each source can be estimated through nonlinearly projecting the outliers on the corresponding manifold in the observation space. However, the nonlinear projection is not unique, and would be more accurate if there is some side-information about the mixing model. Nevertheless, the separated sources in this approach can be written as

y i (t) =      x j (t) x(t) ∈ q i 0 x(t) ∈ q k (k = i) x j (t) o.w. (outlier) (18)
where x j (t) is the chosen nonlinear projection of x(t) onto the j th coordinate (axis of e j ). However, if there is prior information about the sources and/or the mixing model, it can be utilized to improve the performance of the separating algorithm and trying to reconstruct the sources (resolving the ambiguity of the nonlinear function on each source). For example if we know that the transformation preserves the local angles, if its time-derivative is continuous, if it has a parametric model, and so forth, this may lead to a much better estimation of the sources signals (which is left for the future works).

Finally, it is worth noting that although choosing one of the observations for reconstructing the separated sources is the easiest way, one may find other some wiser ways to make the separated signals. For example if the projection of a manifold (q i ) on one of the axes (x j ) is not a invertible transformation, that axis (x j ) is not a good choice for reconstructing the separated signals (e.g. in Fig. 3, projecting the more complicated manifold on both of the axes would not be invertible and would loose information). Because in this case, the reconstructed signal, the unresolved nonlinear transformation of the source is not invertible which looses some information of the source.

Therefore, it is better to take an unfolded version of each manifold as the separated source which is made by a dimension reduction algorithm (e.g. ISOMAP [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF]). Moreover, one may suggest other methods for unfolding the manifold considering the nature of the sources and the mixing model. For example assuming that the nonlinear mixing function preserves distances or angles, leads to specific unfolding (nonlinear dimension reduction) techniques.

III. SIMULATION RESULTS

In order to simulate the proposed algorithm, we have used several nonlinear mixing models and we have applied them on two sparse sources (with the sparsity of 75%) that are uniformly distributed in [-0.5, 0.5] when they are active.

In the simulation, firstly outliers are omitted (they are detected by computing their average distance their few neighbors). Then a simple curve fitting algorithm is applied to perform clustering. It is obvious that in more complicated examples (like Fig. 3), one should employ more powerful clustering techniques and might use kernels.

In each of the figures from Fig. 4 to 6, the scatter plot of the observations and the clustering result are shown in the above plots (where each cluster is indicated with a different color and the outliers are marked as "+"). The two below plots of each figure represent the relation between each of outputs and the corresponding source signal (which is, as expected, Fig. 6. The result of the algorithm for the model x 1 (t) = sin(2s 1 (t)-s 2 (t)) and x 2 (t) = sin(s 2 (t) -s 1 (t)) a nonlinear function that can not be resolved without further information).

IV. DISCUSSION AND FUTURE WORKS

In this work, we proposed a nonlinear BSS algorithm for sparse sources. The algorithm works even if the problem even if the number of observations is less that the number of sources (it only need more than one observation). We also showed that although the separation is done through the proposed method, the sources cannot be exactly reconstructed in the most general case without further information.

The proposed algorithm is also useful when the number of source signals is unknown. It is due to two facts: 1) despite of the number of sources, the separation can be performed as soon as we have more than one observation and 2) if the manifold learning technique which is used for clustering the observation is powerful enough (and does not need the number of clusters in advance), any number of manifolds can be separated in the observation space and hence the generalization of the proposed method would be straightforward.

In addition, the sparsity feature used in the proposed algorithm is satisfied even if each source signal takes a constant value for most of the time (not necessarily zero) and has sparse variations. in this case, all mentioned arguments remain valid but in ( 6) and ( 7), f j (s 1 (t), 0) and f j (0, s 2 (t)) need to be replace with f j (s 1 (t), c 2 ) and f j (c 1 , s 2 (t)) for j ∈ {1, 2}, where c 1 and c 2 are respectively the constant values that s 1 (t) and s 2 (t) are assumed to take most of the time.

For the future work, it would be interesting to develop the proposed framework for the sources that are not sparse in time domain, but in some other domain like frequency domain. In these cases, one has to firstly model the transformation to the sparse domain for both sources and observations in order to be able to cluster the observation from the manifolds in the sparse domain and then separate the observations using the classification which is done.

Moreover, it would be useful to apply the proposed approach to practical applications and to utilize the prior information (related to the real case) for better unfolding the separated manifolds and more precisely estimating the values corresponding to the outliers. It might also need a discussion on different choices of φ(•) in [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF] or the nonlinear projection method in (18).

Fig. 3 .

 3 Fig. 3. A manifold whose projections on the axes are not invertible

Fig. 4 . 2 Fig. 5 .

 425 Fig. 4. The result of the algorithm for the model x 1 (t) = cos(α(t))s 1 (t)sin(α(t))s2(t) and x 2 (t) = sin(α(t))s 1 (t)+cos(α(t))s2(t) where α(t) = π 2 (1 -s 2 1 (t) + s 2 2 (t)) 2

Independent Component Analysis

ACKNOWLEDGEMENT

This work is partly funded by the European project 2012-ERC-Adv-320684 CHESS.


Rafael Assato Ando

Christian Jutten

Bertrand Rivet

Romis Attux

Leonardo Duarte

Nonlinear Blind Source Separation for Chemical Sensor Arrays Based on a Polynomial Representation

In this paper we propose an extension of a blind source separation algorithm that can be used to process the data obtained by an array of ion-selective electrodes to measure the ionic activity of different ions in an aqueous solution. While the previous algorithm used a polynomial approximation of the mixing model and used mutual information as means of estimating the mixture coefficients, but it only worked for a constrained configuration of two sources with the same ionic valence. Our proposed method is able to generalize it to any number of sources and any type of ions, and is therefore able to solve the problem for any configuration. Simulations show good results for the analyzed application.

I. INTRODUCTION

The blind source separation (BSS) problem is an important cornerstone in signal processing theory [START_REF] Comon | Handbook of blind source separation[END_REF], and its applications include audio signal processing, telecommunications, image processing, brain-computer interface design, analysis of seismic data, among others. Essentially, the BSS problem consists of estimating a set of signals, also called sources, using measurements which are effectively a mixture of them, possibly including some a priori information on the nature of the desired signals and mixing process.

The standard approach for solving the BSS problem is based on the hypotheses that the mixing process is linear and that the source signals are statistically independent. In that case, successful algorithms such as the independent component analysis (ICA) [START_REF] Comon | Handbook of blind source separation[END_REF] have been proposed and extensively studied over the last years. However, for certain applications the linear model is insufficient, and a nonlinear mixing model might be required. Among the nonlinear BSS mixing models that have been studied so far, it is important to mention the postnonlinear (PNL) [START_REF] Achard | Criteria based on mutual information minimization for blind source separation in post nonlinear mixtures[END_REF], [START_REF] Duarte | Blind source separation of post-nonlinear mixtures using evolutionary computation and order statistics[END_REF], in which a linear mixture is further distorted by a nonlinear function before the measurement, as well as the linear quadratic (LQ) model, which represents the mixtures as a second-order polynomial on the sources.

In this work, we shall focus on the problem of source separation for ion-selective electrode (ISE) arrays [START_REF] Duarte | A dataset for the design of smart ion-selective electrode arrays for quantitative analysis[END_REF], [START_REF] Duarte | Application of blind source separation methods to ionselective electrode arrays in flow-injection analysis[END_REF], [START_REF] Duarte | Design of smart ion-selective electrode arrays based on source separation through nonlinear independent component analysis[END_REF], [START_REF] Ando | A blind source separation method for chemical sensor arrays based on a second order mixing model[END_REF]. ISEs are simple devices which are used for measuring the ionic The authors would like to thank CNPq, FAPESP, CNRS and ERC project 2012-ERC-AdG-320684 CHESS for the funding. activity (essentially the effective concentration of an ion) in an aqueous solution. An ISE consists of a sensitive membrane for which the electrochemical potential varies according to the concentration of a specific ion. A well-known example of ISE is the glass electrode used for measuring the pH value [START_REF] Mikhelson | Ion-selective electrodes[END_REF].

The problem of the ISEs, however, is that the membranes are not perfectly selective. When trying to estimate the activity of a certain ion, the measurement will include not only the activity from the target ion, but also an interference from different ions which are also present in the same chemical solution. Because the measurements are a mixture of the activities, a BSS paradigm can be used to separate the signals from each ion by using an array of multiple ISEs.

The problem of estimating the ionic activities based on ISE arrays has already been studied in the past, with the use of a Bayesian source separation method [START_REF] Duarte | A bayesian nonlinear source separation method for smart ion-selective electrode arrays[END_REF], sparsitybased methods [START_REF] Duarte | A sparsity-based method for blind compensation of a memoryless nonlinear distortion: Application to ion-selective electrodes[END_REF] and a PNL model and with a polynomial representation as well, but with limiting constraints which made it work only in specific circumstances [START_REF] Ando | A blind source separation method for chemical sensor arrays based on a second order mixing model[END_REF]. In this paper, we shall propose a more general version of the polynomial algorithm, which should be able to account for all possible configurations, at the expense of additional complexity.

In the paper, we will initially present, in Section II, a description of the BSS problem statement and how it relates to measurements with chemical sensor arrays. In Section III, we will briefly present the previous algorithm for the constrained version of the problem, followed by an explanation of the proposed method for its general version. We shall then present experimental results in Section IV, and a conclusion and future plans shall finally be presented in Section V.

II. PROBLEM STATEMENT

A. Blind Source Separation

For the problem of BSS, let us consider a set of N sources s T , where n is the temporal index and F[.] is the mixing function, not necessarily linear. The goal of the source separation problem is to obtain the sources s(n) given the mixtures x(n), and possibly some a priori information on the properties of the sources (e.g. independence, sparsity, non-negativity, etc) and/or on the mixing process (e.g. linearity, PNL model, etc). The problem is said to be non-blind if the mixing function F[.] is known in advance, and it is otherwise said to be blind. If the function F[.] is linear, the mixing function can be represented by a matrix:

(n) = [s 1 (n), s 2 (n), ..., s N (n)] T and M mixtures x(n) = F[s(n)] = [x 1 (n), x 2 (n), ..., x M (t)]
x(n) = As(n) ( 1
)
where A is the N ×M mixing matrix. For the linear case, there are several classical methods in the literature that can solve the BSS problem under different priors for the sources. A very well-known class of methods is the independent component analysis (ICA) [START_REF] Comon | Handbook of blind source separation[END_REF].

For the nonlinear case, however, no general solution exists, and the problem has to be analyzed on a case-by-case basis [START_REF] Comon | Handbook of blind source separation[END_REF], depending on the nature of the mixing function. Among the commonly studied nonlinear models are the post nonlinear (PNL) [START_REF] Duarte | Blind source separation of post-nonlinear mixtures using evolutionary computation and order statistics[END_REF] and the polynomial model [START_REF] Deville | Recurrent networks for separating extractable-target nonlinear mixtures. part i: Non-blind configurations[END_REF].

B. Chemical Sensor Arrays

For the problem of chemical sensor arrays, let us consider an aqueous solution in which we have N different ions and M ion selective electrodes (ISE) to measure their activities s i , which are our sources. However, due to the interference problem, the values measured are actually a nonlinear mixture of the activities of each ion. According to the Nicolsky-Eisenman (NE) equation [START_REF] Bakker | Selectivity of potentiometric ion sensors[END_REF], this interference can be approximately modeled by the equation below:

y i (n) = e i + RT z i F ln   s i (n) + N j=1,j =i a ij s zi/zj j (n)   (2)
where y i are the observed values (our mixtures), e i is the standard electrode potential (a scalar constant), R is the universal gas constant, T is the temperature, F is the Faraday constant, z i is the valence of the ith ion, and a ij are the selectivity coefficients. The RT ziF slope is called the Nernst slope, and at room temperature of 25 • C, it is approximately 26mV, or 59mV if the logarithm is converted to base 10 (which is the standard practice in chemical applications).

III. METHODS AND ALGORITHM

A. Previous Algorithm

In the simplest version of the problem, we have the case of two sources and two mixtures, in which both ions have the same valences, in which case the model can be written as:

y = e + d • log (As) (3)
where A is a 2 × 2 matrix with a 11 = a 22 = 1, d is the vector containing the estimates of the slope for each ISE, and the operator • denotes the Hadamard (entrywise) matrix product 1 . While the slopes represented by d can be 1 The logarithm operation is also performed entrywise.

theoretically calculated directly by ( 2), the empirical nature of the NE equation cannot always accurately model the data, and better results can be obtained by giving the model an additional degree of freedom, allowing it to vary the slope.

In our previous work, we had estimated the d slope values using theoretical values and defined the modified mixture set x, with

x i = 10 y i -e i d i . (4)
If the slope values were exact (and the mixture were perfectly modeled by [START_REF] Achard | Criteria based on mutual information minimization for blind source separation in post nonlinear mixtures[END_REF], which is not the case), it can be shown that x = As, which is a linear mixture of the two sources which can then be solved by classical linear BSS algorithms. This results in the classical PNL approach in which the difficult part is to estimate the parameters that correctly cancel the nonlinearity -in this case the slopes d. However, in our previous work we had simply used the theoretical slope value, and together with the fact that (2) only approximately models the mixture, we can only claim that the modified set x is an approximately linear mixture of the sources. The idea, then, was to represent x as a second degree polynomial, in which the quadratic terms can correct small nonlinearities in the modified mixture. The mixture was represented then as

x = As + Bs •2 + cs 1 s 2 (5)
where a 11 = a 22 = 1 and s

•n = [s n 1 , s n 2 ]
T represents the entrywise power operator. This mixing model encompasses the class known as linear-quadratic model, for which one can find several methods to deal with [START_REF] Deville | Recurrent networks for separating extractable-target nonlinear mixtures. part i: Non-blind configurations[END_REF]. Since the mixing model is not as simple as the linear case, inverting it is no longer trivial, and a recurrent network seems indicated to solve it even in the non-blind case [START_REF] Deville | Recurrent networks for separating extractable-target nonlinear mixtures. part i: Non-blind configurations[END_REF].

This algorithm, however, only works for the simplest case in which we had two mixtures of two sources with the same valence.

B. General case

A more general case could contain any number of sources and mixtures, and the ions could have different valences. The valences of the ions are known in advance. We shall still restrain ourselves to the case in which the number of sources and mixtures are equal -that is, the determined case.

It is important to start by noticing that when the valences are different, the mixing model is no longer PNL, even if the NE equation ( 2) were able to perfectly model the mixture. This is because, in (2), we are actually mixing s zi/zj j , with z i and z j possibly different. A useful change of variable can be made by defining the pseudo-source:

p i (n) = (s i (n)) 1/zi (6)
Then (2) becomes

y i (n) = e i + RT z i F ln   N j=1 a ij p zi j (n)   . (7)
The associated modified mixture, as defined in (4), becomes then polynomial in the pseudo-sources p with degree k = max(z)/gcd(z), where gcd(z) is the greatest common denominator of all the valences of z i . Since we want to add an additional degree of freedom in order to capture the inaccuracies of the NE-model, we can therefore model the modified mixture set as a polynomial of degree k + 1.

Since the modified model is nonlinear, even in the blind case the solution is not trivial, even more so in the general case in which we can have multiple sources and different valences. Therefore, an analytical solution is generally impossible to obtain, we resort to a recurring network to achieve the solution numerically.

C. Recurring Network

Using the definition of our pseudo-sources vector p, we can rewrite the modified mixture vector x, with each component given by ( 4) as

x = F[p] (8)
where F[.] is a polynomial of the appropriate degree, including all possible cross-multiplied terms of the pseudo-sources p i . For example, for the case of 3 sources with valences +1, +2 and +3, our polynomial would have degree 4, and could be written as:

x = 4 i=1 A i p •i + Bp cross (9)
where A i are square matrices, p cross contains all the cross multiplied terms of the pseudo-sources and B is a coefficient matrix of appropriate dimension. While it can be said that the modified mixture set given by (4) should not contain cross multiplied terms, this would only be true if both the estimates of the Nernst slopes and the NE equation ( 2) were perfectly correct. Since they are in fact only empirical approximations, the additional degrees of freedom provided by the cross multiplied terms, as well as the higher order of the polynomial, can be used to correct inaccuracies in the model.

If the coefficients are known (i.e., the nonblind case), we have to solve the nonlinear equation system:

G[p] = F[p] -x = 0. ( 10
)
Applying the Newton-Rhapson root-finding algorithm to [START_REF] Duarte | A sparsity-based method for blind compensation of a memoryless nonlinear distortion: Application to ion-selective electrodes[END_REF], we obtain the following recurrent network:

q t+1 = q t -J -1 G (q t ) G[q t ] (11)
where J G is the Jacobian of G. Since G is a polynomial, this is easy to calculate. As usual, to lower the computational cost of the operation, it is better to solve the linear equation system J G (q t )∆q t = -G[q t ] instead of inverting the matrix. When the network converges, q would eventually converge to the pseudo-sources p.

However, in the problem we are trying to solve, the mixing coefficients are unknown, since the problem is blind. Thus, before being able to use this recurrent network we need a way to estimate the polynomial coefficients first.

D. Estimation of the Polynomial Coefficients

In order to estimate the coefficients of the mixture, we assume that the sources are independent, which automatically implies that our pseudo-sources are also independent. It is important to notice that source independence is a common hypothesis in many BSS problems, and is a good assumption in our application as well.

In order to maximize the independence of the estimated sources, in the first algorithm, for 2 sources, we used the mutual information as a cost function for estimating the parameters, with good results. While different generalizations of the concept of mutual information for more than two sources exists, we can use the total correlation for our problem [START_REF] Watanabe | Information theoretical analysis of multivariate correlation[END_REF]. The total correlation can be defined as the Kullback-Leibner divergence between the joint distribution and the independent distribution, which can be simplified to:

C(X 1 , X 2 , ..., X n ) = n i=1 H(X i ) -H(X 1 , X 2 , ..., X n ) (12) where C(X 1 , X 2 , ..., X n ) is the total correlation, H(X i ) = -p i (x i ) log (p i (x i )) is the information entropy of X i and H(X 1 , X 2 , ..., X n ) = - x1 ... xn p(x 1 , ..., x n ) log (p(x 1 , ..., x n )) (13)
is the joint entropy of the variables. It is important to notice that since the total correlation is a Kullback-Leibner divergence, we have C ≥ 0 with equality if, and only if, the random variables are statistically independent. Therefore, when trying to obtain polynomial coefficients that produces independent source estimates, the total correlation can be used as a cost function, and we can use the gradient method to minimize it. A detailed derivation of the differential of the mutual information (or, for the case of more than two sources, the total correlation) can be found in [START_REF] Babaie-Zadeh | Differential of the mutual information[END_REF], and can be calculated approximated by

∂C ∂w = E ∂q ∂w T β q (q) ( 14
)
where w is the vector containing all the polynomial coefficients, and

β qi (q) = - ∂ log p(q) ∂q i -- d log p(q i ) dq i ( 15
)
is what can essentially be interpreted as "gradient" of C with respect to the source estimates q [START_REF] Babaie-Zadeh | Differential of the mutual information[END_REF]. By defining the score function as the opposite of the logderivative of a density function, i.e.,

φ(x) = - d dx ln p x (x) (16)
it can be seen that β qi (q) is the difference between the ith component of the joint score function of q and the marginal score function of q i . It is interesting to note that the score functions, and, by extension, β q (q) as defined in [START_REF] Pham | Fast algorithms for mutual information based independent component analysis[END_REF], can be accurately computed by efficient methods as described, for instance, in [START_REF] Pham | Fast algorithms for mutual information based independent component analysis[END_REF], which leads to an efficient computation of our cost function.

As a result, the iteration to find the parameters can be written as:

w k+1 = w k -µ ∂C ∂w k (17)
where µ denotes the learning rate 2 . Note that the expression for ∂q ∂w would depend on the mixing model. For our problem, we can start from the modified mixing model ( 8) and apply implicit derivatives, replacing p (our desired result) by q (our estimates). For instance, for the example with 3 sources with valences +1, +2 and +3 shown in ( 9), we would have:

0 = ∂ ∂w 4 i=1 A i q •i + Bq cross = ∂ ∂w 4 i=1 A i q •i + Bq cross A i ,B constant + ∂ ∂w 4 i=1 A i q •i + Bq cross q,qcross constant (18)
where the product rule is used. The first expression, with the constant matrices, is J G ∂q ∂w . For the second expression, with the constant vectors, notice that our coefficient vector w contains all the variables from matrices A 1 , ..., A 4 , B. Let us define w so that the elements are ordered in such a way that all variables from A 1 appear first, followed by the elements of A 2 , and so on until B, and on each matrix the elements are extracted row-by-row.

The second expression, then, would be a matrix K = k ij with k ij = 0 if, and only if, the variable represented by column j appears in the equation for mixture i, and in that case, k ij is the coefficient of that variable. This can be concisely written with the use of the Kronecker product. The Kronecker product of two matrices of arbitrary sizes, denoted by A ⊗ B, is the block matrix defined by:

A ⊗ B =    a 11 B • • • a 1n B . . . . . . . . . a m1 B • • • a mn B    . (19)
The expression with constant vectors from (18) can, therefore, be written as: 2 Smaller learning rates improve the algorithm's robustness at the expense of longer convergence times. In our simulations we used µ ∈ [0.001; 0.01].

K = I 3 ⊗ q T • • • I 3 ⊗ q •4 T I 3 ⊗ q T cross = -J G ∂q ∂w (20)
where I 3 is the identity matrix of order 3 (the number of sources in our example) and, with the Kronecker product, the expression of K is written as block matrix. Using (20), we can finally isolate ∂q ∂w , as necessary for the iterations of the gradient method proposed.

E. Algorithm

The algorithm can then be summarized as follows: 1) From the mixtures y, obtain the modified mixture set x as defined by (4). 2) Start with random parameters w.

3) Estimate the pseudo sources p using the recurrent network [START_REF] Deville | Recurrent networks for separating extractable-target nonlinear mixtures. part i: Non-blind configurations[END_REF] with the weights w. 4) Update w using the gradient of the total correlation of the estimates, calculated via equations ( 15) to (20) 5) If stopping criterion is not met, return to step 3. 6) After convergence, recover the sources s from the pseudo-sources p, using equation [START_REF] Duarte | Design of smart ion-selective electrode arrays based on source separation through nonlinear independent component analysis[END_REF].

It is important to notice that while this paper shows the example of three sources with valences +1, +2 and +3, the same algorithm works with any number of sources and any valences. The number of sources and valences, however, must be known in advance since they are required to calculating the Jacobian J G [.], which is used both on steps 3 and 4.

IV. EXPERIMENTAL RESULTS

In our simulation, we used 3 sources with valences +1, +2 and +3 -a good representation of the more general case in which we have more than 2 sources and different valences, and therefore could not have been solved by the previous algorithm. Our simulation had a sample size of N = 400. In Fig. 1, we can see the results of a simulation.

We can see in the first row of the figure three independent and uniformly distributed sources in the interval [0, 1], but with the samples of the first one being sorted in ascending in order to be easier to visualize, but without affecting the independence of the sources. The sources are then mixed according to Eq. ( 2), as shown in the second row.

The estimates obtained by the algorithm can be seen in the third line, and at the end we can see a plot of the estimates and the sources. Ideally, it should be noted that we want to obtain a straight line, of any inclination, since it only indicates a different scalar factor. We can see that for the three sources, we obtain what can be roughly approximated by a straight line, indicating that the result is good. In order to measure how well the estimates represent the sources, we can use the signal to interference ratio, which can be defined, in decibels, as: where the s is the source vector, y is the estimate vector, and the higher the SIR value the more accurate the estimate is.

SIR = 10 log 10 E{s 2 } E{(s -y) 2 } (21)
For the simulation here, the SIR values obtained were 20.0 dB, 26.6 dB and 21.4 dB, which are typically considered as good values for the application analyzed. It should also be noticed that because of the random nature of the algorithm, both in the initial point for the recurrent network, and the initial set of parameters, the algorithm does not always converge. When the algorithm fails to converge, we can simply rerun it with a different set of initial conditions until an estimate is found.

V. CONCLUSION

In this paper, we have generalized an algorithm proposed in [START_REF] Ando | A blind source separation method for chemical sensor arrays based on a second order mixing model[END_REF] which was able to show promising results when applied to the measurement of ionic concentrations using chemical sensors, but was limited only to two sources that had to be ions of the same valence. Using similar ideas, we were able to generalized the algorithm to any number of sources with possibly different valences.

By allowing the valences to be different, the resulting mixture can no longer be classified as post non-linear (PNL) and, as a result, the method previously found in the literature [START_REF] Duarte | Application of blind source separation methods to ionselective electrode arrays in flow-injection analysis[END_REF] to compare with our proposed method does not work. Nonetheless, while not as good as the simpler case of 2 sources with equal valence, our simulations indicate that the separation can be performed with accuracy, obtaining SIR between 20 and 27dB, which are typically good for this sort of application.

In the future, we would like to further analyze the theoretical stability of the problem in order to improve the selection of initial conditions for the algorithm. For the current version of the algorithm, the conditions are randomly selected and can sometimes lead to the algorithm not converging, in which case we restart the method until it succeeds. While this has not proven an issue for our current application (possibly due to the fact that the mixing coefficients are typically small), a more theoretical analysis would be important to avoid problems with more general mixing models, and possibly improve the accuracy and convergence speed of our proposed method.

Fig. 1 .

 1 Fig. 1. Plots of sources and estimates


Bertrand Rivet

Jeremy E Cohen

Modeling time warping in tensor decomposition

Taking into account subject variability in data mining is one of the great challenges of modern biomedical engineering. In EEG recordings, the assumption that time sources are exactly shared by multiple subjects, multiple recordings of the same subject, or even multiples instances of the sources in one recording is especially wrong. In this paper, we propose to deal with shared underlying sources expressed through time warping in multiple EEG recordings, in the context of ocular artifact removal. Diffeomorphisms are used to model the time warping operators. We derive an algorithm that extracts all sources and diffeomorphism in the model and show successful simulations, giving a proof of concept that subject variability can be tackled with tensor modeling.

I. INTRODUCTION

Electroencephalographic (EEG) recordings of a single subject's neural activity are data sets with two dimensions, namely time and sensors. A wide variety of multilinear methods have been developed in the recent years that try to extract meaningful sources from these data that describe patterns of neurons or systematic responses to certain stimuli. In particular, to perform source separation in a blind context where no a priori knowledge on the cortex nor on the sources is available, tensor decomposition techniques applied on tensorization of two way arrays or multiple recordings of the same subject have proven a powerful tool in EEG [START_REF] Cichocki | Tensor decompositions for signal processing applications: From two-way to multiway component analysis[END_REF]. Indeed the approximate multilinearity assumption made in tensor decompositions seems to fit with underlying complex physical models of brain electrical activity. However a major obstacle to source separation in EEG is the poor signal to noise ratio due to numerous artifacts in the signal. Notably, eye movement is known to generate a signal more powerful than the signal emitted by sources of interest. In this work, source separation is used to estimate these artifacts in the EEG recordings.

For signal generated by eye-movements, it is assumed that the signal x ik (t) recorded by the ith sensor at the kth trial is the summation of R source signals b r (t) up to scaling factors a ir that are independent of the trial. However, this is not necessarily the case since temporal changes can occur due to time shift inherent to differences in speed and amplitude of saccades. The same issue of variability between several measurements of the same phenomenon can be encountered when quantifying the concentration of a chemical substance by measuring a chromatographic signal. In this case, temporal responses of certain chemicals may have various shapes depending on temperature and injection of the mixture which is done by hand. Another example is audio data, where the Doppler effect can cause time-warping.

In this paper, we tackle the problem of tensor factorization when the underlying profiles change shape from trial to trial. In this case, the canonical polyadic (CP) decomposition model is no longer valid and leads to poor estimation results [START_REF] Bro | PARAFAC2-Part II. Modeling chromatographic data with retention time shifts[END_REF]. We introduce a registered CP (rCP) decomposition model which assumes that the second mode is not only shifted from trial to trial as in [START_REF] Harshman | Shifted factor analysis-Part I: Models and properties[END_REF] but can be any time warping modeled by an orientation-preserving diffeomorphism.

The paper is structured as follows. In Section II, the used notations are introduces. Section III briefly recalls the CP decomposition and related methods to deal with profiles change in one mode. The proposed registered CP decomposition is derived in Section IV, next we validate the usefulness of our rCP decomposition (Section V). Finally, Section VI concludes this paper.

II. NOTATIONS

In the following, X 2 R M ⇥L⇥K denotes a third order tensor (here a three way array) whose (m, l, k)th entry is x mlk . The kth matrix slice of X is denoted by X k = X (:, :, k) (i.e. it contains all x mlk by keeping k constant). The column vector y r is the rth column of the matrix Y and y r is the rth element of y. The transposition operator is denoted by • 0 , the pseudoinverse matrix is denoted by • † and ⌦ and ⇥ are the outer and Kronecker products, respectively. k • k 2 and k • k F are the `2 and Frobenius norms, respectively. Diag(•) is the diagonal operator that only keeps the diagonal element of the matrix input. Finally, is the composition operator: x y(t) = x y(t) and ḟ is the derivative of f . III. RELATED TENSOR FACTORIZATIONS Let X be the three way array of the recorded data such that X k is the kth trial or sample. More precisely, on EEG data, the raw recordings are split according to the saccades so that only one saccade is present in each trial, providing a data matrix X k which are then stacked to form X , where M is the number of sensors, L the number of time samples and K the number of trials. Assuming an additive gaussian i.i.d. noise, the recordings are expressed as

X = T + E, (1)
where T and E are the signal and noise arrays, respectively. The signal array T can be factorized to reveal the underlying latent factors. It is assumed that the shapes of sources in the second mode can change from sample to sample. The goal of this paper is to incorporate a priori knowledge on the variations among the shapes into tensor decomposition models, and to pave the way for further research in modeling subject variability in tensor decomposition.

Differences between shapes of second mode sources are modeled as a shape warping from a latent set of 'original' sources b r (t):

bkr (t) = b r rk (t) , (2)
where rk (t) is an orientation-preserving diffeomorphism (see Fig. 1).

In addition to X , a set of P others signals Y (p) 2 R M ⇥K are simultaneously recorded that allows additional a priori information on the latent sources of the second mode to be obtained. In the case of EEG ocular artifacts removal, GAZE signals (i.e. eye motions) are simultaneously recorded. The fusion of these data sets will be the topic of further research, since here the second data set serves only to be used as an a priori on the latent sources b r or diffeomorphisms kr .

A. Canonical polyadic decomposition

The CP decomposition represents an array T by a sum of rank-one tensors as

T = R X r=1 r a r ⌦ b r ⌦ c r . (3)
Equivalently, the kth matrix slice is expressed as

T k = AD k B 0 , (4)
where D k is a diagonal matrix whose entries are the kth row of C⇤, with ⇤ the diagonal matrix whose entries are r . The loading matrices A, B and C and the scaling factors are estimated as

ˆ , Â, B, Ĉ = arg min ,A,B,C 1 2 X T 2 F (5) subject to 8r, ka r k 2 = kb r k 2 = kc r k 2 = 1.
Consequently, the CP decomposition model assumes that each trial can be factorized as a matrix product of spatial map A and temporal shapes B up to scaling factors. An underlying hypothesis of CP decomposition model is therefore that the spatial maps and temporal shapes are the same for each trial: i.e. kr = 0 , with 0 the delta Dirac distribution at 0. However, this latter assumption is not always satisfied if the trials are misaligned and/or if the temporal shapes presents some time warping between them, like temporal shifts or time contractions/dilatations. Especially for outliers among time sources generated by eye movement, variability in eye saccades leads to variability in the temporal responses that are not taken into account by the CP decomposition. The issue is similar to the time retention shift issue encountered with chromatographic data [START_REF] Bro | PARAFAC2-Part II. Modeling chromatographic data with retention time shifts[END_REF].

B. PARAFAC2, SCREAM and Shift-PARAFAC models

To overcome these limitations of the CP decomposition, the PARAFAC2 [START_REF] Harshman | Parafac2: Mathematical and technical notes[END_REF] model was introduced. It is usually defined by

8k, T k = AD k B 0 k (6)
where D k is a diagonal matrix with entries the kth row of C⇤. This model differs from (4) only in the fact that the second loading matrix can be different for each trial. To ensure identifiability, the following constraints are added: B 0 k B k = I, for all k. This model thus assumes that the second loading matrices share the same correlation between trials:

Â, { Bk }, Ĉ = arg min 1 2 X T 2 F , (7)
s.t. 8k, B 0 k B k = I, with T defined by [START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF]. This model better suits the assumed data structure in (2) but it does not take into account the distortion of the time axis. In other words, PARAFAC2 assumes a constant autocorrelation of the sources in B k which is a wrong assumption if the variability among sources is due to time warping.

Another closely related model is the SCREAM [START_REF] Marini | SCREAM: A novel method for multi-way regression problems with shifts and shape changes in one mode[END_REF] model which was introduced to model some data with shifts and shape changes in one mode, but provided auxiliary matrix data for which a linear regression model depending on the latent sources is valid. It is computed as a PARAFAC2 model coupled to a linear regression model with mixing matrix factor C.

Finally when sources in B k are identical but shifted in time with delays ⌧ kr , i.e. kr = ⌧kr , parametrization of the delay in the indexes of columns B r leads to a model called Shift-PARAFAC [START_REF] Harshman | Shifted factor analysis-Part I: Models and properties[END_REF]. For EEG and eye artifact extraction, again this model cannot be applied since a simple time shift between sources is not a satisfactory model.

, T k = AD k B0 k , (8)
where

Bk = [ k1 b 1 , • • • , k,R b R ].
The matrices kr are the warping matrices (i.e. matrix representations of diffeomorphism operators).

A. Registered CP algorithm

The registered CP (rCP) decomposition aims at optimizing the following criterion:

Â, B, Ĉ, ˆ rk r,k = arg min A,B,C,{ rk }r,k 1 2 X T 2 F , (9)
with T defined by [START_REF] Rao | Information and the accuracy attainable in the estimation of statistical parameters[END_REF]. To optimize this criterion, an alternative least square (ALS) algorithm is used [START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF]: it aims at estimating sequentially each parameter with all the other parameters fixed. Consequently, we propose the following step:

1) Update A given B, Ĉ, ˆ kr  K X k=1 X k Bk Dk ! K X k=1 Dk B0 k Bk Dk ! 1 . ( 10
) 2) Update C given Â, B, ˆ kr Dk Diag ⇣ Â † X k B † 0 k ⌘ . (11)
3) i) Estimate Bk , for all k, given Â, Ĉ, ˆ kr

Bk X 0 k  † D † k . (12)
ii) Factorize each matrix

h b1r , • • • , bKr i as r (I K ⇥ b r ), where r = [ 1r , • • • , Kr ] is
8k, Bk = ⇥ ˆ k1 b1 , • • • , ˆ kR bR ⇤ . (14)

B. Factorization of warped loadings

The underlying structure of the warped loadings Bk is recovered by estimating both the generating function b r and the warping diffeomorphisms kr [START_REF] Tucker | Generative models for functional data using phase and amplitude separation[END_REF]. Let F be the set of absolutely continuous functions defined on [0, 1] and G be the set of orientation-preserving diffeomorphisms of the unit interval [0, 1] to [0, 1].

The aim of step 3) defined by ( 13) is to find a function b r 2 F so that it is the best approximation of all functions bkr up to warping functions kr 2 G. To do this, let us use an elastic distance between functions, i.e. a distance in the quotient space F/G which is the set of orbits of the type

[f ] = {f | 8
2 G}, for any f 2 F. The used elastic distance between two functions f 1 and f 2 of F is defined by the Fisher-Rao Riemannian metric [START_REF] Rao | Information and the accuracy attainable in the estimation of statistical parameters[END_REF], [START_REF] Rao | Breakthroughs in Statistics: Foundations and Basic Theory[END_REF]:

d(f 1 , f 2 ) = inf 2G d F R (f 1 , f 2
).

However, as pointed out by [START_REF] Tucker | Generative models for functional data using phase and amplitude separation[END_REF], it has been shown [START_REF] Bhattacharyya | On a measure of divergence between two multinomial populations[END_REF] that this distance can be efficiently computed by transforming the functional space F into L 2 , the set of square integrable functions, using the square-root slope function (SRSF): q = sign ḟ q ḟ , with sign(•) the sign function. Indeed, in L 2 , the Fisher-Rao distance becomes the Euclidian distance:

d F R (f 1 , f 2 ) = kq 1 q 2 k
. Consequently, the latent loading b r is estimated as the Karsher mean of { b1r , • • • , bKr } defines from its SRSF as qr = arg min

q2L 2 K X k=1 ✓ inf k 2G q (q k k ) p ˙ k 2 ◆ , ( 15
)
where q k is the SRSF associated with bkr . The set of warping diffeomorphisms kr are then estimated as

ˆ kr = arg inf qkr (q r ) p ˙ 2 , (16)
by dynamical programming [START_REF] Bertsekas | Dynamic Programming and Optimal Control[END_REF]. Finally, b r is expressed as

br = K X k=1 c 2 kr (a 0 r a r ) ˆ 0 kr ˆ kr ! 1 ⇥ K X k=1 c kr ˆ 0 kr ✓ X 0 k X l6 =r c kl ˆ kl b l a 0 l ◆ a r ! . ( 17
)
It is worth noting that some trivial underdetermination are inherent to the assumption (2). Indeed, one can blindly only recover the latent sources b r and the diffeomorphisms kr up to diffeomorphisms

˜ kr since bkr = b r kr = b r ˜ 1 kr ˜ kr kr .
C. Using priors on b r or kr

The above scheme assumes that no prior information is known about neither the latent factor shapes b r nor the warping diffeomorphisms kr . However, in some experiments, such prior information can be available. This is for instance the case if GAZE&EEG data are recorded synchronously. Indeed, the GAZE recordings can be used either to initialize the iterative ALS algorithm or as prior on either b r or kr .

If all warping diffeomorphisms kr are known, then the third step [START_REF] Kristensen | Comparison between Adjar and xDAWN algorithms to estimate eye-fixation related potentials distorded by overlapping[END_REF] resumes to (17) by substituting ˆ kr by its actual value kr . Similarly, if latent factors b r are known, then the third step [START_REF] Kristensen | Comparison between Adjar and xDAWN algorithms to estimate eye-fixation related potentials distorded by overlapping[END_REF] resumes to (16) by substituting qr with q r , the SRSF of b r .

V. RESULTS

Both numerical simulations and an illustration on actual GAZE&EEG data are presented.

A. Numerical simulations

In each simulation, the number of components R is fixed to 3. The length of the first, second and third loading factors a r , b r and c r are 10, 250 and 50, respectively. The first and third loadings factors A and C are generated randomly from a standard normal distribution. Each component b r of the second loading factor B is generated from a Gaussian process [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] with a zero mean function and a covariance function k r (t, t 0 ) = exp (t t 0 ) 2 /2l 2 r , with l r uniformly distributed in [.02; .22], so that it is smooth enough with still some variations. The warping diffeomorphisms are generated as

kr (t) = Z t 0 ↵ exp w kr (u) du,
where ↵ is a normalisation factor so that kr (1) = 1. Functions w kr are randomly drawn as Gaussian processes with a zero mean value function and a covariance function

k (w) kr (t, t 0 ) = 2 exp (t t 0 ) 2 /2
2 kr , with = .8 and kr uniformly distributed in [.1; .4]. Typical functions b r and warping diffeomorphisms kr are plotted in Fig. 1. Entries of tensor E are drawn from a zero mean normal distribution with a standard deviation equal to .5. For each numerical experiment, 50 configurations are randomly drawn by first regenerating random parameters l r , kr and all warping diffeomorphisms kr and then each factors A, B and C and additive noise E.

The quality of the tensor estimation is quantified by two performance indexes. The overall signal-to-noise ratio improvement (⇧ O ) is defined as

⇧ O T = 10 log T 2 F T T 2 F ! 10 log T 2 F T Tinit 2 F ! ,
where T , T and Tinit are the true, estimated and initial tensors, respectively. The higher, the better the results. The factor performance index (⇧ F ), which measures the quality estimation of each factor, is defined as

⇧ F = 8 > < > : 1 d0 r d r k dr k 2 kd r k 2 , if d r is a non warped factor, inf 2G d F R ( dr , d r ), if d r is a warped factor.
Based on the trivial underdetermination of the warped factor, its performance is measured as the elastic distance in the quotient space F/G, while the other ones are based on the correlation coefficient between the true loading and its estimate. For the factor performance index, the closer to zero, the better the results.

The numerical results are shown on Fig. 2. As one can see, the proposed rCP allows to estimate the underlying factors (Fig. 2(b), 2(c), 2(d)). Indeed, the performance indexes ⇧ F after the estimation are lower than the ⇧ F at the initialization. Both methods rCP and rCP-B provide quite similar performance but lower than rCP-G which assumes known the diffeomorphisms. Quite surprisingly, rCP-B leads to a little bite lower performance than rCP: such results will be more deeply investigated in a future work. The proposed rCP outperforms both the classical CP and PARAFAC2: after convergence, the estimated tensors are worse than the initialized one (the ⇧ O is negative). Indeed, for these two methods, this result is explained by a worse estimation on the second loading factor (Fig. 2(c)). For the CP, one can explain this by the missmodeling which does not assume some shapes modifications in any loadings. This is expected on the warped loading (Fig. 2(c)). Consequently, the first one, which is theoretically not affected by a miss-modeling, and the third one (Fig. 2(d)) which models the amplitude of factors are also less accurately estimated (Fig. 2(b)) than using rCP. The more flexible PARAFAC2 model also fails to estimate the underlying latent factors B since the orthogonality constraint is not valid but its flexibility leads to a better estimation of the first and third loadings than CP.

B. Illustration on GAZE&EEG data

The Gaze&EEG data comes from experiments in visual search where participants had to search a target from a set of distractors [START_REF] Kristensen | Comparison between Adjar and xDAWN algorithms to estimate eye-fixation related potentials distorded by overlapping[END_REF]. GAZE and EEG data were recorded synchronously. From these data, two tensors (space ⇥ time ⇥ trials) are built synchronized on the beginning of the saccades: each tensor is the concatenation in the third dimension of the matrices in which each column is the vector of channel signals at a given time. The number of components is fixed to R = 2, EEG tensor is of size 63 channels, 250 time samples and 47 trials. The initializations of kr have been done from the GAZE signals. As one can see (Fig. 3), after ocular removal by rCP, the remaining EEG is much more stationary than using CP. rCP removes almost all the artefacts while CP does not. VI. CONCLUSION In this article, the factorization of a three-way array has been proposed for the case that one mode presents shape variations. The proposed registered-CP decomposition processes the mode with shape variations as orientation-preserving diffeomorphisms. It is thus a generalization of the Shift-PARAFAC method since it considers a larger class of shape variations and it is related to the PARAFAC2 model, but different constraints are imposed. An ALS scheme has been proposed to estimate the loadings and the latent shapes with the associated diffeomorphisms. The numerical simulations show the potential of the rCP which outperforms both the CP and the PARAFAC2 decompositions.

The current algorithm suffers from its computational cost due to the implementation using dynamic programming. Thus some efficient algorithms must be proposed to increase the efficient use of the rCP algorithm. In perspective, some applications like chromatographic data require additional assumptions on some loadings such as non-negativity, those cases will be considered in the future. Finally, a theoretical study must be made to provide some sufficient/necessary conditions of identifiability of the proposed rCP model.

Fig. 1 .

 1 Fig. 1. Typical functions br, warping diffeomorphisms kr and warped factors bkr . Fig. 1(a) shows 3 latent factors. Fig. 1(b) shows 4 diffeomorphisms kr and Fig. 1(c) shows b 1 k1 , i.e. the application of the 4 diffeomorphisms on the latent factor b 1 in blue in Fig. 1(a).

1 kr = b lr 1 lr

 11 IV. PROPOSED REGISTERED CP DECOMPOSITIONIn this paper, it is assumed that the second loading factors b r come from the sampling of functions b r and are equal from trial to trial up to orientation-preserving diffeomorphisms •r : i.e. 8r and 8(l, k), b kr . Consequently, the kth trial T k is expressed as 8k

d kr b r , bkr 2 ,

 2 the concatenation of warping matrices for the rth component of the second loading factor, by br , ˆ kr k = arg min br, kr k K X k=1 (13) where d(•) is a distance defined in the next paragraph. iii) Then, Bk is updated as

Fig. 2 .Fig. 3 .

 23 Fig.2. Performance of the proposed rCP. Comparison of performance indexes for several assumptions from initialization (Init.), usual CP (CP), PARAFAC2 (PA), proposed registered CP (rCP), rCP with known latent factors B (rCP-B) and rCP with known warping functions kr (rCP-G). On each box, the red mark is the median, the edges of the box are the 25th and 75th percentiles, the whiskers extend to extreme values and the outliers are plotted individually as red crosses.

ACKNOWLEDGMENT

This work was funded by ERC-AdG-2012-320864 CHESS project and ERC-AdG-2013-320594 DECODA project.