Rémi André

Xavier Luciani

Laurent Albera

Eric Moreau

A Two-step Algorithm for Joint EigenValue Decomposition -Application to Canonical Polyadic Decomposition of Fluorescence Spectra

Keywords: Fluorescence spectroscopy, CPD, PARAFAC, Multiway analysis, Joint diagonalization, JEVD, Positive definite matrix

In this paper, we propose a new Joint EigenValue Decomposition (JEVD) algorithm. JEVD problem belongs to the family of joint diagonalization problems. Hence, JEVD algorithms aim at estimating the common basis of eigenvectors of a matrix set. This problem occurs in many signal processing applications. It has notably allowed to develop efficient algorithms for the Canonical Polyadic Decomposition (CPD) of multiway arrays. The proposed JEVD algorithm is based on an original two-step approach. The first step consists in transforming the considered matrix set into a set of positive definite matrices. In this purpose, we introduce an ad hoc joint symmetrization algorithm. This first step allows us to transform the JEVD problem into a simpler orthogonal joint diagonalization problem. The second step is then performed using an efficient orthogonal joint diagonalization algorithm of the literature. Eventually, the performance of the proposed approach is deeply investigated in the CPD context of multidimensional fluorescence data. More particularly, we consider difficult scenarios such as the cases of an overestimated rank and highly correlated factors.

Introduction

In many data sciences or signal processing applications, the collected data has a multidimensional structure and can thus be stored in multiway arrays (tensors). In this context, multilinear algebra provides efficient tools to analyze such data sets. Indeed, tensor decompositions offer a large variety of multidimensional models with some interesting properties. Among all tensor decomposition models, Canonical Polyadic Decomposition (CPD) also known as PARAllel FACtor decomposition (PARAFAC) [START_REF] Hitchcock | Multiple invariants and generalized rank of a p-way matrix or tensor[END_REF][START_REF] Harshman | Foundation of PARAFAC procedure: Models and conditions for an 'explanatory' multi-mode factor analysis[END_REF] decomposes a multiway array as a linear combination of a finite number of rank-one tensors. Moreover, this decomposition is essentially unique under some mild conditions [START_REF]Three-way arrays: Rank and uniqueness of trilinear decompositions[END_REF][START_REF] Sidiropoulos | On the uniqueness of multilinear decomposition of n-way arrays[END_REF]. For these reasons, CPD is largely used to estimate multilinear mixtures of source signals and plays a growing role in signal and image processing such as in brain source localization [START_REF] Hunyadi | Block term decomposition for modelling epileptic seizures[END_REF][START_REF] Becker | Brain source imaging: from sparse to tensor models[END_REF], face recognition [START_REF] Hao | Facial recognition using tensor-tensor decompositions[END_REF] and blind source separation [START_REF] Sidiropoulos | Blind PARAFAC receivers for DS-CDMA systems[END_REF][START_REF] Comon | Handbook of blind source separation[END_REF]. In chemometrics, CPD directly models the data measured from fluorescence spectroscopy [START_REF] Bro | Parafac, tutorial and applications[END_REF]11].

CPD decomposes multiway arrays as a combination of factors matrices. For instance, the rank-N CPD of a three way tensor T ∈ R P×Q×R is given by: T p,q,r N n=1 A p,n B q,n C r,n , ∀(p, q, r)

where matrices A, B and C are the factor matrices. It is noteworthy that the smallest value of N for which equation ( 1) is exact is called the tensor rank.

Many algorithms have been proposed in order to compute the CPD of a given tensor. All of them aim at minimizing iteratively a specific cost function built (in)directly from the data. They can be classified in three main families. First, 25 there are alternating or block coordinate approaches such as Alternating Least Squares (ALS) [START_REF] Harshman | Foundation of PARAFAC procedure: Models and conditions for an 'explanatory' multi-mode factor analysis[END_REF]. These algorithms update different subsets of factors alternatively at each iteration. Secondly, we find the well-known first/second order optimization methods based on (conjugate) descent [12], Quasi-Newton 30 and Levenberg-Marquardt [13] optimization strategies. Finally, semi-algebraic approaches are also available, which generally rephrase the CPD as a Joint Diagonalization (JD) problem [14,15,[START_REF] Luciani | Canonical polyadic decomposition based on joint eigenvalue decomposition[END_REF][START_REF] Hajipour | Canonical polyadic decomposition of complex-valued multi-way arrays based on simultaneous schur decomposition[END_REF].

Joint diagonalization consists in finding a matrix, which transforms (by congruence or by similarity) a given set of square matrices into a set of diagonal matrices. In this paper, we focus on the JD problem by similarity well-known as the Joint EigenValue Decomposition (JEVD) problem. In other words, JEVD consists in jointly diagonalizing a set of non-defective matrices of size N × N sharing the same eigenvectors basis W:

M (k) = WD (k) W -1 , ∀k = 1, . . . , K. (2)
where D (k) is a diagonal matrix containing the eigenvalues of 35 M (k) . It has been shown in [START_REF] Luciani | Canonical polyadic decomposition based on joint eigenvalue decomposition[END_REF], through numerical simulations, that CPD algorithms based on JEVD have several advantages. First, they need few iterations to converge. Second, they are relatively insensitive to correlated factors and to an overestimation of the CPD rank (overfactoring problem). On the 40 other hand, JEVD-based CPD methods may require more re-strictive conditions (see [16, section 3.5] for a comparative analysis).

A naive approach to solve the JEVD problem consists in computing the eigenvalue decomposition of only one matrix, for instance M (k) or the product of several M (k) matrices. However, this approach is very sensitive to the noise and has no unique solution in the presence of a degenerate eigenvalue (i.e. multiplicity strictly greater than one). Indeed, JEVD has less restrictive uniqueness conditions [START_REF] Lathauwer | Computation of the canonical decomposition by means of a simultaneous Schur decomposition[END_REF]. Note that by uniqueness we mean essential uniqueness, that is to say the matrix W of eigenvectors can be estimated only up to a diagonal and a permutation matrix. Consequently, solving the JEVD problem by means of a single EVD is not a good thing. Generally, the JEVD problem is solved by minimizing a diagonality criterion applied simultenously to multiplicative updates of the M (k) matrices. The main differences between the existing JEVD algorithms reside in the parametrization of the updating matrix and the implemented strategy to compute these parameters.

When the eigenvector matrix W is orthogonal, the JEVD problem becomes an Orthogonal Joint Diagonalization (OJD) problem. The latter is a well-know JD problem for which many algorithms have been proposed. One of the most famous algorithm is JAD (Jacobi Angles for simultaneous Diagonalization) [START_REF] Cardoso | Jacobi angles for simultaneous diagonalization[END_REF]. In [START_REF] Flury | An algorithm for simultaneous orthogonal transformation of several positive definite symmetric matrices to nearly 700 diagonal form[END_REF] and [START_REF] Pham | Joint approximate diagonalization of positive definite Hermitian matrices[END_REF], authors focused on the approximate OJD of a set of positive definite matrices by minimizing an original cost function based on the log-likelihood approach. These algorithms are particularly efficient and have good convergence properties. Indeed, the quadratic local convergence of Pham's method is theoretically proven while simulations show that a few iterations are needed to ensure convergence [START_REF] Pham | Joint approximate diagonalization of positive definite Hermitian matrices[END_REF]. In order to take benefit of these attractive properties, we introduce a novel two-step JEVD approach transforming the original matrix set into a set of symmetric positive definite matrices. For this purpose, we introduce the problem of joint symmetrization and design an ad-hoc algorithm. We show the good behavior of the proposed JEVD method when it is combined with DIAG in order to canonically decompose data measured from fluorescence spectroscopy. In addition, we show the ability of our method to face difficult cases encountered in fluorescence measurements.

The paper is organized as follows. We first recall how the CPD problem can be rewritten into a JEVD problem. Then we describe, in section 3, the proposed approach to solve the JEVD problem. Section 4 presents a comparative computational analysis of i) the proposed JEVD method and ii) the DIAG method using this JEVD solution. The last part of the paper is dedicated to a chemometric application. Thus, in section 5, we first recall the CPD model of three-way fluorescence measurements. Then, we investigate the behavior of the proposed approach to decompose the observed data tensors in comparison with some other existing algorithms.

Notations. In the following, scalars are denoted by a lower case (x), vectors by a boldface lower case (x), matrices by a boldface upper case (X) and tensors by an upper case boldface calligraphic (X). x i is the i-th element of vector x, X i, j is the (i, j)-th component of matrix X and X i, j,k is the (i, j, k)-element of third order tensor X. Vectors x i,. and x ., j are the i th row and j th col-umn of matrix X respectively. Matrix X .,.,r is a r th frontal slice of tensor X. Operator [X] + projects the entries of X on R + by forcing all negative entries of X to 0. Diag{X} is an operator that 100 set to zero the off-diagonal entries of matrix X whereas diag{x} creates a diagonal matrix from the entries of vector x. The euclidean norm is denoted by ||.|| 2 . is the Khatri-Rao product, defined as a column-wise Kronecker product. Eventually, X (k) denotes the set of K matrices indexed by k:

X (1) , • • • , X (K) .

From CPD to JEVD

There are several ways to compute CPD by means of JEVD. The first step consists in compressing the considered tensor. Different approaches can be used to do that. Hence, in the SECSI framework [15], JEVD arises from the Higher-Order 110 Singular Value Decomposition (HOSVD) of the tensor while in the DIAG [START_REF] Luciani | Canonical polyadic decomposition based on joint eigenvalue decomposition[END_REF] and SSD-CP [START_REF] Hajipour | Canonical polyadic decomposition of complex-valued multi-way arrays based on simultaneous schur decomposition[END_REF] algorithms, it comes from the Singular Value Decomposition (SVD) of one unfolding matrix. Hereafter, we describe the latter approach as example but note that the proposed JEVD algorithm can be applied to both 115 approaches.

As mentioned in introduction, a fluorescence third order tensor is considered in this paper, for which the rank-N CPD is given by equation [START_REF] Hitchcock | Multiple invariants and generalized rank of a p-way matrix or tensor[END_REF]. By juxtaposing the R frontal slices, the tensor T can be reshaped into the unfolding matrix T of size P × QR where the second and third modes of T are merged. The CPD is then rewritten in matrix form as:

T = A C B T . (3)
Two other unfolding matrices are available. As explained in [START_REF] Luciani | Canonical polyadic decomposition based on joint eigenvalue decomposition[END_REF], it is crucial to use the squarest unfolding matrix or at least the unfolding matrix which is rank-N. Without loss of generality, we consider in this section that it is given by (3).

120

In order to transform the CPD problem into a JEVD problem, a rank-N truncated SVD of T is computed. The three matrices of this truncated SVD are denoted by U, Σ and V T , respectively. Then there exists a non-singular square matrix W of size N × N such that

A = UΣW and C B = VW -T . (4)
Furthermore C B T can be seen as a horizontal block matrix:

C B T = Φ (1) B T , • • • , Φ (R) B T , (5)
where Φ (1) , • • • , Φ (R) are the R diagonal matrices built from the R rows of matrix C. Then, (4) and ( 5) yield

V T = Γ (1)T , • • • , Γ (R)T , (6)
where

Γ (r 1 ) = BΦ (r 1 ) W T for r 1 = 1, • • • , R. Assuming that matrix B is full column rank, we can define for all couples k = (r 1 , r 2 ) with r 1 = 1, • • • , R and r 2 = 1, • • • , R, r 1 r 2 ,
a set of matrices jointly diagonalizable by similarity transformation :

125 M (k) = Γ (r 1 ) Γ (r 2 ) T , (7)
= WD (k) W -1 , (8)
where is the Moore-Penrose inverse and D (k) = Φ (r 2 ) Φ (r 1 ) so that D (k) is a set of diagonal matrices. Hence, performing JEVD of the K = R(R-1) matrices {M (k) } allows to estimate W. This quadratic relationship between the number of matrices in the JEVD step (K) and the dimension of the last mode (R) has two consequences when R is large. On one hand, it achieves a better estimation of W. On the other hand, the numerical complexity quadratically increases with R (see section 4). Thus one has to find a good trade-off between speed and estimation accuracy. As a rule of thumb, if K is very large with respect to the matrix size N, one can reasonably leave out a significant fraction of the matrix set. Here, we apply an arbitrary matrix decimation as a naive approach but better strategies could be used. For instance, we could compute and use a basis of the vector subset spanned by the K matrices to be diagonalized.

At this stage, there are several ways to deduce the factor matrices from matrix W and/or the set D (k) . DIAG resorts to the eigenvector matrix W and equations (4) to immediately estimate matrices A and C B. Then B and C are deduced as explained in [START_REF] Luciani | Canonical polyadic decomposition based on joint eigenvalue decomposition[END_REF]. Conversely to the DIAG procedure, SSD-CP resorts to the eigenvalues contained in matrices {D (k) }. Indeed diagonal entries of these matrices are linearly linked to the entries of C so that a least square estimate of C can be obtained from an overdetermined linear system. B A is then easily deduced, hence B and A. Of course these different schemes could be combined in order to define hybrid versions of DIAG and SSD-CP but this is outside the scope of this paper.

JEVD based on joint symmetrization

A two-step framework

We consider a set of K matrices of size N × N:

M (k) = WD (k) W -1 ∀k. (9)
Any non-singular matrix can be decomposed as a product of a lower triangular matrix and an orthogonal matrix. Let L and Q be the matrices of the QL decomposition of W -1 : L is a lower triangular matrix, Q is an orthogonal matrix and W -1 = QL. Thus the set of equations ( 9) can be rewritten as

M (k) = L -1 Q T D (k) QL, ∀k. (10)
Then, defining matrices S (k) as

S (k) = LM (k) L -1 , ∀k, (11)
leads to a JEVD with orthogonal transformation:

S (k) = Q T D (k) Q, ∀k. (12)
As a consequence, in the noise free case, all matrices S (k) are symmetric and equation (12) defines the OJD of S (k) . We thus propose a two step algorithm. First, we look for a lower triangular matrix L that jointly "symmetrizes" M (k) in order to estimate the set S (k) or an approximation of it in the noisy case. After this Joint Symmetrization (JS) step, JEVD problem is transformed into an OJD problem. An estimate of Q, denoted Q, is then obtained by solving this OJD problem (using, for instance, the JAD algorithm). We call this approach, JDJS for Joint Diagonalization based on Joint Symmetrization. Nevertheless, we aim at using the algorithm introduced by Pham in [START_REF] Pham | Joint approximate diagonalization of positive definite Hermitian matrices[END_REF]. Our motivation is to enjoy the nice convergence properties of this algorithm. Indeed, Pham has proven its convergence when the matrices of the set can be nearly jointly diagonalized. He has also provided a theoretical investigation showing the quadratic convergence of the algorithm close to the diagonalizing solution. Eventually, he has shown, through numerical simulations, that it requires very few iterations to reach the convergence. However, Pham's algorithm is designed for positive definite matrices. This is due to the used diagonality criterion arising from the Hadamard's inequality [START_REF] Cover | Elements of Information Theory (Wiley Series in Telecommunications and Signal Processing[END_REF]. Let S 2 (k) be a set of positive definite matrices, the Hadamard's inequality leads then to

det Diag QS 2 (k) Q T ≥ det QS 2 (k) Q T . ( 13
)
The equality being reached if and only if the matrix QS 2 (k) Q T is diagonal, a natural diagonality criterion based on the likelihood maximum method directly comes as

1 2 K k=1 log det Diag QS 2 (k) Q T -log det QS 2 (k) Q T . (14)
In our approach, we can build the S 2 (k) matrix set directly after the JS step as:

S 2 (k) = S (k) T S (k) ∀k. (15)
Hence, we obtain positive definite matrices jointly diagonalizable in a common orthogonal basis having the following expression:

S 2 (k) = Q T D (k) 2 Q, ∀k. (16)
This second version of JDJS is called JDJS 2 .

155

It is noteworthy that the matrix Q is unique (up to permutation and scaling indeterminacy) in equation ( 12) if and only if the rows of the matrix

Ω =              D (1) 11 • • • D (K) 11 . . . • • • . . . D (1) NN • • • D (K) NN              (17)
are distinct [START_REF] Lathauwer | Computation of the canonical decomposition by means of a simultaneous generalized schur decomposition[END_REF]. Hence, the transformation in equation ( 16) should not impact the identifiability of the matrix Q excepted if and only if one row of the matrix Ω is the opposite of another.

In the next subsection we focus on the JS step which is at the center of the proposed approach.

Joint symmetrization

We call joint symmetrization any process that aims to jointly transform a set of square matrices into a set of symmetric matrices and we call symmetrized matrix any output of a symmetrization process.

The proposed method to build {S (k) } from {M (k) } is a block coordinate method based on multiplicative updates. Block coordinate methods consist in estimating several unknown parameter subsets instead of estimating the whole of them at each optimization step. Here, it means that the matrix L is decomposed as the product of elementary matrices. It is noteworthy that contrarily to JD problems, the JS problem is not invariant to diagonal matrix transformations. That is why, we have to resort to two kinds of elementary matrices:

• Real elementary lower triangular matrices, denoted L (i j) (x), (i > j), are equal to the identity matrix at the exception of the L (i j) i, j term which is equal to x ∈ R.

• Real elementary diagonal matrices, denoted Λ ( j) (a), are equal to the identity matrix at the exception of the D ( j) j, j term which is equal to a ∈ R.

The following proposition gives then a simple factorization of any lower triangular matrix : Proposition 1. Let L be a lower triangular matrix of size N, then there exists N(N -1)/2 elementary lower triangular matrices L (i j) (x (i j) ) and N elementary diagonal matrices Λ ( j) (a ( j) ) such that

L = N-1 j=1 N i= j+1 L (i j) (x (i j) ) Λ ( j) (a ( j) ) Λ (N) (a (N) ). (18)
A constructive proof is trivial. Thus, such a factorization gives us a natural choice for the construction of L by multiplicative updates. Computation of each elementary matrix will be done by successive optimization of a symmetrization criterion with respect of only one parameter. We call a multiplicative update of the set M (k) by matrix X the whole following similarity transformations:

M (k) ← XM (k) X -1 ∀k. (19)
In order to compute the matrix X, we propose to minimize the following symmetrization criterion:

               C(X) = K k=1 N-1 p=1 N q=p+1 (N (k) p,q -N (k) q,p ) 2
where

N (k) = XM (k) X -1 . ( 20
)
This criterion is non-negative and is equal to zero if and only if all matrices N (k) are symmetric. According to proposition 1 and equation ( 11), the joint symmetrization of {M (k) } can be done thanks to a set of elementary lower triangular matrices denoted {L (i j) (x (i j) )} and a set of elementary diagonal matrices denoted {Λ ( j) (a ( j) )}. More precisely, at each couple (i, j), j = 1, . . . , N -1, i = j + 1, . . . , N corresponds a multiplicative update by L (i j) (x (i j) ). In addition , when i = N the update by L (N j) (x (N j) ) is followed by an update by Λ ( j) . Finally the last update is done by Λ (N) (a (N) ). We now describe how the elementary triangular and diagonal updating matrices are computed. The goal of an update is to make the whole set {M (k) } closer to a set of symmetric matrices. In this purpose, each updating matrix is chosen by minimizing the symmetrization criterion [START_REF] Flury | An algorithm for simultaneous orthogonal transformation of several positive definite symmetric matrices to nearly 700 diagonal form[END_REF]. We first consider the update by elementary lower triangular matrix L (i j) (x) i.e. X = L (i j) (x) in [START_REF] Flury | An algorithm for simultaneous orthogonal transformation of several positive definite symmetric matrices to nearly 700 diagonal form[END_REF]. Thus, we have

N (k) = L (i j) (x)M (k) L (i j) (x) -1 . (21)
We can easily check that this update only impacts the j th row and the i th column of M (k) so that

C(L (i j) (x)) = K k=1 p i p j (M (k) p,i -N (k) i,p ) 2 + (N (k) p, j -M (k) j,p ) 2 + (N (k) i, j -M (k) j,i ) 2 + β ( 22
)
where β is a constant independent of x and

               ∀ p j, N (k) i,p = -xM (k) j,p + M (k) i,p , ∀ p i, N (k) p, j = xM (k) p,i + M (k) p, j , N (k) i, j = -x 2 M (k) j,i + x(M (k) i,i -M (k) j, j ) + M (k) i, j . (23)
Thus, we have

C(L (i j) (x)) = P(x) = α 4 x 4 + α 3 x 3 + α 2 x 2 + α 1 x + β (24)
with

195 α 4 = K k=1 (M (k) j,i ) 2 ; α 3 = -2 K k=1 M (k) j,i (M (k) i,i -M (k) j, j ); ( 25
)
α 2 = K k=1 -2M (k) j,i (M (k) i, j -M (k) j,i ) + (M (k) i,i -M (k) j, j ) 2 + p i p j M (k) j,p ) 2 + (M (k) p,i ) 2 ; ( 26
)
α 1 = 2 K k=1 (M (k) i,i -M (k) j, j )(M (k) i, j -M (k) j,i ) + p i p j M (k) p,i (M (k) p, j -M (k) j,p ) -M (k) j,p (M (k) i,p -M (k) p,i ) . ( 27
)
β is a constant term thereby its expression is useless here. P(x) is a coercive 4 degree polynomial, as a consequence we choose x (i j) as the real root of the derivative of P(x) that minimize P(x).

We now explain how to compute matrix Λ ( j) (a ( j) ) for a given j. We follow the same idea as previously. The update by Λ ( j) (a) (i.e. taking X = Λ ( j) (a) in equation ( 20)) gives

N (k) = Λ ( j) (a)M (k) Λ ( j) (a) -1 . ( 28
)
We then easily show that

C(Λ ( j) (a)) = K k=1 p j (N (k) j,p -N (k) p, j ) 2 + β" (29) = K k=1 p j         M (k) j,p a -aM (k) p, j         2 + β" (30) = f (a) + β" (31)
where β" is a constant independent of a. This yields

∀a 0, f (a) = 0 ⇔ a 4 = K k=1 p j (M (k) j,p ) 2 K k=1 p j (M (k) p, j ) 2 . ( 32
)
f (a) being an even function, we can take

a ( j) =         K k=1 p j (M (k) j,p ) 2 K k=1 p j (M (k) p, j ) 2         1 4 . ( 33
)
Performing the N(N + 1)/2 successive updates using the whole set of (i, j) couples is called a sweep. Of course, this op-200 timization scheme is not equivalent to simultaneously estimate all the entries of L. As a consequence, we do not expect that at the end of the first sweep the matrix set {M (k) } is symmetrized. That is why, the sweeping procedure has to be repeated such as one iteration of the symmetrization algorithm corresponds to a 205 sweep. The algorithm is stopped when a suitable stopping criterion reaches a defined threshold or when a maximum number of iterations is reached. A classical stopping criterion is given by the relative decrease of C between two successive iterations.

It is worth mentioning that (i.) The identity matrix belongs to the set of elementary lower triangular matrices and to the set of elementary diagonal matrices. This ensures that criterion C cannot increase after an update and thus after a sweep. (ii.) Elementary lower triangular matrices are invertible by construction, meaning that update ( 21) is always possible. (iii.) During the symmetrization process the numerator and the denominator in equation [START_REF] Mosesov | Adaptive non-negative least squares with applications to non-negative matrix factorization[END_REF] tend to the same value. If the numerator and/or the denominator is equal to 0, one can impose a ( j) = 1 and let the algorithm run.

Computational complexity analysis

Although the ultimate goal of comparing JEVD or CPD techniques is to evaluate the quality of methods as reflected by the computed matrices, it is also interesting to assess the numerical complexity of these methods. Complexity is defined here as the number of floating point operations required to execute an algorithm (flops). A flop corresponds to a multiplication followed by an addition. But, in practice, only the number of multiplications is considered since, most of the time, there are about as many (and slightly more) multiplications as additions. In order to simplify the expressions, the complexity is generally approximated by its asymptotic limit, as the size of the problem tends to infinity. We shall subsequently denote, with some small abuse of notation, the equivalence between two strictly positive functions f and g:

f (x) = O(g(x)) or g(x) = O( f (x)) (34)
if and only if the ratio f (x)/g(x) tends to 1 as x → +∞. In practice, knowing whether an algorithm is computationally heavy is as important as knowing its performance in terms of SNR. This section first addresses the complexity of some JEVD methods needed by the DIAG algorithm to solve the CPD prob-225 lem. Next the numerical complexity of DIAG is given. In order to have a comparative reference, we provide the complexity of the alternating FNNLS (A-FNNLS) technique proposed in [START_REF] Bro | A fast non-negativity-constrained least squares algorithm[END_REF] and implemented in the N-way toolbox [START_REF] Bro | Parafac, tutorial and applications[END_REF][START_REF]The n-way toolbox for matlab[END_REF], which is used in section 5.

Complexity of JEVD methods

The costliest operation of JEVD algorithms is the multiplicative update step of the matrix set defined in equation [START_REF] Cardoso | Jacobi angles for simultaneous diagonalization[END_REF]. Now the numerical complexity of this step strongly depends of the structure of the updating matrix.

235

Regarding the JSJD 2 method, one sweep of the JS step involves N(N -1)/2 elementary lower triangular matrices and N elementary diagonal matrices as updating matrices. The numerical complexity of the JS step is then O(KN 3 ). In a similar way, the JD step based on Pham's algorithm involves Givens 240 matrices. Therefore its numerical complexity is also O(KN 3 ). Consequently, our JEVD method costs O(KN 3 ) flops.

Note that the numerical complexity of classical JEVD methods such as JUST [START_REF] Iferroudjene | A new jacobi-like method for joint diagonalization of arbitrary non-defective matrices[END_REF], SH-RT [START_REF] Fu | Simultaneous diagonalization with similarity transforma-720 tion for non-defective matrices[END_REF], JDTM [START_REF] Luciani | Canonical polyadic decomposition based on joint eigenvalue decomposition[END_REF], JAPAM [START_REF] André | A new class of block coordinate algorithms for the joint eigenvalue decomposition of complex matrices[END_REF], SJDTE [START_REF] André | Joint eigenvalue decomposition algorithms based on first-order taylor expansion[END_REF] and JET-U [START_REF] Luciani | Joint eigenvalue decomposition of non-defective matrices based on the lu factorization with application to ICA[END_REF] is O(KN 3 ) too.

Complexity of two CPD methods: DIAG and A-FNNLS

The computational complexity of DIAG is clearly dominated by the three following computations. First, the truncated SVD of the unfolding matrix of size P × QR requires 2QRP 2 + 5N 2 (P + QR) -2(N 3 + P 3 )/3 multiplications by as-250 suming that QR > P. Then, the computation of the M (k) matrices needs approximately (NR) 2 Q additional multiplications. Finally the cost of the JEVD procedure is O(R 2 N 3 ) according to section 4.1 since K = R(R -1). Additional computations can be neglected. Now, let's consider a low-rank CPD, 255 i.e. N min(P, Q, R), as observed in most of the applications such as the fluorescence one. Thus, the numerical complexity of DIAG is O(P 2 (P 2 + PN 2 + N 3 ).

As far as the A-FNNLS algorithm is concerned, the estimation of the factors is performed by solving alternatively Non-260 Negative Least Squares (NNLS) problems. The first widely used algorithm for solving the NNLS problem is an active-set method published in [START_REF] Lawson | Solving least squares problems[END_REF]. A modification of the latter was proposed in [START_REF] Bro | A fast non-negativity-constrained least squares algorithm[END_REF] giving rise to the FNNLS (Fast NNLS) method. The effective benefit of FNNLS is from forming the normal 265 equations and recognizing that the normal equations on the passive set are simply a subset of these. The computational advantages, then, are only applicable for overdetermined systems, which is generally the case for the CPD problems encountered in applications. Without utilizing the normal equa-270 tions, the original NNLS algorithm has O(MN 2 ) complexity for each iteration where M denotes the second dimension of the unfolding matrix used to estimate some factors (e.g. QR for the unfolding matrix (3)). FNNLS, however, only has an initialization complexity of O(MN 2 ), followed by O(N 3 ) complexity in each iteration. One may recognize that the use of the normal equations comes at the price of numerical stability. This approach, therefore, may not be appropriate when presented with an ill-conditioned problem [START_REF] Mosesov | Adaptive non-negative least squares with applications to non-negative matrix factorization[END_REF]. Under the same assumptions as for DIAG, the numerical complexity of A-FNNLS is

O(P(P 2 N 2 + N 3 )).
This asymptotic analysis shows that when the tensor dimensions of the low-rank CPD tend to infinity, the DIAG method requires more flops than A-FNNLS due to the truncated SVD (O(P 4 ) flops) computed at the beginning of DIAG. Nevertheless, for tensors of smaller size such as that considered in the following section, the truncated SVD is rarely the most costly step of DIAG. Indeed, for such scenarios, by considering cubic tensors, the term 2P 2 N 3 (I S D + 2I JD ) is generally dominant in the complexity of DIAG, where I S D and I JD denote the numbers of iterations of the two steps of JSJD 2 . As far as A-FNNLS is concerned, its dominant term is I A P 3 N 2 , where I A denotes the number of iterations (all factors have been estimated once using the FNNLS method during one iteration). Since the number I A of iterations is generally much larger than I S D + 2I JD and since P > N, the term I A P 3 N 2 appears to be larger than 2P 2 N 3 (I S D + 2I JD ) in such a context. This will be confirmed in the following section by means of the computed CPU time: DIAG will appear to be less time-consuming than A-FNNLS for the considered fluorescence spectroscopy scenarios.

Application to the Canonical Polyadic Decomposition of

Fluorescence spectra

CPD of fluorescence excitation emission matrices

In this section, we show how the proposed approach can be applied in the context of fluorescence spectroscopy.

Fluorescence is defined as the emission of light (fluorescence or emitted light) by a substance (molecule or atom) that has previously absorbed another light (excitation light). A fluorescing chemical component is called a fluorophore. The fluorescence intensity of a diluted solution of several fluorophores varies with the wavelength of the excitation light, the wavelength of the fluorescence light and the fluorophore concentrations. Spectrofluorimeters allow to measure the fluorescence intensity of a diluted solution at various excitation-emission wavelengths couples, providing a data matrix that characterizes the solution. Such a matrix is called Fluorescent Excitation-Emission Matrix (FEEM) or fluorescence 3D spectra.

In practice, chemists often have at their disposal a set of fluorescent solutions from which they measure a set of FEEM. We denote X r the FEEM characterizing the r th chemical solution. Each FEEM can be seen as a frontal slice of a third order tensor. Therefore, measured fluorescence intensities can be stored in a tensor T so that T p,q,r is the fluorescence intensity of solution r at the excitation wavelength p and emission wavelength q. Each solution can be seen as a mixture of N fluorophores.

Then, assuming that non-linear effects such as Rayleigh and Raman scattering or inner filter effects have been corrected, it has been showed in [START_REF] Lakowicz | Principles of Fluorescence Spectroscopy[END_REF] that

T p,q,r N n=1 A p,n B q,n C r,n (35)
where

• A p,n is the capacity of the fluorophore n to absorb light at excitation wavelength number p. Column n of A defines 320 the excitation spectra of the fluorophore n.

• B q,n is the capacity of the fluorophore n to emit light at emission wavelength number q. Column n of B defines the emission spectra of the fluorophore n.

• C r,n is proportional to the concentration of the fluorophore 325 n in solution r. Column n of C defines the concentration profile of the fluorophore n.

Here, we also impose that

||a .,n || 2 = ||b .,n || 2 = 1, ∀n. (36)
As a consequence, c .,n represents the contributions of fluorophore n to the different mixtures. We then define

C n = ||c .,n || 2 , ∀n (37)
as the contribution of fluorophore n to the fluorescence tensor T . In order to gather all the spectral information about fluorophore n in an unique spectral signature, we define its individual 3D spectrum as the rank one matrix a .,n b .,n T . Finally, the CPD model considers that each FEEM X r is modeled as X r = Adiag{c r,. }B T , ∀r.

(38)

Simulations on synthetic data sets

We first propose to evaluate our approach by using several synthetic data sets designed to cover some of the main issues 330 encountered with real unknown fluorescence data:

1. some fluorophores can have similar excitation or emission spectra, 2. the actual number of fluorophores is unknown and it is overestimated (overfactoring). We first build four synthetic excitation and emission spectra defining matrices A and B in equation [START_REF] Luciani | How to correct inner filter effects altering 3d fluorescence spectra by using a mirrored cell[END_REF]. For all data sets, excitation and emission spectra are of size P = 41 and Q = 71 respectively and are norm 1. In the proposed numerical simulations, excitation and emission spectra are unimodal spectra similar to those presented in figure 1. Entries of C are randomly drawn from a uniform distribution. Columns of C are also normalized so that all fluorophores have unit contributions. Then, R synthetic noise-free FEEM are built one after another according to equation [START_REF] Luciani | How to correct inner filter effects altering 3d fluorescence spectra by using a mirrored cell[END_REF]. In order to simulate model errors due to the FEEM acquisition and signal processing, we build from each noise-free FEEM the corresponding slice of the noisy tensor T .,.r as follows:

T .,.r = A{G r X r } (39)
where X r and G r are two random matrices of size P×Q that simulate the photon detection counting and the possible variations 345 of the photomultiplier gain respectively. is the Hadamard product and A is a smoothing operator. More precisely:

• Entry (p, q) of Xr is randomly drawn from a Poisson distribution of mean (X r ) p,q : ∀ (p, q), ( Xr ) p,q ∼ Pois((X r ) p,q ). It is noteworthy that before applying Poisson noise, the matrix X r is normalized in order to set its lowest value to one.

• Entries of G r are randomly drawn from a uniform distribution in the interval [1 -α; 1 + α] : ∀ (p, q), ( Xr ) p,q ∼ unif(1 -α, 1 + α). Hence, parameter α allow us to vary the amplitude of the multiplicative noise.

• Here A is a Savitzky-Golay filter. This filter is commonly used by spectrofluorometer manufacturer softwares.

Other non linear deviations such as diffusion scattering or inner filter effects are not considered here. This process is repeated 360 100 times for six different values of the noise amplitude α from 0.05 to 0.5. Note that α = 0.5 means that the photomultiplier gain varies from 50% to 150% around its mean value (giving a relative standard deviation of 29%). Thus each data set consists in six subsets of 100 noisy tensors. As it will be explained in the dedicated subsections below, data sets will differ in the number of mixtures (R) and in the correlation coefficients between the factors of the two spectral modes. Now regarding the algorithms, we compare three different approaches. Each approach is represented by two algorithms.

370

The first approach resorts to the DIAG algorithm (introduced in section 2) to compute the CPD. Two implementations are considered. The first implementation, referred to as DIAG-JDJS2 , corresponds to the proposed JDJS 2 algorithm: after the symmetrization step, the symmetry of the matrices to diagonalize 375 is imposed by construction and the joint diagonalization step is performed thanks to Pham algorithm. Indeed, thanks to preliminary comparisons (results not shown), we obtained best results with this configuration. The second implementation corresponds to a reference JEVD algorithm: the JEVD step is com-380 puted with the SJDTE algorithm recently introduced in [START_REF] André | Joint eigenvalue decomposition algorithms based on first-order taylor expansion[END_REF]. This algorithm belongs to the family of JDTE algorithm. Its main advantage in the present context of small real matrices is its rapidity. This implementation is referred to as DIAG-SJDTE. Both JEVD algorithms are initialized with the identity 385 matrix. The JEVD step is stopped when a maximum number of iterations is reached or when the relative deviation between two successive values of the algorithms criterion falls under a defined threshold. In JDJS 2 , JS and JD steps are stopped when the relative deviation between two successive values of the cri-390 terion is lower than 10 -6 . Maximum number of iterations are set to 100 for both steps. The same stopping criterion is used for SJDTE but the maximum number of iterations is set to 50. Here, we have the prior knowledge that entries of the factor matrices are nonnegative. Therefore, after the JEVD step in the 395 DIAG procedure, we take A = [USW] + and C B = [VW -T ] + .

The second approach resorts to alternating CPD algorithms with non-negative constraints. The first algorithm of this approach is the A-FNNLS algorithm mentioned in section 4 of the N-way toolbox1 (parafac function) which is commonly used 400 to analyze fluorescent excitation-emission matrices. The second algorithm is the recent AO-ADMM algorithm [START_REF] Huang | A flexible and efficient algorithmic framework for constrained matrix and tensor factorization[END_REF]. This algorithm has proved to be very efficient to deal with overfactoring and correlated factors. We used the AO-ADMM implementation of Jeremy Cohen 2 [START_REF] Cohen | Environmental Multiway Data Mining, Theses[END_REF] (AOadmm function). A-405 FNNLS is initialized using the best-fitting model scheme implemented in the parafac function: algorithm is run from 5 different random starting values and keeps the starting values that has led to the best fit after few iterations. We keep the default values of the convergence criterion (10 -6 ) and the max-410 imal iteration number (2500). The starting values selected by the best-fitting model scheme are reused as starting values for the AO-ADMM algorithm. We keep the default settings (leastsquare loss, non-negative constraint, tolerance of 10 -7 ) of the AOadmm function but the maximal number of iterations is set to 415 1000.

The last approach consists in using the factors estimated by the algorithms of the first approach to initialize an algorithm of the second approach and then perform a very limited number of iterations as proposed in [START_REF] Luciani | Canonical polyadic decomposition based on joint eigenvalue decomposition[END_REF]. The main advantage of this hybrid approach is to improve the results obtained by the DIAG algorithm while keeping a limited computational cost. Here, we use the factors estimated by DIAG-JDJS 2 and DIAG-SJDTE as starting values and we perform a maximum of 10 iterations of AO-ADMM. These algorithms are referred to as JDJS 2 -AO-ADMM and SJDTE-AO-ADMM respectively.

In order to evaluate the robustness of each algorithm with respect to overfactoring, we compute CPD of rank 6 instead of 4. Then, the 6 computed factors are sorted according to their contributions (fluorophores contribution are defined in equation ( 37)). We keep the 4 factors with the highest contributions as the estimation of the actual factors. Scaling and permutation indeterminacy are then removed as explained in [START_REF] Comon | Tensor decompositions, alternating least squares and other tales[END_REF]. We denote the estimated factors linked to fluorophore n by â.,n , b.,n , ĉ.,n . Algorithms are then compared according three criteria:

435

• The mean least square error between the 3D spectra of the 4 fluorophores and their estimation denoted e AB :

e AB = 1 4 4 n=1 ||â .,n b.,n T -a .,n b .,n T || 2 ||a .,n b .,n T || 2 , (40)
• The mean least square error between the concentration profiles of the 4 fluorophores and their estimation denoted e C :

e C = 1 4 4 n=1 ||ĉ .,n -c .,n || 2 ||c .,n || 2 (41)
• The computation time t c .

All the following numerical simulations were performed using matlab R2017a on a linux system with an Intel Core i7-7700T CPU and 8GB RAM.

Synthetic data set 1 : small tensors and uncorrelated

cient between two columns of A or B is 0.63. Figure 2 shows the evolution of the average values of the three criteria with α.

Average values are computed from the 100 tensors of each subset. First of all, recall we are in an over-factoring situation with two extra-components out of 6. In this context, it is noteworthy that all the considered algorithms performed quite well with estimation errors falling under 12 %. Regarding both estimation error criteria (e AB and e C ) DIAG-JDJS 2 provides better results than DIAG-SJDTE, A-FNNLS and AO-ADMM whatever the considered value of α. Hybrid algorithms improve significantly DIAG results. Among them the best results are clearly achieved by JDJS 2 -AO-ADMM. Now regarding the computing time, DIAG-JDJS 2 and JDJS 2 -AO-ADMM are about two to seven times faster than A-FNNLS and AO-ADMM but slower than SJDTE based algorithms. The gap with SJDTE is mainly explained by the number of iterations computed during the symmetrization procedure.

Synthetic data set 2 : large tensors and uncorrelated spectra

We now show what happen on bigger data tensors. FEEM 465 dimensions used for data set 1 are common while the number of mixtures can vary a lot from one real word application to another. Therefore, here we use the same spectra as synthetic data set 1 but we now consider a larger number of mixtures by taking R = 200. The JEVD step in data set 1 involved 870 ma-470 trices of size 6. Thus, we had K/N = 145 and KN 3 2.10 5 .

Here, the situation is different since we have now : K = 39800, K/N 6633 and KN 3 9.10 6 . In other words, the numerical complexity has been multiplied by 45 whereas the size of the matrix set was already significantly greater than the matrix 475 size. We thus decimate the matrix set by a factor 10 (we leave out 90% of the matrices) in order to speed up the DIAG algorithm as explained in section 2. The same decimation will be applied to the next two synthetic data sets. Results are plotted on figure 3. These results are similar to data set 1. In particular, 480 JDJS 2 -AO-ADMM still provides the lowest estimation errors while keeping a very limited computing time. Indeed thanks to the decimation JDJS 2 -AO-ADMM is now about ten and fifteen times faster than A-FNNLS and AO-ADMM respectively. in term of estimation error. The correlations affect the estimation performances of JDJS 2 -AO-ADMM when α ≥ 0.4. In this range, A-FNNLS and AO-ADMM provide now the best results. Under this critical value JDJS 2 -AO-ADMM performs slightly better than A-FNNLS and AO-ADMM. It is worth mentioning 500 that at the lowest values of α (0.05 and 0.1) DIAG-JDJS 2 competes with those algorithms. Regarding the computing times, t c (JDJS 2 -AO-ADMM) is about 20 to 30 times lower than t c (A-FNNLS) and t c (AO-ADMM) whatever the value of α. This indicate that we could decrease the decimation factor for α ≥ 0.4 505 in order to improve the performances of JDJS 2 -AO-ADMM. with those of data set 3, it seems that two correlations in one mode (data set 3) affect more severely the proposed approach than one correlation in two modes (data set 4). Finally, as for data set 3, the main difference between the different approaches is the computation speed. Indeed, mean values of t c are 1.35s for JDJS 2 -AO-ADMM, 32.6s for AO-ADMM and 30.6s for A-FNNLS.

Synthetic data set 3 : large tensors and correlated exci-

Simulations on real data sets

Real data set 1

This data set is composed of 11 FEEM, measured from diluted mixtures of fluorescein, quinine sulfate and tryptophane. Fluorophore spectra and concentration profiles are plotted on figure 6. These FEEM were used as reference spectra (no inner filter effect) in a previous study [START_REF] Luciani | How to correct inner filter effects altering 3d fluorescence spectra by using a mirrored cell[END_REF]. Excitation an emission wavelengths vary from 270 to 550 nm with a step of 5 nm. The data set is thus a tensor of size (57 × 57 × 11). Rayleigh and Raman scatters were removed numerically from each measured FEEM using the method proposed by Zepp [START_REF] Zepp | Dissolved organic fluorophores in southeastern us coastal waters: correction method for eliminating rayleigh and raman scattering peaks in excitation-emission matrices[END_REF]. We referrer to the above mentioned paper for details about sample preparation and FEEM acquisition. We compare the results obtained from various CPD computed by the previous algorithms. Indeed, in order to evaluate the influence of possible over-factoring, outperform the alternating algorithms. Finally at F = 4, all algorithms fail. If we focus on e AB , JDJS 2 algorithm helps to improve DIAG performances. However, here, when the hybrid approach is used, we cannot discriminate between SJDTE and the proposed algorithm. Thus, these results mainly highlight the 555 ability of the DIAG and hybrid approaches to deal with overfactoring on a real data set. These algorithms are also clearly faster than A-FNNLS and AO-ADMM with random initialization in particular for F > 0.

Real data set 2 560

We performed the same comparison study on an other data set which is publicly available at http://www.models.life. ku.dk/Fluorescence [START_REF] Rinnan | Application of parafac on spectral data[END_REF][START_REF] Bro | Standard error of prediction for multilinear pls -2. practical implementation in fluorescence spectroscopy[END_REF]. This data set consists in 16 mixtures of resorcinol, tryptophane and tyrosine and corresponds to dataset #1 in the web page 3 . Fluorescence inten-565 sity was measured for 19 different excitation wavelengths and 3 We chose this data set because, according to the publisher of the data, in all other data sets one of the components "contains some impurities, and thus gives rise to an extra component" 136 different emission wavelengths. We removed the data corresponding to excitation wavelength shorter than 250 nm in order to focus on the main excitation peaks of the 3 fluorophores.We also removed the first sample because it corre-570 sponds to a pure fluorophore. Thus our final data set is a tensor of size 16 × 136 × 15. Rayleigh and Raman scatters were removed as explained previously. Fluorophore spectra and concentration profiles are plotted on figure 8. Fluorophore spectra are highly correlated both in excitation and emission mak-575 ing this data set quite difficult. Results are plotted on figure 9. First of all, performances of DIAG and hybrid approaches are remarkably stable with respect to the number of extracomponents, whatever the considered criterion. In contrast, A-FNNLS and AO-ADMM error plots are quite fluctuating and 580 A-FNNLS t c plot shows a clear increase of the computing time with the number of extra-components. Above all, DIAG-JDJS 2 and JDJS 2 -AO-ADMM clearly provide the lowest estimation errors whatever the number of extra-components. For instance at F = 0, e AB is about 0.04 for DIAG-JDJS 2 and JDJS 2 -AO-585 ADMM, about 0.18 for DIAG-SJDTE and SJDTE-AO-ADMM and about 0.3 for A-FNNLS and AO-ADMM. Finally, DIAG and hybrid algorithms are much faster than A-FNNLS and AO-ADMM (at F = 0 we have t c (JDJS 2 -AO-ADMM)=0.22s, t c (A-FNNLS)=1.6s and t c (AO-ADMM)=3.9s). 590

Conclusion

In this study, we have first recalled how the canonical polyadic decomposition of a multiway array can be computed by means of joint eigenvalue decomposition. At this occasion, we also explain how this kind of CPD algorithm can be eas-595 ily speed-up by matrix decimation in order to obtain fast algorithms. We have then showed how the JEVD can be rewritten into the simpler problem of orthogonal joint diagonalization of a set of positive definite matrices. This is done by joint symmetrization of the original matrix set. We have thus proposed 600 an ad hoc joint symmetrization algorithm based on a sweeping procedure and multiplicative updates. This allowed us to introduce an original JEVD algorithm called JDJS 2 based on a two step process: joint diagonalization after joint symmetrization. This approach can thus take benefit of efficient exist-605 ing OJD algorithms. We have then proposed to apply this approach to the CPD of data tensors containing real or synthetic fluorescence measurements by including JDJS 2 into a JEVDbased CPD algorithm. We call the resulting algorithm algorithm DIAG-JDJS 2 . Both synthetic and real data sets were built 610 to be challenging. Two main difficulties were considered: the number of component (i.e. the rank of the decomposition) is unknown and overestimated and several factors can be highly correlated. The proposed DIAG-JDJS 2 algorithm was compared to an other version of DIAG using a recent reference JEVD al-615 gorithm and with two reference CPD algorithms based on an alternating approach namely A-FNNLS and AO-ADMM. Our results first show that in most situations, the proposed algorithm significantly improves DIAG performances in term of factor estimation accuracy. Second, DIAG-JDJS 2 sometimes outperforms the A-FNNLS and AO-ADMM algorithms in the simplest situations (in terms of noise level and factor correlations) while having a computing time significantly lower. We then show that A-FNNLS and AO-ADMM are more robust to the noise level and correlations between factors. How-625 ever, we also show that DIAG-JDJS 2 provides a very good starting point for alternating algorithms. This hybrid approach not only achieves the best performances in terms of estimation accuracy (in most of the considered situations) but it is also very fast in comparison to a classical random initialization strategy 630 such as the best fitting scheme. This study also opens some perspectives. In particular, the proposed JEVD algorithm allows us to impose a non-negative constraint on the eigenvalues. This feature is not really exploited here. An extension to the complex case is also worth consid-635 ering. Regarding the CPD, here we use matrix decimation in the DIAG algorithm to build a smaller JEVD problem. Better selection strategies may be investigated. Finally it may be interesting to include JDJS 2 in other CPD algorithms based on JEVD. • Cannonical Polyadic Decomposition (CPD) also known as PARAFAC has become an important chemometric tool to decompose multiway data sets, notably in fluorescence spectroscopy.

• Joint eigenvalue decomposition and joint diagonalization are at the center of efficient CPD algorithms.

•

We propose an original joint eigenvalue decomposition algorithm based on joint symmetrization.

• This algorithm is succesfully applied to the CPD of synthetic fluorescence data sets.

• We show that our approach improves the estimation of the CPD factors in difficult but realistic cases (correlated factors and unknown number of factors.)

J o u r n a l P r e -p r o o f

  105

  160

  230

  245

335 5 . 2 . 1 .

 521 Data sets and algorithm settingsOur data sets correspond here to different synthetic fluorescence tensors of size (P × Q × R). These tensors are built from a noisy CPD model of rank 4. Hence, each tensor gathers simulated fluorescence measurements of R random mixtures of 4 340 synthetic fluorophores and each data set contains 100 realizations of these synthetic fluorescence tensors for each considered noise level. Synthetic data sets are built as follows.

Figure 1 :

 1 Figure 1: Synthetic excitation and emission spectra used to build data set 1.

  keep R = 200 but we modify the position (in the wavelength range) of the excitation spectra of fluorophore 2 and 4 in order to increase the overlap with the excitation spectra of fluorophore 1 and 3 (the shapes of the spectra remain un-490 changed). Correlation coefficients between fluorophores 1 and 2 on one side and between fluorophores 3 and 4 on the other side are now equal to 0.98. Results are plotted on figure 4. If we consider the whole range of α values, JDJS 2 -AO-ADMM, A-FNNLS and AO-ADMM outperform the three other algorithms 495

5. 2 . 5 .Figure 2 :

 252 Figure 2: Synthetic data set 1: evolution of the average values of the three comparison criteria with respect to the noise amplitude.

Figure 3 :Figure 4 :

 34 Figure 3: Synthetic data set 2: evolution of the average values of the three comparison criteria with respect to the noise amplitude.

Figure 5 :

 5 Figure 5: Synthetic data set 4: evolution of the average values of the three comparison criteria with respect to the noise amplitude.

Figure 6 :Figure 7 :

 67 Figure 6: Real data set 1: Fluorophore spectra and concentration profiles.

Figure 8 :

 8 Figure 8: Real data set 2: Fluorophore spectra and concentration profiles.

Figure 9 :

 9 Figure 9: Real data set 2: evolution of the average values of the three comparison criteria with respect to the number of extra-components.

http://www.models.life.ku.dk/nwaytoolbox

http://www.gipsa-lab.grenoble-inp.fr/ pierre.comon/TensorPackage/tensorPackage.htm

spectra For this first data set, we consider a (quite) small number of mixtures by taking R = 30. Synthetic spectra used to build the data set are plotted on figure 1. These spectra are moderately correlated. More precisely the higher correlation coeffi-

Declaration of interests

☐ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: