CP (CANDECOMP/PARAFAC) Decomposition of Order-3 Tensors

Keywords: tensor decomposition; canonical polyadic; CANDECOMP/PARAFAC; alternating least squares; Khatri-Rao product; tensor rank

1. Introduction

A great deal of experimental and simulated data is naturally three-way. A flow stress dataset indexed by strain, strain rate, and temperature is one example; an excitation-emission-sample fluorescence dataset is the classical chemometrics example [1]. Flattening such data into a matrix and applying principal component analysis or the singular value decomposition (see the companion note on SVD) collapses two of the three modes together, discards the coupling structure, and yields factors that are only determined up to an arbitrary rotation.

The polyadic form of a tensor, a sum of outer products, was introduced by Hitchcock in 1927 [2]. It was rediscovered and turned into a practical data-analysis tool twice in 1970: as PARAFAC (parallel factor analysis) by Harshman [3], and as CANDECOMP (canonical decomposition) by Carroll and Chang [4]. The two are the same model, and the neutral modern name canonical polyadic (CP) decomposition is now common following the survey of Kolda and Bader [5]. Its defining property, in contrast to matrix factorization, is that under Kruskal's condition [6] the decomposition is unique up to the unavoidable permutation and scaling of components, without any orthogonality assumption being imposed. This uniqueness is why the individual CP components can be read as physical modes rather than as an arbitrary basis.

The rest of this note is organized as follows. Section 2 defines the outer product, rank-1 tensors, tensor rank, the CP model, matricization, and the Khatri-Rao product. Section 3 derives the CP-ALS algorithm and its efficient implementation. Section 4 discusses rank selection. Section 5 states Kruskal's uniqueness condition. Section 6 works a numerical example by hand. Section 7 connects the model to the constitutive-modeling work on this site, and Section 8 concludes.

2. Basic Definitions

2.1 Outer product and rank-1 tensors

Let \( \mathcal{X} \in \mathbb{R}^{I \times J \times K} \) be an order-3 tensor with entries \( x_{ijk} \). Given three vectors \( a \in \mathbb{R}^{I} \), \( b \in \mathbb{R}^{J} \), \( c \in \mathbb{R}^{K} \), their outer product is the \( I \times J \times K \) rank-1 tensor

A tensor is rank-1 if it can be written as a single such outer product. This is the direct three-way analogue of a rank-1 matrix \( a\,b^{\mathsf T} \).

2.2 The CP model and tensor rank

The CP decomposition approximates \( \mathcal{X} \) as a sum of \( R \) rank-1 tensors:

The columns of the factor matrices \( A = [a_1,\dots,a_R] \in \mathbb{R}^{I \times R} \), \( B \in \mathbb{R}^{J \times R} \), and \( C \in \mathbb{R}^{K \times R} \) are usually normalized to unit length, and the scalar weights are collected into \( \lambda = (\lambda_1,\dots,\lambda_R) \). The tensor rank of \( \mathcal{X} \) is the smallest \( R \) for which Eq. (2) holds exactly. Unlike the matrix case, computing the exact tensor rank is NP-hard [5,7], and the maximum possible rank can exceed the largest dimension, so in practice one fixes a small target \( R \) and computes a low-rank approximation.

2.3 Matricization and the Khatri-Rao product

Unfolding (matricizing) \( \mathcal{X} \) along its three modes gives the matrices \( X_{(1)} \in \mathbb{R}^{I \times JK} \), \( X_{(2)} \in \mathbb{R}^{J \times IK} \), and \( X_{(3)} \in \mathbb{R}^{K \times IJ} \). The Khatri-Rao product of two matrices \( U \in \mathbb{R}^{I \times R} \) and \( V \in \mathbb{R}^{J \times R} \), written \( U \odot V \in \mathbb{R}^{IJ \times R} \), is the column-wise Kronecker product: its \( r \)-th column is \( u_r \otimes v_r \). With this notation the CP model of Eq. (2) has the three equivalent matricized forms [5]

where \( \Lambda = \mathrm{diag}(\lambda) \); when the weights are absorbed into one factor matrix, \( \Lambda \) can be dropped. The Khatri-Rao product obeys the Gram identity

where \( * \) is the elementwise (Hadamard) product. Because \( \langle u_r \otimes v_r,\, u_s \otimes v_s \rangle = \langle u_r, u_s \rangle \langle v_r, v_s \rangle \), the identity collapses a \( JK \times R \) inner product into the Hadamard product of two \( R \times R \) matrices. This is the single fact that makes CP-ALS cheap.

3. The Alternating Least Squares Algorithm (CP-ALS)

3.1 The subproblem

CP fitting minimizes the Frobenius-norm reconstruction error

This objective is non-convex in the three factor matrices jointly, but it is quadratic, hence convex, in any single factor matrix when the other two are held fixed. Fixing \( B \) and \( C \) and using the mode-1 matricization from Eq. (3), the problem for \( A \) becomes an ordinary linear least squares problem

whose normal equations are \( A\,(C \odot B)^{\mathsf{T}}(C \odot B) = X_{(1)}(C \odot B) \). Applying the Gram identity of Eq. (4) to the coefficient matrix gives the CP-ALS update in closed form:

where \( {}^{\dagger} \) is the Moore-Penrose pseudoinverse of the \( R \times R \) Gram matrix. The term \( X_{(1)}(C \odot B) \) is the matricized-tensor-times-Khatri-Rao-product, abbreviated MTTKRP; it is the dominant cost of the algorithm and never requires forming the large \( JK \times R \) matrix \( C \odot B \) explicitly. By symmetry, the updates for the other two modes are

3.2 Normalization and stopping

After each update the columns of the just-updated factor matrix are rescaled to unit norm and the removed norms are stored in \( \lambda \). This prevents one factor from growing while its partner shrinks (a scaling drift that leaves the product unchanged but degrades conditioning). One outer sweep performs the three updates of Eqs. (7) and (8), then evaluates the relative error \( \|\mathcal{X} - [\![A,B,C]\!]\|_F / \|\mathcal{X}\|_F \). Iteration stops when the change in this error between sweeps falls below a tolerance, or a maximum sweep count is reached.

3.3 The general algorithm

Here (.) denotes the Khatri-Rao product and [ ]^+ the pseudoinverse. Each sweep costs \( O(IJK \cdot R) \), dominated by the three MTTKRPs; the \( R \times R \) Gram products and their pseudoinverses are negligible by comparison. This is the algorithm implemented as cp_als in the MATLAB Tensor Toolbox [8] and as parafac in the N-way toolbox, and it is the exact routine used in the Part 2 flow stress paper below.

When the factors are additionally constrained to be non-negative, the closed-form pseudoinverse update of Eq. (7) is replaced by a non-negative least squares solve or a multiplicative update; that variant is the subject of the separate non-negative CP note.

4. Choosing the Rank

Because the true rank is not known in advance and is expensive to compute, \( R \) is chosen empirically. Three practical guides are used together. First, the residual error curve: fit the tensor at \( R = 1, 2, 3, \dots \) and plot the relative reconstruction error; the rank is taken at the knee, beyond which extra components buy little accuracy. Second, the core consistency diagnostic (CORCONDIA) of Bro and Kiers [9], which measures how close the fitted model is to an ideal trilinear one and drops sharply once \( R \) exceeds the trilinear rank. Third, interpretability and stability: an over-large \( R \) causes degeneracy, in which two components become nearly collinear with large equal-and-opposite weights, a symptom that the model is fitting noise rather than structure. In the flow stress work the rank is justified by the residual-error curve together with an accuracy target expressed in MAPE.

5. Uniqueness: Kruskal's Condition

The reason CP is preferred to matrix factorization for interpretation is its essential uniqueness. Let \( k_M \) denote the Kruskal rank (k-rank) of a matrix \( M \), the largest number \( k \) such that every set of \( k \) columns of \( M \) is linearly independent. Kruskal's sufficient condition [6] for the rank-\( R \) CP decomposition to be unique is

When Eq. (9) holds, any two exact decompositions \( [\![A,B,C]\!] = [\![A',B',C']\!] \) are related by a single common permutation of the components and a rescaling of each component that leaves the product unchanged. There is no rotational freedom, in sharp contrast to the SVD, whose factors can be mixed by any orthogonal rotation once the singular values are equal. The unavoidable permutation and scaling mean that two decompositions must be aligned (for example by sorting on \( \lambda \)) before their components are compared.

6. Worked Example

Take \( I = 3 \), \( J = 2 \), \( K = 2 \), and build an exact rank-2 tensor from

with both weights equal to 1, so \( \mathcal{X} = a_1 \circ b_1 \circ c_1 + a_2 \circ b_2 \circ c_2 \). Its two mode-3 slices are

Running CP-ALS at \( R = 2 \) from a random start, the sweeps converge to a relative error below \( 10^{-10} \) and recover, after column normalization and component alignment,

The recovered factors match the construction to within numerical noise, but the algorithm output them in a data-driven order that need not match the construction order, and each component is recovered only up to the sign-free scaling absorbed into \( \hat\lambda \). Both facts are direct manifestations of the permutation-and-scaling indeterminacy of Section 5. Fitting the same tensor at \( R = 1 \) leaves a large residual, and at \( R = 3 \) the third weight collapses toward zero, so the residual-error curve of Section 4 correctly identifies the rank as 2.

7. Use in the Flow Stress Work

8. Conclusions

1. CP decomposition writes an order-3 tensor as a sum of rank-1 outer products. Its matricized forms \( X_{(1)} = A(C \odot B)^{\mathsf T} \), with the Khatri-Rao Gram identity \( (U \odot V)^{\mathsf T}(U \odot V) = (U^{\mathsf T}U)*(V^{\mathsf T}V) \), turn the trilinear fit into a sequence of ordinary least squares subproblems.
2. The CP-ALS algorithm solves each subproblem in closed form via one MTTKRP and one \( R \times R \) pseudoinverse, followed by column normalization, at a per-sweep cost of \( O(IJK\,R) \).
3. Tensor rank is NP-hard to compute exactly, so \( R \) is chosen empirically from the residual-error knee, the core consistency diagnostic, and the absence of degeneracy.
4. Under Kruskal's condition the decomposition is unique up to component permutation and scaling, with no orthogonality imposed, which is what allows its components to be interpreted as physical modes.

References

1. Bro, R. (1997). PARAFAC. Tutorial and applications. Chemometrics and Intelligent Laboratory Systems, 38(2), 149-171.
2. Hitchcock, F. L. (1927). The expression of a tensor or a polyadic as a sum of products. Journal of Mathematics and Physics, 6(1-4), 164-189.
3. Harshman, R. A. (1970). Foundations of the PARAFAC procedure: Models and conditions for an explanatory multi-modal factor analysis. UCLA Working Papers in Phonetics, 16, 1-84.
4. Carroll, J. D., & Chang, J. J. (1970). Analysis of individual differences in multidimensional scaling via an N-way generalization of the Eckart-Young decomposition. Psychometrika, 35(3), 283-319.
5. Kolda, T. G., & Bader, B. W. (2009). Tensor decompositions and applications. SIAM Review, 51(3), 455-500.
6. Kruskal, J. B. (1977). Three-way arrays: rank and uniqueness of trilinear decompositions, with application to arithmetic complexity and statistics. Linear Algebra and Its Applications, 18(2), 95-138.
7. Hastad, J. (1990). Tensor rank is NP-complete. Journal of Algorithms, 11(4), 644-654.
8. Bader, B. W., & Kolda, T. G. (2007). Efficient MATLAB computations with sparse and factored tensors. SIAM Journal on Scientific Computing, 30(1), 205-231.
9. Bro, R., & Kiers, H. A. L. (2003). A new efficient method for determining the number of components in PARAFAC models. Journal of Chemometrics, 17(5), 274-286.

Back to Technical Notes