Tensor Package


Introduction

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 sometimes 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 papers below for further details:

[1] P. Comon, ``Tensors: a brief introduction,''   IEEE Sig. Proc. Magazine, 31(3):44-53, May 2014. Preprint: hal-00923279

[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:
expression of CP decomposition
where one can impose the columns of matrices A, B, and C to have unit norm columns. In compact notation, this is noted as
Compact form of CP
The minimal number of such rank-1 terms that are necessary to reconstruct tensor T exactly is called the rank of T. Suppose 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.
objective to 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:
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, see:

[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

Downloads
Creative Commons License
Tensor package by Pierre Comon is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.
Based on a work at P.Comon's Tensor Package home page.

General-purpose codes

Most of the following codes have been written or rewritten by X.Luciani during his post-doc with P.Comon at I3S in 2009:
Download all above m-files in a single zip archive.

References:

[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
[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

Codes for rank-1 deflation

• These codes have been developed by Alex P. da Silva under supervision of P.Comon:
Download all above m-files in a single zip archive.

References:

  • [3] A. P. da Silva, P. Comon and A. L. F. de Almeida, ``An Iterative Deflation Algorithm for Exact CP Tensor Decomposition'', ICASSP, April 19-24, 2015, Brisbane, Australia.  hal-01071402 DOI: 10.1109/ICASSP.2015.7178714
  • A. P. da Silva, P. Comon and A. L. F. de Almeida, ``Rank-1 Tensor Approximation Methods and Application to Deflation'', arxiv:1508.05273, August 21, 2015.
    Shorter version published in:
    [4] A. Pereira da Silva, P. Comon, and A. Lima Ferrer de Almeida, ``A finite algorithm to compute rank-1 tensor approximations,'' Sig. Proc. Letters, 23(7):959-963, July 2016. DOI: 10.1109/LSP.2016.2570862

  • Codes for real nonnegative tensors

    • Other codes for tensors with nonnegative entries, and decomposed into nonnegative rank-1 terms. These codes have been developed by J-P.Royer under supervision of N.Thirion and P.Comon, and include:
    Download the codes dedicated to tensors with nonNegative entries in a single zip archive.

    Reference:
    [5] 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.

    • When tensors have a large size, we developed special codes able to compress the data (low rank) while at the same time maintain nonnegativity. These codes have been developed by J.E.Cohen and R.Cabral-Farias under the supervision of P.Comon:

    Download the codes dedicated to big nonnegative tensors in a single zip archive

    Reference:
    [6] J. E. Cohen, R. Cabral-Farias and P. Comon, ``Fast Decomposition of Large Nonnegative Tensors'',  IEEE Sig. Proc. Letters, vol.22, no.7, pp.862-866, July 2015. Preprint:  hal-01069069.

    Codes for coupled decompositions

    • When several data sets are available in the form of tensors of different sizes with latent vectors related to each other through (possibly noisy) linear relations, joint flexible CP decomposition can be performed. In the contribution below, A Bayesian (MAP) approch is developed :

    Download the codes to compute CP decompositions with flexible couplings.

    Reference:
    [8] R. Cabral Farias, J. E. Cohen and P. Comon , ``Exploring multimodal data fusion through joint decompositions with flexible couplings'',  IEEE Trans. Signal Processing, vol.64, no.18, pp.4830-4844, Sept. 2016. Preprints:  hal-01158082, arxiv:1505.07717.

    Codes for special problems

    • When factor matrices are structured (such as Toeplitz or Hankel), the CP decomposition can be computed in a noniterative manner.  Computer experiments presented in the ICASSP reference below had been partly written in Maple and run by M.Sorensen and E.Trigaridas. But the codes downloadable below have been written later in Matlab by P.Comon:

    Download the codes to compute the CP decomposition when a factor matrix is structured, in a single zip archive.

    Reference:
    [9] P. Comon, M. Sorensen, and E. Tsigaridas , ``Decomposing tensors with structured matrix factors reduces to rank-1 approximations'',  ICASSP, Dallas, , March 14-19, 2010. Preprint:  hal-00490248.

    • 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:
    [10] 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.

    Codes used in the PhD thesis of Jeremy E. Cohen

    • Personal codes: • Modified toolboxes: Reference:
    [11] J. E. Cohen, ``Environmental multiway data mining '', PhD thesis, Sept. 15, 2016, Univ. Grenoble Alpes. Reprint: tel-01371777

    Codes for computing generic and typical ranks

    • This code provides the rank(s) that a tensor will have with probability one, as a function of order and dimensions. In the complex field, this rank is unique and is called generic. In the real field, a tensor can have several ranks with non-zero probability; these ranks are referred to as typical. The calculation of these ranks is based on the Jacobian of the map.
    A main file calls one of the 3 functions below, depending of the type of tensor: symmetric, square non symmetric, or 3rd order with symmetric matrix slices:

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

    References:
    [12] P. Comon and J. ten Berge, ``Generic and Typical Ranks of Three-Way Arrays,'' Icassp'08, Las Vegas, March 30 - April 4, 2008. Reprint: hal-00327627
    arxiv:0802.2371
    [13] P. Comon, J. M. F. ten Berge, et al., ``Generic and Typical Ranks of Multi-Way Arrays,'' Linear Algebra Appl., vol.430, no.11-12, pp.2997-3007, June 2009. Reprint: hal-00410058

    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 BroThree-mode companyDecotes project, Decoda project.