Tensor Package

We consider the decomposition of a tensor in a sum of rank-1 terms. This decomposition is generally referred to as "Canonical Polyadic" decomposition, and sometinmes as "Parafac". In order to fix this discrepancy, many authors have adopted the "CP" acronym, which stands both for "Candecomp-Parafac" and for "Canonical Polyadic". See the introductory paper below for further details:

[1] P. Comon, X. Luciani, and A. L. F. de Almeida, ``Tensor Decompositions, Alternating Least Squares and other Tales,''  Jour. Chemometrics, 23:393--405, August 2009. Preprint: hal-00410057 or local preprint

In order to fix the ideas, let's take the case of 3rd order tensors. Once the three bases are chosen in each of the three vector spaces, the tensor is defined by a 3-way array. Its CP decomposition takes the form:
where one can impose the columns of matrices A, B, and C to have unit norm columns. In compact notation, this is denoted as
The minimal number of such rank-1 terms that are necessary to reconstruct tensor T exactly is called the rank of T. Supose we are given the array T, and we wish to compute the CP decomposition. If the rank is smaller than the so-called expected dimension [1], then the CP is generally unique. However, noise may be present in measurements, so that the array may contain undesired contributions, which we may prefer to remove. In addition, noise may increase the rank beyond the latter bound, so that uniqueness may not be ensured. Hence, before computing the CP, it is often useful to reduce the rank.

• Rank reduction. Removing noise is standard in many applications, and can be performed by truncating the multilinear rank. This rank is a triplet of numbers, which correspond to the matrix rank of three possible flattening matrices (obtained by storing the 3rd order array values in matrix form); see [1] for complementary explanations and bibliographical references. The reduction of the multilinear rank implies (indirectly) tensor rank reduction. Denote G the new tensor obtained from T by multilinear rank reduction.
  
• Optimization. The next step consists of minimizing , i.e. a norm of the difference between actual tensor and model, for increasing values of r until the resulting difference is null. If the Frobenius norm is chosen, we minimize:

The difficulty in minimizing this objective function lies in the fact  it is multimodal (there are local minima) and that there are many variables. Several algorithms have been developed to minimize this criterion:

• Alternating Least Squares (ALS). It is the least attractive as far as performances are concerned. But it is the easiest to program.
• Gradient descent. It can be implemented with fixed or variable step size. In practice, computing the locally optimal step size is unuseful and unnecessarily costly.
• Quasi-Newton. The BFGS implementation is quite efficient.
• Levenberg-Marquardt. The implementation presented in [1] performs quite well.
• Conjugate gradient. This is the implementation of Fletcher-Reeves, with regular reinitializations.
  

Each of these (standard) algorithms defines a direction of search, and a stepsize. In order to permit the algorithm to escape from local minima at a reasonable computational price, the Enhanced Line Seach (ELS) function has been developed; it computes the global minimum of the objective in any given direction. Because the objective function is polynomial, efficient routines are available and computationally little demanding. The largest part of the load consists of computing the polynomial coefficients. All this is done by the els.m function. For further details on the ELS algorithm, described in the case of ALS, see (for the use of ELS in other algorithms, see [4]):

[2] M. Rajih, P. Comon, R. Harshman, ``Enhanced Line Search: A Novel Method to Accelerate Parafac,'' SIAM Journal on Matrix Analysis Appl., 30(3):1148--1171, September 2008. Preprint: hal-00327595 or local preprint

---

• als.m   (Alternating  Least Squares)
  
• grad.m   (Gradient descent)
  
• bfgs.m   (Quasi Newton)
  
• lm.m   (Levenberg-Marquardt)
  
• cg.m   (Conjugate gradient)
  
• els.m     (Enhanced Line Search, for any of the previous iterations)
• main.m  (a main code, showing examples of how to call the previous functions)
  

Download all above m-files in a single zip archive.

• Other codes for tensors with nonnegative entries, and decomposed into positive rank-1 terms. These codes have been developed by J-P.Royer under supervision of N.Thirion and P.Comon, and include:

• gradP.m  (gradient descent)
  
• bfgsP.m  (BFGS implementation of quasi-Newton)
  
• cgP.m  (conjugate gradient)

Download the codes dedicated to tensors with nonNegative entries in a single zip archive.

Reference:
[3] J.P. Royer, N. Thirion and P. Comon, ``Computing the polyadic decomposition of nonnegative third order tensors'',  Signal Processing, Elsevier, vol.91, no.9, pp. 2159-2171, Sept. 2011. Preprint:  hal-00618729.

• Algorithms based on the joint characteristic function aiming at identifying blindly a linear system with more inputs than outputs: the CAF_toolbox, developed by X.Luciani, A.deAlmeida and P.Comon.

Reference:
[4] X. Luciani, A. L. F. de Almeida and P. Comon, ``Blind Identification of Underdetermined Mixtures Based on the Characteristic Function: The Complex Case'', IEEE Trans. Signal Processing, vol.59, no.2, pp.540-553, February 2011. Preprint: hal-00537838.

• Note: all these codes may be used freely, provided the above authors are acknowledged. In addition, any modification performed in the source codes should also be accurately specified by appending comments at the beginning of the source files.

---

See also web sites of:  Rasmus Bro,  Three-mode company,  Decotes project,