Pierre Comon

Xavier Luciani

André L F De Almeida

Tensor Decompositions, Alternating Least Squares and other Tales

come   L'archive ouverte pluridisciplinaire

Introduction

Richard Harshman liked to explain Multi-Way Factor Analysis (MFA) as one tells a story: with words, sentences, appealing for intuition, and with few formulas. Sometimes, the story turned to a tale, which required the belief of participants, because of the lack of proof of some strange -but correct-results.

MFA has been seen first as a mere nonlinear Least Squares problem, with a simple objective criterion. In fact, the objective is a polynomial function of many variables, where some separate. One could think that this kind of objective is easy because smooth, and even infinitely differentiable. But it turns out that things are much more complicated than they may appear to be at first glance. Nevertheless, the Alternating Least Squares (ALS) algorithm has been mostly utilized to address this minimization problem, because of its programming simplicity. This should not hide the inherently complicated theory that lies behind the optimization problem.

Note that the ALS algorithm has been the subject of much older tales in the past. In fact, it can be seen as a particular instance of the nonlinear Least Squares problem addressed in [START_REF] Golub | Differentiation of pseudo-inverses, separable nonlinear least square problems and other tales[END_REF], where variables separate. The latter analysis gave rise in particular to the Variable Projection algorithm, developed by a great figure of Numerical Analysis, who also passed away recently, barely two months before Richard.

Tensors play a wider and wider role in numerous application areas, much beyond Chemometrics. Among many others, one can mention Signal Processing for Telecommunications [START_REF] Sidiropoulos | Blind PARAFAC receivers for DS-CDMA systems[END_REF] [START_REF] Delathauwer | Tensor-based techniques for the blind separation of DS-CDMA signals[END_REF] [START_REF] Almeida | PARAFAC-based unified tensor modeling for wireless communication systems with application to blind multiuser equalization[END_REF], and Arithmetic Complexity [START_REF] Kruskal | Three-way arrays: Rank and uniqueness of trilinear decompositions[END_REF] [62] [START_REF] Ürgisser | Algebraic Complexity Theory[END_REF], which are two quite different frameworks. One of the reasons of this success is that MFA can often replace Principal Component Analysis (PCA), when available data measurements can be arranged in a meaningful tensor form [START_REF] Smilde | Multi-Way Analysis[END_REF]. When this is not the case, that is, when the observation diversity is not sufficient (in particular when a 3rd order tensor has proportional matrix slices), one can sometimes resort to High-Order Statistics (HOS), allowing to build symmetric tensors of arbitrarily large order from the data [START_REF] Comon | Blind identification of under-determined mixtures based on the characteristic function[END_REF].

In most of these applications, the decomposition of a tensor into a sum of rank-1 terms (cf. Section 2) allows to solve the problem, provided the decomposition can be shown to be essentially unique (i.e. unique up to scale and permutation). Necessary and sufficient conditions ensuring existence and essential uniqueness will be surveyed in Section 2.3. The main difficulty stems from the fact that actual tensors may not exactly satisfy the expected model, so that the problem is eventually an approximation rather than an exact decomposition.

Richard Harshman was very much attached to the efficiency of numerical algorithms when used to process actual sets of data. He pointed out already in 1989 [START_REF] Kruskal | How 3-MFA data can cause degenerate Parafac solutions, among other relationships[END_REF] that one cause of slow convergence (or lack of convergence) can be attributed to "degeneracies" of the tensor (cf. Section 2.2). Richard Harshman proposed the following classification:

1. Bottleneck. A bottleneck occurs when two or more factors in a mode are almost collinear [START_REF] Rajih | Enhanced line search : A novel method to accelerate PARAFAC[END_REF].

2.

Swamp. If a bottleneck exists in all the modes, then we have what has been called by R. Harshman and others a "swamp" [START_REF] Mitchell | Slowly converging Parafac sequences: Swamps and two-factor degeneracies[END_REF] [52] [START_REF] Rajih | Enhanced line search : A novel method to accelerate PARAFAC[END_REF].

3. CP-degeneracies may be seen as particular cases of swamps, where some factors diverge and at the same time tend to cancel each other as the goodness of fit progresses [START_REF] Kruskal | How 3-MFA data can cause degenerate Parafac solutions, among other relationships[END_REF] [47] [START_REF] Harshman | An annoted bibliogrtaphy of articles on degenerate solutions or decompositions[END_REF].

CP-degeneracies occur when one attempts to approximate a tensor by another of lower rank, causing two or more factors to tend to infinity, and at the same time to almost cancel each other, giving birth to a tensor of higher rank (cf. Section 2.2). The mechanism of these "degeneracies" is now rather well known [START_REF] Paatero | Construction and analysis of degenerate Parafac models[END_REF] [63] [START_REF] Stegeman | Degeneracy in Candecomp/Parafac explained for p × p × 2 arrays of rank p + 1 or higher[END_REF] [59] [START_REF] Comon | Symmetric tensors and symmetric tensor rank[END_REF], and has been recently better understood [START_REF] De | Tensor rank and the ill-posedness of the best low-rank approximation problem[END_REF] [START_REF] Krijnen | On the nonexistence of optimal solutions and the occurrence of degeneracy in the Candecomp/Parafac model[END_REF]. If one dimension is equal to 2, special results can be obtained by viewing third order tensors as matrix pencils [START_REF] Stegeman | Low-rank approximation of generic pxqx2 arrays and diverging components in the Candecomp/Parafac model[END_REF]. According to [START_REF] De | Tensor rank and the ill-posedness of the best low-rank approximation problem[END_REF], the first example of a rank-r tensor sequence converging to a tensor of rank strictly larger than r was exhibited by D.Bini as early as the seventies [START_REF] Bini | Border rank of a p × q × 2 tensor and the optimal approximation of a pair of bilinear forms[END_REF]. CP-degeneracies can be avoided if the set in which the best approximate is sought is closed. One solution is thus to define a closed subset of the set of lower-rank tensors; another is to define a superset, yielding the notion of border rank [START_REF] Strassen | Rank and optimal computation of generic tensors[END_REF] [4] [START_REF] Ürgisser | Algebraic Complexity Theory[END_REF].

On the other hand, swamps may exist even if the problem is well posed. For example, Paatero reports in [START_REF] Paatero | Construction and analysis of degenerate Parafac models[END_REF] precise examples where the trajectory from a well conditioned initial point to a well-defined solution has to go around a region of higher rank, hence a "swamp". In such examples, the convergence to the minimum can be slowed down or even compromised.

One solution against bottlenecks that we developed together with R. Harshman just before he passed away was the Enhanced Line Search (ELS) principle [START_REF] Rajih | Enhanced line search : A novel method to accelerate PARAFAC[END_REF], which we subsequently describe. This principle, when applied to the Alternating Least Squares (ALS) algorithm, has been shown in [START_REF] Rajih | Enhanced line search : A novel method to accelerate PARAFAC[END_REF] to improve significantly its performances in the sense that it would decrease the risk to terminate in a local minimum, and more importantly also need an often much smaller number of iterations. Also note that an extension to the complex case has been recently proposed [START_REF] Nion | An enhanced line search scheme for complex-valued tensor decompositions. application in DS-CDMA[END_REF].

However, the additional computational cost per iteration was shown to be negligible only when the rank was sufficiently large compared to the dimensions; see condition [START_REF] Carroll | Analysis of individual differences in multidimensional scaling via n-way generalization of Eckart-Young decomposition[END_REF]. Nevertheless, we subsequently show that this condition is not too much restrictive, since one can always run a prior dimensional reduction of the tensor to be decomposed with the help of a HOSVD truncation, which is a compression means that many users practice for several years [START_REF] Delathauwer | A multilinear singular value decomposition[END_REF] [55] [START_REF] Bro | Improving the speed of multiway algorithms. part ii: Compression[END_REF].

The goal of the paper is two-fold. First we give a summary of the state of the art, and make the distinction between what is known but attributed to usual practice, and thus belongs to the world of conjectures, and what has been rigorously proved. Second, we provide details on numerical complexities and show that ELS can be useful and efficient in a number of practical situations.

Section 2 defines a few notations and addresses the problems of existence, uniqueness, and genericity of tensor decompositions. Section 3 gathers precise update expressions appearing in iterative algorithms used in Section 5. We shall concentrate on Alternating Least Squares (ALS), Gradient descent with variable step, and Levenberg-Marquardt algorithms, and on their versions improved with the help of the Enhanced Line Search (ELS) feature, with or without compression, the goal being to face bottlenecks more efficiently. Section 4 gives orders of magnitude of their computational complexities, for subsequent comparisions in Section 5. In memory of Richard, we shall attempt to give the flavor of complicated concepts with simple words.

Tensor decomposition 2.1 Tensor rank and the CP decomposition

Tensor spaces. Let V (1) , V (2) , V (3) be three vector spaces of dimension I, J and K respectively. An element of the vector space V (1) ⊗ V (2) ⊗ V (3) , where ⊗ denotes the outer (tensor) product, is called a tensor. Let us choose a basis {e

(ℓ)
i } in each of the three vector spaces V (ℓ) . Then any tensor T of that vector space of dimension IJK has coordinates T ijk defined by the relation ℓ) , for some given matrices Q (ℓ) , 1 ≤ ℓ ≤ d then the new coordinates T ′ pqr of tensor T may be expressed as a function of the original ones as

T = I i=1 J j=1 K k=1 T ijk e (1) i ⊗ e (2) j ⊗ e (3) k Multi-linearity. If a change of basis is performed in space V (ℓ) , defined by e (ℓ) i = Q (ℓ) e ′ i (
T ′ pqr = I i=1 J j=1 K k=1 Q (1) pi Q (2) qj Q (3) rk T ijk
This is a direct consequence of the construction of tensors spaces, which can be readily obtained by merely pugging the exporession of e ′ i (ℓ) in that of T . By convention, this multilinear relation between arrays of coordinates is chosen to be written as [START_REF] De | Tensor rank and the ill-posedness of the best low-rank approximation problem[END_REF] [14]:

T ′ = (Q (1) , Q (2) , Q (3) ) • T CP decomposition. Any tensor T admits a decomposition into a sum of rank-1 tensors. This decomposition takes the form below in the case of a 3rd order tensor:

T = F p=1 a(p) ⊗ b(p) ⊗ c(p) (1)
where a(p), b(p) and c(p) denote vectors belonging to spaces V (1) , V (2) , V (3) . Denote A ip (resp. B jp and C kp ) the coordinates of vector a(p) in basis {e

(1) i , 1 ≤ i ≤ I}, (resp. b(p) in {e (2) j , 1 ≤ j ≤ J} and c(p) in {e (3) k , 1 ≤ k ≤ K}).
Then one can rewrite this decomposition as

T = F p=1 ( I i=1 A ip e (1) i ) ⊗( J j=1 B jp e (2) j ) ⊗( K k=1 C kp e (3) k ) or equivalently T = ijk ( F p=1 A ip B jp C kp ) e (1) i ⊗ e (2)
j ⊗ e

(3) k which shows that arrays of coordinates are related by

T ijk = F p=1 A ip B jp C kp (2)
One may see this equation as a consequence of the multi-linearity property. In fact, if one denotes by I the F × F × F three-way diagonal array, having ones on its diagonal, then (2) can be rewritten as

T = (A, B, C) • I
Note that, contrary to (1), equation ( 2) is basis-dependent.

Tensor rank. The tensor rank is the minimal number F of rank-1 terms such that the equality (1), or equivalently (2), holds true. Tensor rank always exists and is well defined. This decomposition can be traced back to the previous century with the works of Hitchcock [START_REF] Hitchcock | Multiple invariants and generalized rank of a p-way matrix or tensor[END_REF]. It has been later independently referred to as the Canonical Decomposition (Candecomp) [START_REF] Carroll | Analysis of individual differences in multidimensional scaling via n-way generalization of Eckart-Young decomposition[END_REF] or Parallel Factor (Parafac) [29] [31]. We shall simply refer to it via the acronym CP, as usual in this journal.

Points of terminology. The order of a tensor is the number of tensor products appearing in its definition, or equivalently, the number of its ways, i.e. the number of indices in its array of coordinates [START_REF] Alexander | La methode d'Horace eclatee: application a l'interpolation en degre quatre[END_REF]. Hence a tensor of order d has d dimensions.

In Multi-Way Factor Analysis, the columns a(p) of matrix A are called factors [57] [34]. Another terminology that is not widely assumed is that of modes [START_REF] Kroonenberg | Three mode Principal Component Analysis[END_REF]. A mode actually corresponds to one of the d linear spaces V (ℓ) . For instance, according to the classification sketched in Section 1, a bottleneck occurs in the second mode if two columns of matrix B are almost collinear.

Existence

In practice, one often prefers to fit a multi-linear model of lower rank, F < rank{T }, fixed in advance, so that we have to deal with an approximation problem. More precisely, it is aimed at minimizing an objective function of the form

Υ(A, B, C) = ||T - F p=1 a(p) ⊗ b(p) ⊗ c(p)|| 2 (3)
for a given tensor norm defined in V (1) ⊗ V (2) ⊗ V (3) . In general, the Euclidean norm is used. Unfortunately, the approximation problem is not always well posed. In fact, as pointed out in early works of R.Harshman [START_REF] Kruskal | How 3-MFA data can cause degenerate Parafac solutions, among other relationships[END_REF] [30], the infimum may never be reached, for it can be of rank higher than F .

CP-degeneracies. According to the classification proposed by Richard Harshman and sketched in Section 1, CP-degeneracies can occur when one attempts to approximate a tensor by another of lower rank. Let us explain this in this section, and start with an example.

Example. Several examples are given in [56] [14], and in particular the case of a rank-2 tensor sequence converging to a rank-4 tensor. Here, we give a slightly more tricky example. Let x, y and z be three linearly independent vectors, in R 3 for instance. Then the following sequence of rank-3 symmetric tensors

T (n) = n 2 (x + 1 n 2 y + 1 n z) ⊗ 3 + n 2 (x + 1 n 2 y - 1 n z) ⊗ 3 -2 n 2 x ⊗ 3
converges towards the tensor below, which may be proved to be of rank 5 (despite its 6 terms) as n tends to infinity:

T (∞) = x ⊗ x ⊗ y + x ⊗ y ⊗ x + y ⊗ x ⊗ x + x ⊗ z ⊗ z + z ⊗ x ⊗ z + z ⊗ z ⊗ x
Rank 5 is the maximal rank that can achieve third order tensors of dimension 3. To conclude, tensor T (∞) can be approximated arbitrarily well by rank-3 tensors T (n), but the limit is never reached since it has rank 5. Note that this particular example cannot be directly generated by Leibniz tensors described in [START_REF] De | Tensor rank and the ill-posedness of the best low-rank approximation problem[END_REF], which could have been used equally well for our demonstration (note however that T (∞) can probably be expressed as a Leibnitz operator after a multilinear change of coordinates).

Classification of third order tensors of dimension 2. Contrary to symmetric tensors in the complex field, whose properties have been studied in depth [START_REF] Brambilla | On the Alexander-Hirschowitz theorem[END_REF] [2] [20] [START_REF] Comon | Symmetric tensors and symmetric tensor rank[END_REF] (see also Section 2.3) and utilized in various applications including Signal Processing [START_REF] Comon | Decomposition of quantics in sums of powers of linear forms[END_REF], much less is known concerning non symmetric tensors [START_REF] Strassen | Rank and optimal computation of generic tensors[END_REF] [8], and even less in the real field. However, the geometry of 2-dimensional 3rd order asymmetric tensors is fully known (symmetric or not, real or not). Since the works of Bini [START_REF] Bini | Border rank of a p × q × 2 tensor and the optimal approximation of a pair of bilinear forms[END_REF], other authors have shed light on the subject [START_REF] Berge | Kruskal's polynomial for 2 × 2 × 2 arrays and a generalization to 2 × n × n arrays[END_REF] [58] [START_REF] Paatero | Construction and analysis of degenerate Parafac models[END_REF], and all the proofs can now be found in [26] [56]. Since the case of 2-dimensional tensors is rather well understood, it allows to figure out what probably happens in larger dimensions, and is worth mentioning. The decomposition of 2 × 2 × 2 tensors can be computed with the help of matrix pencils if one of the 2d = 6 matrix slices is invertible (one can show that if none of them is invertible, then the tensor is necessarily of rank 1). But hyperdeterminants can avoid this additional assumption [START_REF] De | Tensor rank and the ill-posedness of the best low-rank approximation problem[END_REF], which may raise problems in higher dimensions. Instead of considering the characteristic polynomial of the matrix pencil, det(T 1 + λ T 2 ), one can work in the projective space by considering the homogeneous polynomial p(λ 1 , λ 2 ) = det(λ 1 T 1 + λ 2 T 2 ). This does not need invertibility of one matrix slice. It can be shown [START_REF] De | Tensor rank and the ill-posedness of the best low-rank approximation problem[END_REF] to take the simple form:

p(λ 1 , λ 2 ) = λ 2 1 det T 1 + λ 1 λ 2 2 (det(T 1 + T 2 ) -det(T 1 -T 2 )) + λ 2 2 det T 2
It turns out that the discriminant ∆(T 1 , T 2 ) of polynomial p(λ 1 , λ 2 ) is nothing else but Kruskal's polynomial mentioned in [START_REF] Berge | Kruskal's polynomial for 2 × 2 × 2 arrays and a generalization to 2 × n × n arrays[END_REF] [58] [START_REF] Paatero | Construction and analysis of degenerate Parafac models[END_REF], which itself is due to Cayley [START_REF] Cayley | On the theory of linear transformation[END_REF]. As pointed out in [START_REF] De | Tensor rank and the ill-posedness of the best low-rank approximation problem[END_REF], the sign of this discriminant is invariant under changes of bases, and tells us about the rank of T : if ∆ < 0, rank{T } = 3, and if ∆ > 0, rank{T } = 2. Thus the variety defined by ∆ = 0 partitions the space into two parts, which hence have a non zero volume. The closed set of tensors of rank 1 lies on this hypersurface ∆ = 0. But two other rare cases of tensors of rank 2 and 3 share this hypersurface: in particular, this is where tensors of rank 3 that are the limit of rank-2 tensors are located.

How to avoid degeneracies. In [START_REF] De | Tensor rank and the ill-posedness of the best low-rank approximation problem[END_REF] [58] [START_REF] Krijnen | On the nonexistence of optimal solutions and the occurrence of degeneracy in the Candecomp/Parafac model[END_REF], it has been shown that CPdegeneracies necessarily involve at least two rank-1 terms tending to infinity. In addition, if all rank-1 terms tend to infinity, then the rank jump is the largest possible. In our example, all the three tend to infinity, and almost cancel with each other. This is the residual that gives rise to a tensor of higher rank.

If cancellation between diverging rank-1 tensors cannot occur, then neither CPdegeneracy can. This is in fact what happens with tensors with positive entries. Thus, imposing the non negativity during the execution of successive iterations is one way to prevent CP-degeneracy, if tensors with non negative entries are concerned [START_REF] Lim | Nonnegative approximations of nonnegative tensors[END_REF].

In order to guarantee the existence of a minimum, the subset of tensors of rank F must be closed, which is not the case of tensors with free entries, except if F ≤ 1 [START_REF] Comon | Symmetric tensors and symmetric tensor rank[END_REF]. Should this not be the case, one can perform the minimization either over a closed subset (e.g. the cone of tensors with positive entries [START_REF] Lim | Nonnegative approximations of nonnegative tensors[END_REF]) or over a closed superset. The most obvious way to define this superset is to define the smallest closed subset containing all tensors of rank F , i.e. their adherence; such a superset is that of tensors with border rank at most F [62] [4] [START_REF] Ürgisser | Algebraic Complexity Theory[END_REF]. Another idea proposed in [START_REF] Navasca | Swamp reducing technique for tensor decomposition[END_REF] is to add a regularization term in the objective function, which avoids the permutation-scale ambiguities to change when iterations progress. The effect is that swamps are better coped with.

However, no solution is entirely satisfactory, as pointed out in [START_REF] De | Tensor rank and the ill-posedness of the best low-rank approximation problem[END_REF], so that this problem remains open in most practical cases.

Uniqueness

The question that can be raised is then: why don't we calculate the CP decomposition for F = rank{T }? The answer generally given is that the rank of T is unknown. This answer is not satisfactory. By the way, if this were the actual reason, one could choose F to be larger than rank{T }, e.g. an upper bound. This choice is not made, actually for uniqueness reasons. In fact, we recall hereafter that tensors of large rank do not admit a unique CP decomposition.

Typical and generic ranks. Suppose that the entries of a 3-way array are drawn randomly according to a continuous probability distribution. Then the typical ranks are those that one will encounter with non zero probability. If only one typical rank exists, it will occur with probability 1 and is called generic.

For instance, in the 2 × 2 × 2 case addressed in a previous paragraph, there are two typical ranks: rank 2 and rank 3. In the complex field however, there is only one: all 2 × 2 × 2 tensors generically have rank 2.

In fact the rank of a tensor obtained from real measurements is typical with probability 1 because of the presence of measurement errors with a continuous probability distribution. And we now show that the smallest typical rank is known, for any order and dimensions.

With the help of the so-called Terracini's lemma [START_REF] Terracini | Sulla rappresentazione delle forme quaternarie mediante somme di potenze di forme lineari[END_REF], sometimes attributed to Lasker [START_REF] Lasker | Kanonische formen[END_REF] [START_REF] Comon | Decomposition of quantics in sums of powers of linear forms[END_REF], it is possible to compute the smallest typical rank for any given order and dimensions, since it coincides with the generic rank (if the decomposition were allowed to be computed in the complex field) [START_REF] Mella | Singularities of linear systems and the Waring problem[END_REF] [18] [START_REF] Brambilla | On the Alexander-Hirschowitz theorem[END_REF]. The principle consist of computing the rank of the Jacobian of the map associating the triplet (A, B, C) to tensor T defined in (2) [START_REF] Comon | Generic and typical ranks of three-way arrays[END_REF]. The same principle extends to more general tensor decompositions [START_REF] Comon | Generic and typical ranks of multi-way arrays[END_REF]. The resulting ranks are reported in Tables 1 and2 One could believe that this generic rank could be computed by just counting the number of degrees of freedom in both sides of decomposition [START_REF] Alexander | La methode d'Horace eclatee: application a l'interpolation en degre quatre[END_REF]. If we do this job we find IJK free parameters in the lhs, and F (I + J + K -2) in the rhs. So we can say that the "expected rank" of a generic tensor would be

F = IJK I + J + K -2 (4)
Now notice that the values given by ( 4) do not coincide with those given by Tables 1 and2; generic ranks that differ from expected values are underlined. These values are exceptions to the ceil rule. For instance, the row d = 3 in Table 2 shows that there is a single exception, which is consistent with [39, p.110]. We conjecture that the "ceil rule" is correct almost everywhere, up to a finite number of exceptions. In fact, this conjecture has been proved in the case of symmetric tensors [START_REF] Comon | Symmetric tensors and symmetric tensor rank[END_REF], thanks to Alexander-Hirschowitz theorem [START_REF] Alexander | La methode d'Horace eclatee: application a l'interpolation en degre quatre[END_REF] [5]. In the latter case, there are fewer exceptions than in the present case of tensors with free entries. Results of [START_REF] Abo | Induction for secant varieties of Segre varieties[END_REF] also tend to support this conjecture.

Another open problem is the following: what are the other typical ranks, when they exist? For what values of dimensions can one find several typical ranks, and indeed how many can there be? Up to now, only two typical ranks have been proved to exist for some given dimensions [START_REF] Berge | Partial uniqueness in CANDECOMP/PARAFAC[END_REF] [64] [START_REF] Berge | Symmetry transformations for square sliced three way arrays, with applications to their typical rank[END_REF] [66] [START_REF] De | Tensor rank and the ill-posedness of the best low-rank approximation problem[END_REF], not more. Number of decompositions. These tables not only give the most probable ranks, but allow to figure out when the decomposition may be unique. When the ratio in ( 2) is an integer, one can indeed hope that the solution is unique, since there are as many unknowns as free equations. This hope is partially granted: if the rank F is not an exception, then it is proved that there is a finite number of decompositions with probability 1. These cases are shown in a framebox in Tables 1 and2. In all other cases there are infinitely many decompositions. In the symmetric case, stronger results have been obtained recently and prove that the decomposition is essentially unique with probability 1 if the dimension does not exceed the order [START_REF] Mella | Singularities of linear systems and the Waring problem[END_REF]. In the case of non symmetric tensors, such proofs are still lacking.

So why do we go to the approximation problem? In practice, imposing a smaller rank is a means to obtain an essentially unique decomposition (cf. Section 2.4). In fact, it has been proved that symmetric tensors of any order d ≥ 3 and of dimension 2 or 3 always have an essentially unique decomposition with probability 1 if their rank is strictly smaller than the generic rank [START_REF] Chiantini | On the concept of k-sectant order of a variety[END_REF]. Our conjecture is that this also holds true for non symmetric tensors: all tensors with a sub-generic rank have a unique decomposition with probability 1, up to scale and permutation indeterminacies.

Another result developed in the context of Arithmetic Complexity has been obtained by Kruskal in 1977 [START_REF] Kruskal | Three-way arrays: Rank and uniqueness of trilinear decompositions[END_REF]. The proof has been later simplified [START_REF] Stegeman | On Kruskal's uniqueness condition for the Candecomp/Parafac decomposition[END_REF] [START_REF] Landsberg | Kruskal's theorem[END_REF], and extended to tensors of order d higher than 3 [START_REF] Sidiropoulos | On the uniqueness of multilinear decomposition of N-way arrays[END_REF]. The theorem states that essential uniqueness of the CP decomposition (equation (2) for d = 3) is ensured if the sufficient condition below is satisfied:

2 F + d -1 ≤ d ℓ=1 min(I ℓ , F )
where I ℓ denotes the ℓth dimension. It is already known that this condition is not necessary. Our conjecture goes in the same direction, and our tables clearly show that Kruskal's condition on F above is much more restrictive than the condition that F is strictly smaller than generic. Another interesting result has been obtained in [START_REF] Delathauwer | A link between canonical decomposition in multilinear algebra and simultaneous matrix diagonalization[END_REF] for tensors having one large dimension. The condition of uniqueness is then less restrictive than the above.

Dimension reduction

When fitting a parametric model, there is always some loss of information; the fewer parameters, the larger the loss. Detecting the optimal number of free parameters to assume in the model, for a given data length, is a standard problem in parameter estimation [START_REF] Rissanen | Stochastic Complexity in Statistical Inquiry[END_REF]. If it is taken too small, there is a bias in the model; but if it is too large, parameters suffer from a large estimation variance. At the limit, if the variance of some parameters are infinite, it means that uniqueness is lost for the corresponding parameters. In our model [START_REF] Abo | Induction for secant varieties of Segre varieties[END_REF], the number of parameters is controlled by the rank F and the dimensions. If the rank F is taken too small, the factors are unable to model the given tensor with enough accuracy; if it is too large, factors are poory estimated, or can even become non unique even if the permutation-scale ambiguity is fixed.

It is clear that assuming that the rank F is generic will lead either to large errors or to non unique decompositions. Reducing the rank is a means to eliminate part of the estimation variance, noise effects and measurement errors. Restoring the essential uniqueness of the model by rank reduction may be seen as a limitation of the estimation variance, and should not be seen as a side effect.

In this section, we propose to limit the rank to a subgeneric value and to reduce the dimensions. Statistical analyses leading to the optimal choices of dimensions and rank are out of our scope, and only numerical aspects are subsequently considered.

Unfolding matrices and multilinear rank. One way to reduce the dimensions of the problem, and at the same time reduce the rank of T , is to truncate the Singular Value Decomposition (SVD) of unfolding matrices of array T. Let's see how to do that. An order d array T admits d different unfolding matrices, also sometimes called flattening matrices (other forms exist but are not fundamentally different). Take a I × J × K array. For k fixed in the third mode, we have a I × J matrix that can be denoted as T ::k . The collection of these K such matrices can be arranged in a KI × J block matrix:

T KI×J =         T ::1 . . .
T ::k . . .

T ::K        
The two other unfolding matrices, T IJ×K and T JK×I , can be obtained in a similar manner as T KI×J , by circular permutation of the modes; the latter contain blocks of size J × K and K × I, respectively. The rank of the unfolding matrix in mode p may be called the p-mode rank. It is known that tensor rank is bounded below by every mode-p rank [START_REF] Delathauwer | On the best rank-1 and rank-(R1,R2,...RN) approximation of high-order tensors[END_REF]. The d-uple containing all p-mode ranks is sometimes called the multilinear rank.

Kronecker product. Let two matrices A and H, of size I × J and K × L respectively. One defines the Kronecker product A ⊗ H as the IK × JL matrix [START_REF] Rao | Linear Statistical Inference and its Applications[END_REF] [51]:

A ⊗ H def =    a 11 H a 12 H • • • a 21 H a 22 H • • • . . . . . .   
Now denote a j and h ℓ the columns of matrices A and H. If A and H have the same number of columns, on can define the Khatri-Rao product [START_REF] Rao | Linear Statistical Inference and its Applications[END_REF] as the IK × J matrix:

A ⊙ H def = a 1 ⊗ h 1 a 2 ⊗ h 2 • • • .
The Khatri-Rao product is nothing else but the column-wise Kronecker product. Note that the Kronecker product and the tensor product are denoted in a similar manner, which might be confusing. In fact, this usual practice has some reasons: A ⊗ H is the array of coordinates of the tensor product of the two associated linear operators, in some canonical basis.

High-Order SVD (HOSVD). The Tucker3 decomposition has been introduced by Tucker [START_REF] Tucker | Some mathematical notes for three-mode factor analysis[END_REF] and studied in depth by L.DeLathauwer [24] [25]. It is sometimes referred to as the High-Order Singular Value Decomposition (HOSVD). A third order array of size I × J × K can be decomposed into an all-orthogonal core array T c of size R 1 × R 2 × R 3 by orthogonal changes of coordinates, where R p denote the mode-p ranks, earlier defined in this section. One can write:

T = (Q (1) , Q (2) , Q (3) ) • T c ⇔ T c = (Q (1)T , Q (2)T , Q (3)T ) • T (5)
In order to compute the three orthogonal matrices Q (ℓ) , one just computes the SVD of the three unfolding matrices of T. The core tensor T c is then obtained by the rightmost equation in [START_REF] Brambilla | On the Alexander-Hirschowitz theorem[END_REF]. Now the advantage of the SVD is that it allows to compute the best low-rank approximate of a matrix, by just truncating the decomposition (setting singular values to zero, below some chosen threshold); this is the well-known Eckart-Young theorem. So one could think that truncating the HOSVD in every mode would yield a best lower multilinear rank approximate. It turns out that this is not true, but generally good enough in practice.

By truncating the HOSVD to some lower multilinear rank, say (r 1 , r 2 , r 3 ) smaller than (R 1 , R 2 , R 3 ), one can reduce the tensor rank. But this is not under control, since tensor rank may be larger than all mode-ranks: F ≥ r i . In other words, HOSVD truncation may yield a sub-generic tensor rank, but the exact tensor rank will not be known. To conclude, and contrary to what may be found in some papers in the literature, there does not exist a generalization of the Eckart-Young theorem to tensors.

Algorithms

Various iterative numerical algorithms have been investigated in the literature, and aim at minimizing the fitting error (3). One can mention in particular the Conjugate Gradient [START_REF] Paatero | The multilinear engine: A table-driven, least squares program for solving multilinear problems, including the n-way parallel factor analysis model[END_REF], Gauss-Newton and Levenberg-Marquardt [START_REF] Tomasi | A comparison of algorithms for tting the Parafac model[END_REF] algorithms.

Notations. It is convenient to define the vec operator that maps a matrix to a vector: m = vec(M), meaning that m (i-1)J+j = M ij . Also denote by ⊡ the Hadamard entry-wise matrix product between matrices of the same size. With this notation, the gradient of Υ with respect to vec(A) is given by:

g A = [I A ⊗ (C T C ⊡ B T B)]vecA T -[I A ⊗ (C ⊙ B)] T vecT JK×I
and other gradients are deduced by circular permutation:

g B = [I B ⊗ (A T A ⊡ C T C)]vecB T -[I B ⊗ (A ⊙ C)] T vecT KI×J g C = [I C ⊗ (B T B ⊡ A T A)]vecC T -[I C ⊗ (B ⊙ A)] T vecT IJ×K Also define p =   vecA T vecB T vecC T   , and g =   g A g B g C  
Lastly, we have the compact forms below for Jacobian matrices:

J A = I A ⊗ (C ⊙ B) J B = Π 1 [I B ⊗ (A ⊙ C)] J C = Π 2 [I C ⊗ (B ⊙ A)]
where Π 1 and Π 2 are permutation matrices that put the entries in the right order.

The goal is to minimize the Frobenius norm of tensor T -(A, B, C) • I, which may be written in a non unique way as a matrix norm:

Υ(p) def = ||T JK×I -(B ⊙ C)A T || 2 (6)

Gradient descent

The gradient descent is the simplest algorithm. The iteration is given by: p(k +1) = p(k)µ(k) g(k), where µ(k) is a stepsize that is varied as convergence progresses.

Even with a good strategy of variation of µ(k), this algorithm has often very poor convergence properties. But they can be dramatically improved with the help of ELS, as will be demonstrated in Section 5.4.

Levenberg-Marquardt (LM)

Taking into account the second derivative (Hessian) allows faster local convergence. An approximation using only first order derivatives is given by the iteration below, known as Levenberg-Marquardt:

p(k + 1) = p(k) -[J(k) T J(k) + λ(k) 2 I] -1 g(k)
where J(k) denotes the Jacobian matrix [J A , J B , J C ] at iteration k and λ(k) 2 a positive regularization parameter. There exist several ways of calculating λ(k) 2 , and this has an important influence on convergence. In the sequel, we have adopted the update described in [START_REF] Madsen | Methods for Non-Linear Least Squares Problems[END_REF]. Updates of p and λ 2 are controlled by the gain ratio γ:

γ def = Υ(k) -Υ(k + 1) • ( Υ(0) -Υ(∆p(k))) -1
where Υ(∆p(k)) is a second order approximation of Υ(p(k) + ∆p(k)). More precisely, each new iteration of the algorithm follows this simple scheme:

1. Find the new direction ∆p(k) = -[J(k) H J(k) + λ(k) 2 I] -1 g(k), 2

Alternating Least Squares (ALS)

The principle of the ALS algorithm is quite simple. We recall here the various versions that exist, since there are indeed several. For fixed B and C, this is a quadratic form in A, so that there is a closed form solution for A. Similarly for B and C. Consquently, we can write:

A T = (B ⊙ C) † T JK×I def = f A (B, C) B T = (C ⊙ A) † T KI×J def = f B (C, A) (7)
C T = (A ⊙ B) † T IJ×K def = f C (A, B)
where M † denotes the pseudo-inverse of M. When M is full column rank, we have M † = (M T M) -1 M T , which is the Least Squares (LS) solution of the overdetermined linear system. But it may happen that M is not full column rank. In that case, not only there are more equations than unknowns, so that the system is impossible, but some equations are linked to each other so that there are also infinitely many solutions for some unknowns. The pseudo-inverse then gives the solution of minimal norm for those variables (i.e. null component in the kernel).

In some cases, this may not be what we want: we may want to obtain full rank estimates of matrices (A, B, C).

Another well-known approach to singular systems is that of regularization. It consists of adding a small diagonal matrix D, which allows to restore invertibility. Matrix D is often taken to be proportional to the identity matrix D = ηI, η being a small real number. This corresponds to define another inverse of the form M † = (M T M + ηD) -1 M T . With this definition, we would have three other relations, with obvious notation:

A T = fA (B, C; D), B T = fB (C, A; D), C T = fC (A, B; D) (8)
Thus from ( 7) one can compute the updates below

ALS1. B n+1 = f B (C n , A n ), C n+1 = f C (A n , B n ), A n+1 = f A (B n , C n ) ALS2. B n+1 = f B (C n , A n ), C n+1 = f C (A n , B n+1 ), A n+1 = f A (B n+1 , C n+1 ), ALS3. B n+1 = fB (C n , A n ; ηI), C n+1 = fC (A n , B n+1 ; ηI), A n+1 = fA (B n+1 , C n+1 ; ηI)
Algorithm ALS1 is never utilized in practice except to compute a direction of search at a given point, and it is preferred to replace the values of the loading matrices by the most recently computed, as shown in ALS2. Algorithm ALS3 is the regularized version of ALS2. If the tensor T to be decomposed is symmetric, there is only one matrix A to find. At least two implementations can be thought of (more exist for higher order tensors):

ALS4. Soft constrained:

A n+1 = f A (A n , A n-1 )
ALS5. Hard constrained:

A n+1 = f A (A n , A n )
There does not exist any proof of convergence of the ALS algorithms, which may not converge, or converge to a local minimum. Their extensive use is thus rather unexplainable.

Enhanced Line Search (ELS)

The ELS enhancement is applicable to any iterative algorithm, provided the optimization criterion is a polynomial or a rational function. Let ∆p(k) be the direction obtained by an iterative algorithm, that is, p(k + 1) = p(k) + µ ∆p(k). ELS ignores the value of µ that the iterative algorithm has computed, and searches for the best stepsize µ that corresponds to the global minimum of (6):

Υ(p(k) + µ ∆p(k)).
Stationary values of µ are given by the roots of a polynomial in µ, since Υ is itself a polynomial. It is then easy to find the global minimum once all the roots have been computed, by just selecting the root yielding the smallest value of Υ [START_REF] Rajih | Enhanced line search: A novel method to accelerate Parafac[END_REF]. Note that the polynomial is of degree 5 for third order tensors. The calculation of the ELS stepsize can be executed every every p iteration, p ≥ 1, in order to spare some computations (cf. Section 4).

Numerical complexity

The goal of this section is to give an idea of the cost per iteration of various algorithms. In fact, to compare iterative algorithms only on the basis on the number of iterations they require is not very meaningful, because the costs per iteration are very different.

Nevertheless, only orders of magnitudes of computational complexities are given. In fact, matrices that are the Kronecker product of others are structured. Yet, the complexity of the product between matrices, the solution of linear systems or other matrix operations, can be decreased by taking into account their structure. And such structures are not fully taken advantage of in this section.

Singular Value Decomposition (SVD).

Let the so-called reduced SVD of a m × n matrix M of rank r, m ≥ n ≥ r, be written as:

M = U Σ V T = U r Σ r V r T where U and V are orthogonal matrices, Σ is m×n diagonal, U r is a m×r submatrix of U, U r T U r = I r , Σ r is r×r diagonal and V r is n×r, V r T V r = I r .
The calculation of the diagonal matrix Σ needs 2mn 2 -2n 3 /3 multiplications, if matrices U and V are not explicitly formed, but kept in the form of a product of Householder symmetries and Givens rotations [START_REF] Golub | Matrix computations[END_REF]. Calculating explicitly matrices U r and V r requires additional 5mr 2r 3 /3 and 5nr 2r 3 /3 multiplications, respectively.

For instance if r = n, the total complexity for computing Σ, U n and V is of order 7mn 2 + 11n 3 /3. If m ≫ n, this complexity can be decreased to O(3mn 2 ) by resorting to Chan's algorithm [START_REF] Chan | An improved algorithm for computing the singular value decomposition[END_REF] [START_REF] Comon | Tracking of a few extreme singular values and vectors in signal processing[END_REF]. Even if in the present context this improvement could be included, it has not yet been implemented.

High-Order SVD (HOSVD).

Computing the HOSVD of a I × J × K array requires three SVD of dimensions I × JK, J × KI and K × IJ respectively, where only the left singular matrix is required explicitly. Assume that we want to compute in mode i the reduced SVD of rank R i . Then the overall computational complexity for the three SVDs is

2IJK(I +J +K)+5(R 2 1 JK +IR 2 2 K +IJR 2 3 )-2(I 3 +J 3 +K 3 )/3-(R 3 1 +R 3 2 +R 3 3 )/3
In addition, if the core tensor needs to be computed, singular matrices need to be contracted on the original tensor. The complexity depends on the order in which the contractions are executed. If we assume I ≥ J ≥ K > R i , then it is better to contract the third mode first and the first mode last. This computation

requires IR 1 JK + R 1 JR 2 K + R 1 R 2 KR 3 additional multiplications. If R i ≈ F ≪ min(I, J, K
) one can simplify the overall complexity as:

2IJK(I + J + K) -2(I 3 + J 3 + K 3 )/3 + IJKF
Alternating Least Square (ALS). Let's look at the complexity of one iteration for the first mode, and consider the linear system below, to be solved for A in the LS sense:

M A T = T
where the dimensions of M, A and T are m × r, r × q and m × q, respectively. First, matrix M = B ⊙ C needs to be calculated, which requires one Khatri-Rao product, that is, JKF multiplications. Then the SVD of matrix M of size F × JK needs 7JKF 2 + 11 F 3 /3 multiplications, according to the previous paragraph. Last, the product

V F Σ -1 F U F T T represents F JKI + F I + F 2 I multiplications.
The cumulated total of one iteration of ALS2 for the three modes is hence

(JK + KI + IJ)(7F 2 + F ) + 3F IJK + (I + J + K)(F 2 + F ) + 11 F 3 ,
where the two last terms are often negligible. The complexity of ALS3 is not detailed here, but could be calculated following the same lines as LM in the next section.

Levenberg-Marquardt (LM). One iteration consists of calculating ∆p

def = [J T J + λ 2 I] -1 g where J is of size m × n, m def = IJK, n def = (I + J + K)F .
Of course, it is not suitable to compute the product J T J, nor the inverse matrix. Instead, one shall work with the (m + n) × n matrix S def = [J, λI] T , and consider the linear system S T S ∆p = g.

As in ALS, the first step is to compute the Khatri-Rao products. Then each of the three gradients requires (I + J + K)F 2 + IJKF multiplications. Next, the orthogonalization (QR factorization) of matrix S requires mn 2 multiplications, since making the orthogonal matrix explicit is not necessary. Finally, g is obtained by solving two triangular systems, each requiring an order of n 2 /2 multiplications. Thus the overall complexity of one iteration of the algorithm is dominated by the orthogonalization step, which costs:

IJK(I + J + K) 2 F 2 multiplications.
Enhanced Line Search (ELS). Given a direction of search, the complexity of one ELS iteration is dominated by the calculation of the 6 coefficients of the polynomial in µ to be rooted. This complexity may be shown to be of order (8F + 10)IJK [START_REF] Rajih | Enhanced line search : A novel method to accelerate PARAFAC[END_REF]. But the constant coefficient is not needed here so that we may assume that (8F + 9)IJK multiplications are required.

It is interesting to compare the complexity of ELS with those of ALS or LM. This gives simple rules of thumb. If the rank is smaller than the bound below, then the additional complexity of ELS is negligible compared to ALS2 [START_REF] Rajih | Enhanced line search : A novel method to accelerate PARAFAC[END_REF]:

1 I + 1 J + 1 K F ≫ 1 (9)
Similarly, if 8 ≪ (I + J + K)F , ELS is negligible compared to LM, which in practice always happens. These rules are extremely simplified, but they still show that the larger the rank F , the more negligible ELS; they also have the advantage to give quite reliable orders of magnitude when dimensions are large enough (i.e. at least 10).

5 Computer results

Performance measure

Performance evaluation is carried out in terms of error between the calculated and actual factor matrices. The problem would be easy if there were not an indetermination up to scale and permutation. In fact, it is well known that

(A, B, C) • I = (AΛ A P, BΛ B P, CΛ C P) • I For any F × F invertible diagonal matrices Λ A , Λ B , Λ C satisfying Λ A Λ B Λ C = I
and for any F × F permutation matrix P. It is thus desirable that the performance measure be invariant under the action of these four matrices. First, it is convenient to define a distance between two column vectors u and v, which is invariant to scale. This is given by the expression below:

δ(u, v) = ||u - v T u v T v v|| (10)
Now we describe two ways of computing a scale-permutation invariant distance between two matrices M and M.

Optimal algorithm. For every F × F permutation P, one calculates

∆(P) = K k=1 δ(u(k), v(k))
where u(k) and v(k) denote the kth column of matrices

  A B C   P et M def =   Â B Ĉ  
respectively. The error retained is the minimal distance obtained over all F ! possible permutations. Of course, this optimal algorithm is usable only for moderate values of F , say smaller than 9. For larger ranks, one must resort to a suboptimal algorithm, as the one described below.

Greedy algorithm. The goal is to avoid to describe all permutations. For doing this, one selects an hopefully relevant subset of all possible permutations. One first selects the column m(j 1 ) of M having the largest norm and one calculates its distance δ(•) with the closest column m(k 1 ) of M:

∆ 1 = min k1 δ(m(k 1 ), m(j 1 ))
Then one deletes from matrices M and M the columns k 1 and j 1 previously selected, and one keeps going with the remaining columns. At iteration n, one has to compute Fn + 1 distances δ(m(k n ), m(j n )). After F iterations, one eventually obtains a suboptimal distance ∆, defined either as

∆ = K k=1 ∆ k or as ∆ 2 = K k=1 ∆ 2 k .
If M is L × F , this greedy calculation of the error ∆ requires an order of 3F 2 L/2 multiplications.

-0.4 -0. Paatero gave in his paper [START_REF] Paatero | Construction and analysis of degenerate Parafac models[END_REF] a nice example of a 2 × 2 × 2 tensor, fully parameterizable, hence permitting reproducible experiments. We use here the same parameterization. Its form is as follows:

T = 0 1 e 0 1 d 0 h
The discriminant defined in Section 2.2 is in the present case ∆ = 4 h + d 2 e. As in [START_REF] Paatero | Construction and analysis of degenerate Parafac models[END_REF], we represent the decompositions computed in the hyperplane e = 30. The tensor to be decomposed are chosen to have a discriminant of ∆ > 0, so that they have a rank 2, but they lie close to the variety ∆ = 0, represented by a parabola in the subsequent figures.

Case 1. First consider the tensor defined by (e, d, h) = (30, 0.26, -0.34), and start with the initial value (e, d, h) = (30, -0.3, -0.12). In that case, neither the Gradient nor even LM converge towards the solution, and are stuck on the left side of the parabola. The ALS algorithm does not seem to converge, or perhaps does after a very large number of iterations. With the ELS enhancement executed every 3 iterations, ALS2 converges to a neighborhood of the solution (error of 10 -6 ) within only 60 iterations, as shown in Figure 1. Then, the algorithm is slowed down for a long time before converging eventually to the solution; the beginning of this plateau may be seen in the right graph of Figure 1. Note that it does not correspond to a plateau of the objective function. In this particular example, the gain brought by ELS is very significant at the beginning. On the other hand, no gain is brought at the end of the trajectory towards the goal, where ELS cannot improve ALS2 because the directions given by ALS2 are repeatedly bad. the most efficient and converges within 100 iterations, as shown in Figure 3. The fastest ALS2 algorithm is this time the one using the ELS enhancement at every iteration, and converges within 400 iterations (cf. Figure 2).

In the second case, ALS2 with the ELS enhancement could avoid being trapped in a kind of "cycling" about the goal, whereas it could not do it in the first example. This is a limitation of ELS: if the directions provided by the algorithm are bad, ELS cannot improve much on the convergence. The same has been be observed with case 1: ELS does not improve LM since directions chosen by LM are bad (LM is not represented for case 1). The other observation one can make is that executing ELS at every iteration is not always better than at every p iterations, p > 1: this depends on the trajectory, and there is unfortunately no general rule. This will be also confirmed by the next computer experiments.

Random 30 × 30 × 30 tensor of rank 20

In the previous section, the computational complexity per iteration was not an issue, because it is so small. Now we shall turn our attention to this issue, with a particular focus on the effect of dimension reduction by HOSVD truncation. As in the 2 × 2 × 2 case, the initial value has an influence on convergence; one initial point is chosen randomly, the same for all algorithms to be compared. The tensor considered here is generated by three 30 × 20 factor matrices whose entries have been independently randomly drawn according to a zero-mean unit variance Gaussian distribution. In the sequel, we shall denote such a random drawing as A ∼ N 30×20 (0, 1), in short. The correlation coefficient (absolute value of the cosine) between any pair of columns of these matrices was always smaller than 0.8, which is a way to check out that the goal to reach does not lie itself in a strong bottleneck.

The ALS algorithm is run either without ELS, or with ELS enhancement at every iteration. Figure 4 shows the clear improvement brought by ELS without dimension reduction (the computational complexity is 3 times smaller). If the dimensions are reduced to 20 by HOSVD truncation, we also observe a similar improvement: ALS runs better, but the additional complexity of ELS also becomes smaller (the rank is of same order as dimensions). Also notice by comparing top and bottom of Figure 4 that curves 2 and 3 do not appear in the same order: the dimension reduction gives the advantage to ALS in dimension 20 without ELS, compared to ALS with ELS in dimension 30. The advantage was reversed in the top part of the figure .  A reduction of the dimensions to a value smaller than 20, while keeping the rank equal to 20, is theoretically possible. However, we know that truncating the HOSVD is not the optimal way of finding the best [START_REF] Cox | Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra[END_REF][START_REF] Cox | Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra[END_REF][START_REF] Cox | Ideals, Varieties, and Algorithms: An Introduction to Computational Algebraic Geometry and Commutative Algebra[END_REF] multilinear rank approximation. In practice, we have observed that truncating the HOSVD yields reconstruction errors if the chosen dimension is smaller than the rank (e.g. 19 in the present case or smaller).

Random 5 × 5 × 5 tensor of rank 4

The goal here is to show the interest of ELS, not only for ALS, but also for the Gradient algorithm. In fact in simple cases, the Gradient algorithm may be attractive, if it is assisted by the ELS enhancement. 1: Gradient + ELS with period 2, 2: Gradient + ELS with period 5, 3: ALS + ELS with period 1 (every iteration), 4: ALS + ELS with period 5, 5: ALS + ELS with period 7, 6: ALS.

We consider the easy case of a 5 × 5 × 5 tensor of rank 4. The latter has been built from three 5 × 4 Gaussian matrices drawn randomly as before, each N 5×4 (0, 1). Results are reported in figure 5, and show that ELS performs similarly if executed every 2 iterations, or every 5. Results of the Gradient without ELS enhancement are not shown, but are as poor as ALS. Convergence of the Gradient with ELS is reached earlier than ALS with ELS, but this is not always the case in general. The Gradient with ELS appears to be especially attractive in this example if computational complexity is considered.

5.5 Double bottleneck in a 30 × 30 × 30 tensor Now more difficult cases are addressed, and we consider a tensor with two double bottlenecks (with R. Harshman terminology). The tensor is built from three Gaussian matrices A, B and C of size 30 × 4, and we have created two bottlenecks in the two first modes as follows. The first column a(1) of A is drawn randomly as N 30 (0, 1). Then the second column is set to be a(2) = a(1) + 0.5 v(2), where v(2) ∼ N 30 (0, 1). The third and fourth columns of A are generated in the same manner, that is, a(4) = a(3)+0.5 v(4), where a(3) ∼ N 30 (0, 1) and v(4) ∼ N 30 (0, 1). Matrix B is independently generated exactly in the same manner. Vectors belonging to a bottleneck generated this way have a correlation coefficient of 0.89. On the other hand, matrix C is independently N 30×4 (0, 1), and its columns are not close to dependent.

Results reported in Figure 6 show that even when the rank is moderate compared to dimensions, ELS brings improvements to ALS, despite the fact that its additional complexity is not negligible according to Equation [START_REF] Carroll | Analysis of individual differences in multidimensional scaling via n-way generalization of Eckart-Young decomposition[END_REF]. Surprisingly, a dimension reduction down to the rank value does not change this conclusion. As previously, ELS with period 1 cannot be claimed to be always better than ELS with period p > 1.

ALS with ELS enhancement at every 5 iterations (solid line in Figure 7) benefits from an important computational complexity improvement thanks to the dimension reduction down to 4: its convergences needs 6.10 6 multiplications instead of 8.10 8 . The improvement is even more impressive for ALS with ELS enhancement at every iteration (dashed line).

5.6 Swamp in a 30 × 30 × 30 tensor of rank 4

In this last example, we have generated a tensor with two triple bottlenecks, that is, a tensor that lies in a swamp. Matrices A and B are generated as in Section 5.5, but matrix C is now also generated the same way. This tensor is difficult to decompose, and ALS does not converge after 5000 iterations. With the help of ELS every 5 iterations, the algorithm converges with about a thousand of iterations. This was unexpected since ELS is not always able to escape from swamps. The reduction of dimensions to 4 allows to decrease the computational complexity by a factor of 50, even if the number of iterations has slightly increased.

Our paper was three-fold. Firstly we have surveyed theoretical results, and pointed out which are still conjectures. Secondly we have described several numerical algorithms and given their approximate complexities for large dimensions. Thirdly, we have compared these algorithms in situations that Richard Harshman had classified, and in particular multiple bottlenecks and swamps. Several conclusions can be drawn.

First, one can claim that dimension reduction by HOSVD truncation is always improving the results, both in terms of speed of convergence and computational complexity, provided the new dimensions are not smaller than the rank of the tensor to be decomposed.

Second, ELS has been shown to be useful to improve ALS2 or the Gradient descent, either as an help to escape from a local minimum or to speed up convergence. ELS has been often efficient in simple or multiple bottlenecks, and even in swamps. Its help ranges from excellent to moderate. On the other hand, ELS has never been able to help LM when it was stuck.

Third, ALS2 can still progress very slowly in some cases where the objective function is not flat at all. This is a particularity of ALS, and we suspect that it is "cycling". In such a case, ELS does not help because the directions defined by ALS2 are not good. One possibility would be to call for ALS1 to give other directions; this is currently being investigated.

Last, no solution has been proposed in cases of CP-degeneracies. In the later case, the problem is ill-posed, and the solution should be searched for by changing the modeling, i.e., either by increasing the rank or by imposing additional constraints that will avoid factor cancellations. Future works will be focussed on the two latter conclusions.

3 , 1 -

 31 . Compute p(k + 1) and deduce Υ(k + 1) 3. Compute Υ(∆p(k)) 4. Compute γ. 5. If γ ≥ 0, then p(k + 1) is accepted, λ(k + 1) 2 = λ(k) 2 * max 1 (2γ -1)3 and ν = 2. Otherwise, p(k + 1) is rejected, λ(k + 1) 2 = νλ(k) 2 and ν ← 2ν.

Figure 1 :

 1 Figure 1: First case treated by ALS with ELS at every 3 iterations. Left: symbols '+', 'o' and ' * ' denote respectively the initial value, the successive iterations, and the global minimum, in the (d, h) plane. Right: value of the error as a function of iterations.

5. 2 A 2 × 2 × 2

 2222 tensor of rank 2

Figure 2 :

 2 Figure 2: Second case treated by ALS with ELS at every iteration. Left: symbols '+', 'o' and ' * ' denote respectively the initial value, the successive iterations, and the global minimum, in the (d, h) plane. Right: value of the error as a function of iterations.

Case 2 .

 2 Consider now a second (easier) case given by (e, d, h) = (30, 0.26, -0.29), with the initial value (e, d, h) = (30, -0.26, -0.19). This time, the LM algorithm is

Figure 3 :

 3 Figure 3: Second case treated by LM. Left: symbols '+', 'o' and ' * ' denote respectively the initial value, the successive iterations, and the global minimum, in the (d, h) plane. Right: value of the error as a function of iterations.

Figure 4 :

 4 Figure 4: Reconstruction error of a 30 × 30 × 30 random tensor of rank 20 as a function of the number of iterations (top) or the number of multiplications (bottom), and for various algorithms. 1: ALS + HOSVD + ELS, 2: ALS + ELS, 3: ALS + HOSVD, 4: ALS.

Figure 5 :

 5 Figure 5: Reconstruction error of a 5 × 5 × 5 tensor of rank 4 as a function of the number of iterations (top) or the number of multiplications (bottom), and for various algorithms.1: Gradient + ELS with period 2, 2: Gradient + ELS with period 5, 3: ALS + ELS with period 1 (every iteration), 4: ALS + ELS with period 5, 5: ALS + ELS with period 7, 6: ALS.

Figure 7 :

 7 Figure 7: Reconstruction error of a 30 × 30 × 30 of rank 4 with 2 triple bottlenecks as a function of the number of iterations (top) or the number of multiplications (bottom), and for various algorithms. Solid: ALS + ELS with period 5, dotted: ALS. Top: without dimension reduction, bottom: after reducing dimensions down to to 4.

Table 1 :

 1 . Values reported in bold correspond to smallest typical ranks that are now known only via this numerical computation. Smallest typical rank for some third order unconstrained arrays. Frame: finite number of solutions. Underlined: exceptions to the ceil rule.

	K			2				3		4
	J		2	3	4	5	3	4	5	4	5
		2	2	3	4	4	3	4	5	4	5
		3	3	3	4	5	5	5	5	6	6
		4	4	4	4	5	5	6	6	7	8
		5	4	5	5	5	5	6	8	8	9
		6	4	6	6	6	6	7	8	8	10
	I	7	4	6	7	7	7	7	9	9	10
		8	4	6	8	8	8	8	9	10	11
		9	4	6	8	9	9	9	9	10	12
	10	4	6	8	10	9	10 10	10	12
	11	4	6	8	10	9	11 11	11	13
	12	4	6	8	10	9	12 12	12	13
	I		2	3	4	5		6	7	8	9
	d 3	2	5	7	10	14	19	24	30
	4	4	9	20	37	62	97

Table 2 :

 2 Smallest typical rank for unconstrained arrays of order d = 3 and d = 4 and with equal dimensions I. Frame: finite number of solutions. Underlined: exceptions to the ceil rule.

Acknowledgment. This work has been partly supported by contract ANR-06-BLAN-0074 "Decotes".