## Efficient MATLAB computations with sparse and factored tensors (2007)

Venue: | SIAM JOURNAL ON SCIENTIFIC COMPUTING |

Citations: | 47 - 13 self |

### BibTeX

@ARTICLE{Bader07efficientmatlab,

author = {Brett W. Bader and Tamara G. Kolda},

title = {Efficient MATLAB computations with sparse and factored tensors},

journal = {SIAM JOURNAL ON SCIENTIFIC COMPUTING},

year = {2007},

volume = {30},

number = {1},

pages = {205--231}

}

### OpenURL

### Abstract

In this paper, the term tensor refers simply to a multidimensional or $N$-way array, and we consider how specially structured tensors allow for efficient storage and computation. First, we study sparse tensors, which have the property that the vast majority of the elements are zero. We propose storing sparse tensors using coordinate format and describe the computational efficiency of this scheme for various mathematical operations, including those typical to tensor decomposition algorithms. Second, we study factored tensors, which have the property that they can be assembled from more basic components. We consider two specific types: A Tucker tensor can be expressed as the product of a core tensor (which itself may be dense, sparse, or factored) and a matrix along each mode, and a Kruskal tensor can be expressed as the sum of rank-1 tensors. We are interested in the case where the storage of the components is less than the storage of the full tensor, and we demonstrate that many elementary operations can be computed using only the components. All of the efficiencies described in this paper are implemented in the Tensor Toolbox for MATLAB.

### Citations

1518 | Iterative methods for sparse linear systems - Saad - 2003 |

307 |
Analysis of individual differences in multidimensional scaling via an n-way generalization of Eckart-Young decomposition
- Carroll, Chang
- 1970
(Show Context)
Citation Context ...eing used in an increasing variety of fields in science, engineering, and mathematics. Tensor decompositions date back to the late 1960s with work by Tucker [51], Harshman [19], and Carroll and Chang =-=[9]-=-. Recent decades have seen tremendous growth in this area with a focus towards improved algorithms for computing the decompositions [13, 12, 58, 50]. Many innovations in tensor decompositions have bee... |

262 |
Foundations of the PARAFAC procedure: Models and conditions for an explanatory multimodal factor analysis. UCLA Working Papers Phonet
- Harshman
- 1970
(Show Context)
Citation Context ... decompositions, which are being used in an increasing variety of fields in science, engineering, and mathematics. Tensor decompositions date back to the late 1960s with work by Tucker [51], Harshman =-=[19]-=-, and Carroll and Chang [9]. Recent decades have seen tremendous growth in this area with a focus towards improved algorithms for computing the decompositions [13, 12, 58, 50]. Many innovations in ten... |

231 | A multilinear Singular Value Decomposition
- Lathauwer, Moor, et al.
- 2000
(Show Context)
Citation Context ...s with work by Tucker [51], Harshman [19], and Carroll and Chang [9]. Recent decades have seen tremendous growth in this area with a focus towards improved algorithms for computing the decompositions =-=[13, 12, 58, 50]-=-. Many innovations in tensor decompositions have been motivated by applications in chemometrics [3, 31, 8, 44]. More recently, these methods have been applied to signal processing [10, 11], image proc... |

202 |
Some mathematical notes on three-mode factor analysis
- Tucker
- 1966
(Show Context)
Citation Context ...iency of tensor decompositions, which are being used in an increasing variety of fields in science, engineering, and mathematics. Tensor decompositions date back to the late 1960s with work by Tucker =-=[51]-=-, Harshman [19], and Carroll and Chang [9]. Recent decades have seen tremendous growth in this area with a focus towards improved algorithms for computing the decompositions [13, 12, 58, 50]. Many inn... |

170 |
Three-way arrays: rank and uniqueness of trilinear decompositions with application to arithmetic complexity and statistics
- Kruskal
- 1977
(Show Context)
Citation Context ...∈ RIn for n =1,...,N. The scalar multiplier λr is optional and can be absorbed into one of the factors, e.g., v (r) 1 . The rank of X is defined as the minimal R such that X can be exactly reproduced =-=[28]-=-. The right-hand side of (2.5) is a Kruskal tensor, which is discussed in more detail in section 5. The CP decomposition can also be computed via an ALS algorithm; see, e.g., [44, 50]. Here we briefly... |

147 | Multilinear analysis of image ensembles: Tensorfaces - Vasilescu, Terzopoulos - 2002 |

130 | Sparse matrices in MATLAB: Design and Implementation
- Gilbert, Moler, et al.
- 1992
(Show Context)
Citation Context ...support mathematical operations on them beyond conversion to unfactored format. MATLAB cannot store sparse tensors except for sparse matrices which are stored in compressed sparse column (CSC) format =-=[16]-=-. Mathematica, an alternative to MATLAB, also supports multidimensional arrays, and there is a Mathematica package for working with tensors that accompanies the book [41]. In terms of sparse arrays, M... |

111 |
S.K.Mitra, Generalized inverses of matrices and its applications
- Rao
- 1971
(Show Context)
Citation Context ...⎢ A ⊗ B ≡ ⎢ ⎣ a11B a12B ··· a1JB a21B a22B ··· a2JB . . Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. . .. . aI1B aI2B ··· aIJB ⎤ ⎥ ⎦ ∈ RIK×JL . The Khatri–Rao product =-=[35, 40, 8, 44]-=- of matrices A ∈ R I×K and B ∈ R J×K is A ⊙ B ≡ � � IJ×K a1 ⊗ b1 a2 ⊗ b2 ··· aK ⊗ bK ∈ R . The Hadamard (elementwise) product of matrices A and B is denoted by A ∗ B. See, e.g., [44] for properties of... |

93 | Face transfer with multilinear models
- Vlasic, Brand, et al.
(Show Context)
Citation Context ...innovations in tensor decompositions have been motivated by applications in chemometrics [3, 31, 8, 44]. More recently, these methods have been applied to signal processing [10, 11], image processing =-=[52, 54, 57, 53]-=-, data mining [43, 46, 1, 45], scientific computing [6], and elsewhere [26, 36]. Though this work can be applied in a variety of contexts, we concentrate on operations that are common to tensor decomp... |

84 | Orthogonal Tensor Decompositions
- Kolda
- 2001
(Show Context)
Citation Context ...s an N-way tensor, defined elementwise as � a (1) ◦ a (2) ◦···◦a (N)� i1i2...iN = a (1) i1 a(2) i2 ···a(N) iN for 1 ≤ in ≤ In,n∈ N. Sometimes the symbol ⊗ is used rather than the symbol ◦ (see, e.g., =-=[24]-=-). 2.3. Matricization of a tensor. Matricization, sometimes called unfolding or flattening, is the rearrangement of the elements of a tensor into a matrix. Let X ∈ RI1×I2×···×IN be an order-N tensor. ... |

78 | CubeSVD: a novel approach to personalized Web search
- Sun, Zeng, et al.
(Show Context)
Citation Context ...itions have been motivated by applications in chemometrics [3, 31, 8, 44]. More recently, these methods have been applied to signal processing [10, 11], image processing [52, 54, 57, 53], data mining =-=[43, 46, 1, 45]-=-, scientific computing [6], and elsewhere [26, 36]. Though this work can be applied in a variety of contexts, we concentrate on operations that are common to tensor decompositions, such as the Tucker ... |

75 |
Principal component analysis of three-mode data by means of alternating least squares algorithms
- Kroonenberg, Leeuw
- 1980
(Show Context)
Citation Context ... all n =1,...,N.IfJn = rankn(X) for all n, then the approximation is exact and the computation is trivial. More typically, an alternating least squares (ALS) approach is used for the computation; see =-=[27, 47, 13]-=-. The Tucker decomposition is not unique, but measures can be taken to correct this [20, 21, 22, 48]. Observe that the right-hand side of (2.4) is a Tucker tensor, to be discussed in more detail in se... |

71 | Beyond streams and graphs: dynamic tensor analysis - Sun, Tao, et al. |

69 | On the best rank-1 and rank-(r1,r2,. . .,rn) approximation of higher-order tensors
- Lathauwer, Moor, et al.
(Show Context)
Citation Context ...s with work by Tucker [51], Harshman [19], and Carroll and Chang [9]. Recent decades have seen tremendous growth in this area with a focus towards improved algorithms for computing the decompositions =-=[13, 12, 58, 50]-=-. Many innovations in tensor decompositions have been motivated by applications in chemometrics [3, 31, 8, 44]. More recently, these methods have been applied to signal processing [10, 11], image proc... |

67 |
Multi-way Analysis
- Smilde, Bro, et al.
- 2004
(Show Context)
Citation Context ... in this area with a focus towards improved algorithms for computing the decompositions [13, 12, 58, 50]. Many innovations in tensor decompositions have been motivated by applications in chemometrics =-=[3, 31, 8, 44]-=-. More recently, these methods have been applied to signal processing [10, 11], image processing [52, 54, 57, 53], data mining [43, 46, 1, 45], scientific computing [6], and elsewhere [26, 36]. Though... |

62 |
Algorithm 862: MATLAB tensor classes for fast algorithm prototyping
- Bader, Kolda
(Show Context)
Citation Context ...n 2.1 [5]. 1.1. Related work and software. MATLAB (version 2006a) provides dense multidimensional arrays and operations for elementwise and binary operations. Version 1.0 of our MATLAB Tensor Toolbox =-=[4]-=- extends MATLAB’s core capabilities to support operations such as tensor multiplication and matricization. The previous version of the toolbox also included objects for storing Tucker and Kruskal fact... |

51 |
Facial expression decomposition
- Wang, Ahuja
- 2003
(Show Context)
Citation Context ...innovations in tensor decompositions have been motivated by applications in chemometrics [3, 31, 8, 44]. More recently, these methods have been applied to signal processing [10, 11], image processing =-=[52, 54, 57, 53]-=-, data mining [43, 46, 1, 45], scientific computing [6], and elsewhere [26, 36]. Though this work can be applied in a variety of contexts, we concentrate on operations that are common to tensor decomp... |

50 |
2000, The N-way toolbox for MATLAB
- Andersson, Bro
(Show Context)
Citation Context ...the other hand, is generally considered more suited to numerical computations. There are three well known packages for (dense) tensor decompositions. The N-way toolbox for MATLAB by Andersson and Bro =-=[2]-=- provides a suite of efficient functions and alternating least squares algorithms for decomposing dense tensors into a variety of models including Tucker and CANDECOMP/PARAFAC. The Multilinear Engine ... |

50 | Multi-way analysis in the food industry: Models, algorithms and applications
- Bro
- 1998
(Show Context)
Citation Context ... in this area with a focus towards improved algorithms for computing the decompositions [13, 12, 58, 50]. Many innovations in tensor decompositions have been motivated by applications in chemometrics =-=[3, 31, 8, 44]-=-. More recently, these methods have been applied to signal processing [10, 11], image processing [52, 54, 57, 53], data mining [43, 46, 1, 45], scientific computing [6], and elsewhere [26, 36]. Though... |

50 |
Rank-one approximation to high order tensors
- Zhang, Golub
(Show Context)
Citation Context ...s with work by Tucker [51], Harshman [19], and Carroll and Chang [9]. Recent decades have seen tremendous growth in this area with a focus towards improved algorithms for computing the decompositions =-=[13, 12, 58, 50]-=-. Many innovations in tensor decompositions have been motivated by applications in chemometrics [3, 31, 8, 44]. More recently, these methods have been applied to signal processing [10, 11], image proc... |

45 |
Towards a standardized notation and terminology in multiway analysis
- Kiers
(Show Context)
Citation Context ...pes of tensors. Finally, in section 7, we conclude with a discussion on the Tensor Toolbox, our implementation of these concepts in MATLAB. 2. Notation and background. We follow the notation of Kiers =-=[23]-=-, except that tensors are denoted by boldface Euler script letters, e.g., X, rather than using underlined boldface X. Matrices are denoted by boldface capital letters, e.g., A; vectors are denoted by ... |

39 | Multilinear analysis of image ensembles - Vasilescu, Terzopoulos - 2002 |

35 | Algorithms for numerical analysis in high dimensions
- Beylkin, Mohlenkamp
(Show Context)
Citation Context ...ions in chemometrics [3, 31, 8, 44]. More recently, these methods have been applied to signal processing [10, 11], image processing [52, 54, 57, 53], data mining [43, 46, 1, 45], scientific computing =-=[6]-=-, and elsewhere [26, 36]. Though this work can be applied in a variety of contexts, we concentrate on operations that are common to tensor decompositions, such as the Tucker [51] and canonical decompo... |

35 |
The multilinear engine: A table-driven, least squares program for solving multilinear problems, including the n-way parallel factor analysis model
- Paatero
- 1999
(Show Context)
Citation Context ... suite of efficient functions and alternating least squares algorithms for decomposing dense tensors into a variety of models including Tucker and CANDECOMP/PARAFAC. The Multilinear Engine by Paatero =-=[38]-=- is a FORTRAN code based on on the conjugate gradient algorithm that also computes a variety of multilinear models. The commerical PLS Toolbox [55] provides a number of multidimensional models with an... |

35 | Multi-way analysis: application in chemical sciences - Smilde, Bro, et al. - 2004 |

32 |
Rank, decomposition, and uniqueness for 3-way and N-way arrays
- Kruskal
- 1989
(Show Context)
Citation Context ...R and U =[u 1 ··· u (n) R ] ∈ RIn×R . This is the format that results from a PARAFAC decomposition [19, 9], and we refer to it as a Kruskal tensor due to the work of Kruskal on tensors of this format =-=[28, 29]-=-. We use the shorthand notation from [25]: (5.1) X = �λ ; U (1) ,...,U (N) �. In some cases, the weights λr are not explicit, and we write X = �U (1) ,...,U (N) �. Other notation can be used. For inst... |

31 | Multilinear operators for higherorder decompositions. United States
- Kolda
- 2006
(Show Context)
Citation Context ...ypes of factored tensors that correspond to the Tucker [51] and CANDECOMP/PARAFAC [9, 19] models. The Tucker format stores a tensor as the product of a core tensor and a factor matrix along each mode =-=[25]-=-. For example, if X is a third-order tensor that is stored as the product of a core tensor G of size R × S × T with corresponding factor matrices, then we express it as X = �G ; A, B, C�, which means ... |

30 |
A comparison of algorithms for fitting the PARAFAC model
- Tomasi, Bro
(Show Context)
Citation Context |

30 | Interactive tensor field design and visualization on surfaces - ZHANG, HAYS, et al. - 2007 |

29 | Matlab tensor classes for fast algorithm prototyping - Bader, Kolda - 2004 |

26 |
Some design features of a sparse matrix code
- DUFF, REID
- 1979
(Show Context)
Citation Context ...end to higher dimensions. The simplest storage format is coordinate format, which stores each nonzero along with its row and column index in three separate one-dimensional arrays, which Duff and Reid =-=[14]-=- called “parallel arrays.” For a matrix A of size I × J with nnz(A) nonzeros, the total storage is 3 ·nnz(A), and the indices are not necessarily presorted. More common are CSR and CSC formats, which ... |

26 | Singular values and eigenvalues of tensors: a variational approach
- Lim
(Show Context)
Citation Context ...lts from a Tucker decomposition [51] and is therefore termed a Tucker tensor. We use the shorthand notation �G ; U (1) , U (2) ,...,U (N) � from [25], but other notation can be used. For example, Lim =-=[32]-=- proposes that the covariant aspect of the multiplication be made explicit by expressing (4.1) as � U (1) , U (2) ,...,U (N)� · G. Copyright © by SIAM. Unauthorized reproduction of this article is pro... |

19 |
Tensor decompositions, state of the art and applications
- COMON
- 2002
(Show Context)
Citation Context ...ons [13, 12, 58, 50]. Many innovations in tensor decompositions have been motivated by applications in chemometrics [3, 31, 8, 44]. More recently, these methods have been applied to signal processing =-=[10, 11]-=-, image processing [52, 54, 57, 53], data mining [43, 46, 1, 45], scientific computing [6], and elsewhere [26, 36]. Though this work can be applied in a variety of contexts, we concentrate on operatio... |

18 |
A simple comprehensive model for the analysis of covariance structures
- McDonald
- 1978
(Show Context)
Citation Context ...⎢ A ⊗ B ≡ ⎢ ⎣ a11B a12B ··· a1JB a21B a22B ··· a2JB . . Copyright © by SIAM. Unauthorized reproduction of this article is prohibited. . .. . aI1B aI2B ··· aIJB ⎤ ⎥ ⎦ ∈ RIK×JL . The Khatri–Rao product =-=[35, 40, 8, 44]-=- of matrices A ∈ R I×K and B ∈ R J×K is A ⊙ B ≡ � � IJ×K a1 ⊗ b1 a2 ⊗ b2 ··· aK ⊗ bK ∈ R . The Hadamard (elementwise) product of matrices A and B is denoted by A ∗ B. See, e.g., [44] for properties of... |

17 |
Simplicity of core arrays in three-way principal component analysis and the typical rank of p×q×2 arrays. Linear Algebra and its Applications
- Berge, Kiers
- 1999
(Show Context)
Citation Context ...s trivial. More typically, an alternating least squares (ALS) approach is used for the computation; see [27, 47, 13]. The Tucker decomposition is not unique, but measures can be taken to correct this =-=[20, 21, 22, 48]-=-. Observe that the right-hand side of (2.4) is a Tucker tensor, to be discussed in more detail in section 4. The CANDECOMP/PARAFAC decomposition was simultaneously developed as the canonical decomposi... |

14 |
Multilinear models: applications in spectroscopy
- Leurgans, Ross
- 1992
(Show Context)
Citation Context ... in this area with a focus towards improved algorithms for computing the decompositions [13, 12, 58, 50]. Many innovations in tensor decompositions have been motivated by applications in chemometrics =-=[3, 31, 8, 44]-=-. More recently, these methods have been applied to signal processing [10, 11], image processing [52, 54, 57, 53], data mining [43, 46, 1, 45], scientific computing [6], and elsewhere [26, 36]. Though... |

11 |
Strategies for analyzing data from video fluorometric monitoring of liquid chromatographic effluents
- Appellof, Davidson
- 1981
(Show Context)
Citation Context |

11 |
Decomposing the time-frequency representation of EEG using non-negative matrix and multi-way factorization,” tech
- Mørup, Hansen, et al.
- 2006
(Show Context)
Citation Context ...s [3, 31, 8, 44]. More recently, these methods have been applied to signal processing [10, 11], image processing [52, 54, 57, 53], data mining [43, 46, 1, 45], scientific computing [6], and elsewhere =-=[26, 36]-=-. Though this work can be applied in a variety of contexts, we concentrate on operations that are common to tensor decompositions, such as the Tucker [51] and canonical decomposition (CANDECOMP), and ... |

10 | Collective sampling and analysis of high order tensors for chatroom communications
- ACAR, ÇAMTEPE, et al.
(Show Context)
Citation Context ...itions have been motivated by applications in chemometrics [3, 31, 8, 44]. More recently, these methods have been applied to signal processing [10, 11], image processing [52, 54, 57, 53], data mining =-=[43, 46, 1, 45]-=-, scientific computing [6], and elsewhere [26, 36]. Though this work can be applied in a variety of contexts, we concentrate on operations that are common to tensor decompositions, such as the Tucker ... |

9 | Analyses and tests of handwritten digit recognition algorithms
- Savas
- 2003
(Show Context)
Citation Context ...itions have been motivated by applications in chemometrics [3, 31, 8, 44]. More recently, these methods have been applied to signal processing [10, 11], image processing [52, 54, 57, 53], data mining =-=[43, 46, 1, 45]-=-, scientific computing [6], and elsewhere [26, 36]. Though this work can be applied in a variety of contexts, we concentrate on operations that are common to tensor decompositions, such as the Tucker ... |

8 | Lathauwer, Blind identification of convolutive MIM systems with 3 sources and 2 - Chen, Petropolu, et al. - 2002 |

8 |
Tensor displacement structures and polyspectral matching,” in Fast Reliable Algorithms for Matrices with Structure
- GRIGORASCU, REGALIA
- 1998
(Show Context)
Citation Context ...sing (4.1) as � U (1) , U (2) ,...,U (N)� · G. Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.sSPARSE AND FACTORED TENSORS 219 As another example, Grigorascu and Regalia =-=[17]-=- emphasize the role of the core tensor in the multiplication by expressing (4.1) as X = U (1) G ⋆ U (2) G ⋆ ··· G ⋆ U (N) , which is called the weighted Tucker product; the unweighted version has G = ... |

8 |
Joint orthomax rotation of the core and component matrices resulting from three-mode principal component analysis
- Kiers
- 1998
(Show Context)
Citation Context ...s trivial. More typically, an alternating least squares (ALS) approach is used for the computation; see [27, 47, 13]. The Tucker decomposition is not unique, but measures can be taken to correct this =-=[20, 21, 22, 48]-=-. Observe that the right-hand side of (2.4) is a Tucker tensor, to be discussed in more detail in section 4. The CANDECOMP/PARAFAC decomposition was simultaneously developed as the canonical decomposi... |

8 |
Use of the properties of the Khatri-Rao product for the computation of Jacobian, Hessian, and gradient of the PARAFAC model under MATLAB, 2005. Private communication
- Tomasi
(Show Context)
Citation Context ... (n) =[v (n) 1 ··· v (n) R ] for n =1,...,N. If we fix everything but V(n) , then solving for it is a linear least squares problem. The pseudoinverse of the Khatri–Rao product W has special structure =-=[7, 49]-=-: W † � = V (N) ⊙···⊙V (n+1) ⊙ V (n−1) ⊙···⊙V (1)� Z † , where � Z = V (1)T V (1)� � ∗...∗ V (n−1)T V (n−1)� � ∗ V (n+1)T V (n+1)� � ∗...∗ V (N)T V (N)� . The least squares solution is given by V (n) ... |

7 | Efficient Representation Scheme for Multidimensional Array Operations
- Lin, Liu, et al.
- 2002
(Show Context)
Citation Context ...ach frontal slice X::k as a sparse matrix in, say, CSC format. The entries are consequently sorted first by the third index and then by the second index. Another idea, proposed by Lin, Liu, and Chung =-=[34, 33]-=-, is to use extended Karnaugh map representation (EKMR). In this case, a three- or four-dimensional tensor is converted to a matrix (see section 2.3) and then stored using a standard sparse matrix sch... |

6 |
From vectors to tensors, Universitext
- Ruíz-Tolosa, Castillo
- 2005
(Show Context)
Citation Context ...ed sparse column (CSC) format [16]. Mathematica, an alternative to MATLAB, also supports multidimensional arrays, and there is a Mathematica package for working with tensors that accompanies the book =-=[41]-=-. In terms of sparse arrays, Mathematica stores its SparseArray’s in compressed sparse row (CSR) format and claims that its format is general enough to describe arbitrary r=1sSPARSE AND FACTORED TENSO... |

6 |
Some additional results on principal components analysis of three-mode data by means of alternating least squares algorithms
- Berghe, Leeuw, et al.
- 1987
(Show Context)
Citation Context ... all n =1,...,N.IfJn = rankn(X) for all n, then the approximation is exact and the computation is trivial. More typically, an alternating least squares (ALS) approach is used for the computation; see =-=[27, 47, 13]-=-. The Tucker decomposition is not unique, but measures can be taken to correct this [20, 21, 22, 48]. Observe that the right-hand side of (2.4) is a Tucker tensor, to be discussed in more detail in se... |

5 |
A.: MultiArray: a c++ library for generic programming with arrays
- Garcia, Lumsdaine
- 2005
(Show Context)
Citation Context ...ld be sorted. It does not appear that HTL supports general tensor multiplication, but it does support inner product, addition, elementwise multiplication, and more. We also briefly mention MultiArray =-=[15]-=-, which provides a general array class template that supports multiarray abstractions and can be used to store dense tensors. Because it directly informs our proposed data structure, related work on s... |

5 |
Some basic techniques for solving sparse systems of linear equations
- Gustavson
- 1971
(Show Context)
Citation Context ...r a matrix A of size I × J with nnz(A) nonzeros, the total storage is 3 ·nnz(A), and the indices are not necessarily presorted. More common are CSR and CSC formats, which appear to have originated in =-=[18]-=-. The CSR format stores three one-dimensional arrays: an array of length nnz(A) with the nonzero values (sorted by row), an array of length nnz(A) with corresponding column indices, and an array of le... |