Non-negative CP Decomposition of Order-3 Tensors

Keywords: tensor decomposition; CP/PARAFAC; non-negativity constraint; alternating least squares; Khatri-Rao product

1. Introduction

Many experimental or simulated datasets naturally have a three-way (or higher) multi-factor structure. A flow stress tensor measured over temperature × strain × strain rate is one example; fluorescence spectroscopy data over wavelength × time × space is another [1]. Simply flattening such data into a matrix and applying principal component analysis (PCA) or non-negative matrix factorization (NMF) discards the coupling structure among the factors, and matrix-factorization factors are generally non-unique: extra constraints such as orthogonality are needed to fix them.

CP decomposition, proposed independently by Harshman [2] and by Carroll and Chang [3], performs low-rank additive modeling directly on the tensor's original three-way structure. Under mild conditions its decomposition is inherently unique, with no orthogonality assumption required [4], which makes it better suited to explaining the data-generating mechanism through physically meaningful components. When the data itself is non-negative, as concentration, intensity, and stress all are, further imposing a non-negativity constraint (NN-CP) yields additive, interpretable factors free of cancelling negative components. This is the same advantage NMF has over PCA [5], and this three-way generalization of the non-negativity constraint is the direct motivation behind early work such as Welling and Weber [6].

The remainder of this note is organized as follows. Section 2 gives the general mathematical definition of the CP model, outer-product notation, matricization, and the Khatri-Rao product. Section 3 formulates NN-CP decomposition as a non-negativity-constrained least squares problem and states its KKT optimality conditions. Section 4 derives the alternating non-negative least squares algorithm based on multiplicative updates, writes out explicit matrix forms for Rank-1, Rank-2, and Rank-3, and gives the general algorithm for arbitrary Rank-R. Section 5 reviews Kruskal's uniqueness condition. Section 6 gives hand-verifiable numerical examples. Section 7 describes the reference implementation and discusses convergence and complexity. Section 8 concludes.

2. Basic Definitions of CP Decomposition

2.1 Order-3 tensors and the outer product

Let \( \mathcal{X} \in \mathbb{R}_{\ge 0}^{I \times J \times K} \) be an order-3 (three-way) non-negative tensor with entries \( x_{ijk} \), where \( i = 1,\dots,I \), \( j = 1,\dots,J \), \( k = 1,\dots,K \) correspond respectively to the tensor's 1st, 2nd, and 3rd modes. Given three vectors \( a \in \mathbb{R}^{I} \), \( b \in \mathbb{R}^{J} \), \( c \in \mathbb{R}^{K} \), their outer product is defined as an \( I \times J \times K \) rank-1 tensor:

2.2 The general form of the CP model (Rank-R)

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

where \( A = [a_1,\dots,a_R] \in \mathbb{R}^{I \times R} \), \( B = [b_1,\dots,b_R] \in \mathbb{R}^{J \times R} \), and \( C = [c_1,\dots,c_R] \in \mathbb{R}^{K \times R} \) are the mode-1, mode-2, and mode-3 factor matrices, and \( \lambda = (\lambda_1,\dots,\lambda_R) \) is the weight vector. When \( R \) is the smallest value for which \( \mathcal{X} \) can be represented exactly, \( R \) is the (non-negative) rank of the tensor. Non-negative CP decomposition requires

Without loss of generality, an implementation may absorb \( \lambda_r \) directly into the columns of one factor matrix, and separate \( \lambda \) out again after decomposition by column normalization.

2.3 Matricization (unfolding) and the Khatri-Rao product

To turn the trilinear model of Eq. (2) into a least squares problem that matrix algebra can handle, we introduce the tensor's matricization. Unfolding \( \mathcal{X} \) along mode-1, mode-2, and mode-3 gives 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} \) respectively.

Next, define the Khatri-Rao product (column-wise Kronecker product) of two matrices \( U \in \mathbb{R}^{I \times R} \) and \( V \in \mathbb{R}^{J \times R} \) as \( U \odot V \in \mathbb{R}^{IJ \times R} \), whose \( r \)-th column is \( u_r \otimes v_r \). Using it, the three matricized forms of the CP model can be written, following the convention of Kolda and Bader [4], as:

The Khatri-Rao product satisfies a key identity that makes the algorithm efficient. For any \( U, V \) with the same number of columns,

where \( * \) denotes the elementwise (Hadamard) product. This holds 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 \). Equation (4) reduces an inner-product computation of size \( IJ \times IJ \) to the Hadamard product of two small \( R \times R \) matrices, which is the core reason the alternating least squares algorithm below can be implemented efficiently.

3. The Optimization Problem Under Non-negativity Constraints

Expressing Eq. (2) as a least squares objective under the Frobenius norm, with non-negativity constraints imposed:

Equation (5) is non-convex jointly in the three factor matrices, but fixing any two of them and solving only for the third is a convex non-negative least squares (NNLS) problem. This is exactly the origin of the alternating least squares (ALS) idea. Taking the case of fixing \( B, C \) and solving for \( A \), Eq. (5) is rewritten in matrix form using Eq. (3):

This is a standard matrix NNLS problem. Its KKT optimality conditions are that the gradient satisfies \( \nabla_A f \ge 0 \) elementwise, and that the complementary slackness condition \( A_{ir}\,[\nabla_A f]_{ir} = 0 \) holds for all \( (i,r) \). The gradient of Eq. (6) with respect to \( A \) is:

where the identity of Eq. (4) has been used to reduce \( (C \odot B)^{\mathsf{T}} (C \odot B) \) to an \( R \times R \) Hadamard product, avoiding explicit construction of a matrix product of size \( JK \times R \).

4. The Alternating Non-negative Least Squares (NN-ALS) Algorithm

4.1 Derivation of the multiplicative update rule

The multiplicative update technique proposed by Lee and Seung [5,7] for non-negative matrix factorization generalizes directly to Eq. (7). Split the gradient into two non-negative parts:

Both terms are elementwise non-negative whenever \( A, B, C, \mathcal{X} \ge 0 \). Define the multiplicative update:

where \( \circledast \) and the division are both elementwise operations, and a small positive \( \varepsilon \) is added to the denominator to avoid division by zero. Equation (9) has two important properties:

1. Non-negativity is preserved automatically. As long as the initial \( A^{(0)} \ge 0 \), the ratio \( \nabla_A^{-} / \nabla_A^{+} \ge 0 \), so the product remains non-negative. No extra projection or clipping step is needed.
2. Fixed points satisfy the KKT conditions. If the update converges to a fixed point \( A = A \circledast (\nabla_A^{-}/\nabla_A^{+}) \), then for every \( (i,r) \) we have \( A_{ir}\nabla^{+}_{A,ir} = A_{ir}\nabla^{-}_{A,ir} \), that is, \( A_{ir}[\nabla_A f]_{ir} = 0 \), which is exactly the complementary slackness condition of the NNLS problem. It can also be shown, via an auxiliary-function argument, that this update monotonically decreases the objective in Eq. (5); the proof follows the same reasoning as Lee and Seung's original NMF multiplicative update and is omitted here.

The updates for \( B \) and \( C \) follow by symmetry, based on the matricizations \( X_{(2)} \) and \( X_{(3)} \) respectively:

In each outer sweep, the algorithm performs the three updates of Eqs. (9) and (10) in turn, then recomputes the reconstructed tensor \( \hat{\mathcal{X}} = [\![A,B,C]\!] \) and the relative error \( \|\mathcal{X} - \hat{\mathcal{X}}\|_F / \|\mathcal{X}\|_F \), until the change in this error falls below a threshold or the maximum number of iterations is reached.

Below, the \( R \times R \) Gram-matrix portion of Eqs. (9) and (10) is expanded explicitly for \( R = 1, 2, 3 \), to show what the abstract rank-R formula looks like concretely at small ranks.

4.2 Explicit derivation for Rank-1

When \( R = 1 \), the factor matrices \( A = a \), \( B = b \), \( C = c \) degenerate to column vectors, \( C^{\mathsf{T}}C = \|c\|^2 \) and \( B^{\mathsf{T}}B = \|b\|^2 \) are both scalars, and the Khatri-Rao product \( C \odot B \) degenerates to the ordinary Kronecker product of vectors \( c \otimes b \). Equation (9) simplifies to an elementwise scalar update:

That is,

where the numerator is the bilinear contraction of \( \mathcal{X} \) with \( b \) and \( c \) over modes 2 and 3, and the denominator is the product of the two squared norms. This is exactly the form of tensor power iteration for updating \( a \). For an exactly rank-1 non-negative \( \mathcal{X} \), Perron-Frobenius theory guarantees the existence of non-negative \( a, b, c \) such that the iteration converges to the exact solution, and in principle a single alternating update together with normalization suffices, with no crosstalk between components to handle as in the \( R \gt 1 \) case.

4.3 Explicit derivation for Rank-2

When \( R = 2 \), \( C^{\mathsf{T}}C \) and \( B^{\mathsf{T}}B \) are \( 2 \times 2 \) symmetric matrices:

The \( 2 \times 2 \) Hadamard product inside the denominator matrix \( A[(C^{\mathsf{T}}C) * (B^{\mathsf{T}}B)] \) is

Its off-diagonal entry \( (c_1 \cdot c_2)(b_1 \cdot b_2) \) is precisely the source of crosstalk between the two rank-1 components. When the \( b \) and \( c \) factors of the two components are mutually orthogonal, this matrix degenerates to a diagonal matrix and Eq. (9) reduces to two independent Rank-1 updates as in Eq. (11). In general the off-diagonal entry is non-zero, the two components are coupled during the iteration, and many alternations are required to converge, which is the slower convergence observed in the numerical example of Section 6.2.

4.4 Explicit derivation for Rank-3

For \( R = 3 \), the Gram matrices above extend to \( 3 \times 3 \):

The \( 3 \times 3 \) Hadamard product \( (C^{\mathsf{T}}C) * (B^{\mathsf{T}}B) \) in the denominator of Eq. (9) now has six independent off-diagonal entries, since the matrix is symmetric, each corresponding to the coupling strength between one pair of components. Computationally, relative to Rank-2 and Rank-1, only the Gram-matrix size grows from \( 1 \times 1 \) to \( 2 \times 2 \) to \( 3 \times 3 \). The cost of the MTTKRP (matricized tensor times Khatri-Rao product, the term \( X_{(1)}(C \odot B) \)) grows linearly with \( R \), so the good scalability of the ALS algorithm with rank is preserved overall.

4.5 The general Rank-R algorithm

Combining Eqs. (9) and (10), the NN-CP alternating non-negative least squares algorithm for an arbitrary rank \( R \) can be summarized as follows.

Here (*) denotes the elementwise product and (.) denotes the Khatri-Rao product. The rank parameter is exactly \( R \), and the same code applies to \( R = 1, 2, 3 \) or any larger \( R \) without a separate implementation for each rank.

5. Uniqueness: Kruskal's Condition

Unlike matrix decompositions such as SVD or NMF, CP decomposition is essentially unique under mild conditions, up to permutation and per-component rescaling, without requiring extra constraints such as orthogonality. The classical sufficient condition of Kruskal [8] is:

where \( k_M \) denotes the Kruskal rank (k-rank) of matrix \( M \): the largest integer \( k \) such that every \( k \) columns of \( M \) are linearly independent. When Eq. (15) is satisfied, if \( [\![A,B,C]\!] = [\![A',B',C']\!] \), then there must exist a permutation matrix \( \Pi \) and diagonal scaling matrices \( D_1, D_2, D_3 \) satisfying \( D_1 D_2 D_3 = I \) such that \( A' = A \Pi D_1 \), \( B' = B \Pi D_2 \), \( C' = C \Pi D_3 \).

This property is the most important theoretical advantage of CP decomposition over matrix decomposition: when Eq. (15) holds, the decomposed components are determined by the data-generating mechanism itself, rather than by an algorithmic or artificially chosen constraint. The component-order permutation observed in the Rank-3 numerical example of Section 6.3 is a direct manifestation of this permutation indeterminacy.

6. Numerical Examples

The three examples below are all generated and verified by first constructing an exact rank-R tensor from given non-negative vectors, then decomposing it with the NN-ALS algorithm of Section 4.5. All tensors are constructed synthetically and rely on no external database or dataset.

6.1 Rank-1 example

Take \( I = J = K = 2 \) and construct the vectors

Using Eq. (1) to construct the exact rank-1 tensor \( \mathcal{X} = a \circ b \circ c \), its two slices, fixing \( k \) along mode-3, are:

Iterating Eq. (11) from a random non-negative initialization of \( a, b, c \) at rank 1, after convergence and column normalization with weight \( \hat{\lambda} \), we obtain:

The final relative reconstruction error is \( \|\mathcal{X} - \hat{\mathcal{X}}\|_F / \|\mathcal{X}\|_F \approx 2.4 \times 10^{-14} \), that is, exact recovery to floating-point precision. This is consistent with the theoretical analysis of Section 4.2: decomposing an exactly rank-1 non-negative tensor is essentially a power-iteration problem, which converges quickly and uniquely, with no permutation ambiguity between components since there is only one component.

6.2 Rank-2 example

Take \( I = 3 \), \( J = 2 \), \( K = 2 \), with the vectors of the two rank-1 components given by

with both component weights taken as 1. Constructing \( \mathcal{X} = a_1 \circ b_1 \circ c_1 + a_2 \circ b_2 \circ c_2 \), its two mode-3 slices are:

After 8000 iterations at rank 2, the recovered results, component-aligned and column-normalized, are:

The recovered values agree fully with the true normalized values, with numerical noise below \( 10^{-4} \) and a final relative reconstruction error of about \( 4.1 \times 10^{-5} \). Compared with the Rank-1 example, the number of iterations required for convergence is markedly larger, consistent with the crosstalk between components discussed via Eq. (13) in Section 4.3.

6.3 Rank-3 example

Take \( I = 4 \), \( J = 3 \), \( K = 3 \), with three mutually non-orthogonal rank-1 components carrying different weights 2.0, 1.5, and 1.0:

After 15000 iterations, the rank-3 decomposition converges with a relative reconstruction error of about \( 2.7 \times 10^{-5} \). The recovered weights are \( \hat{\lambda} = (2.8284,\, 5.6568,\, 4.2426) \), and the internal component ordering produced by the algorithm differs from the construction order. As discussed under permutation indeterminacy in Section 5, the actual correspondence is:

Checking the three sets of normalized factor vectors one by one confirms numerical agreement with the corresponding true components, with differences below \( 10^{-4} \). This behaviour, components reordered but numerically exactly recovered, is a direct manifestation in a numerical experiment of the theoretical property that CP decomposition is unique only up to permutation and scaling. It also shows that comparing two decomposition results requires component alignment first, for example sorting by \( \hat{\lambda} \) or matching by inter-factor correlation, rather than comparing directly by component index.

7. Implementation and Convergence

The reference implementation is a pure-NumPy routine that relies on no dedicated tensor-decomposition library and loads no ready-made dataset. It consists of the following parts:

• Khatri-Rao product of an arbitrary number of matrices.
• Unfolding: matricization of an order-3 tensor along mode-1, mode-2, or mode-3.
• Reconstruction: rebuilding the full tensor \( [\![A,B,C]\!] \) from the factor matrices.
• Normalization: absorbing column norms into the weight \( \lambda \), leaving the factor matrices with unit column norm.
• NN-ALS solver: a direct implementation of the pseudocode in Section 4.5, with the rank parameter supporting any \( R \ge 1 \).
• Synthetic tensor generator: builds an exact rank-R tensor from given vectors and calls the solver to verify the decomposition result. The three examples in Section 6 and one additional \( R = 6 \) example are all generated this way.

Convergence and complexity. The multiplicative updates of Eqs. (9) and (10) can be shown to monotonically decrease the objective in Eq. (5), via a Lee-Seung-type auxiliary-function argument, but the convergence rate is first-order (linear), generally slower than implementations based on exact non-negative least squares such as the Lawson-Hanson active-set method [10], or on hierarchical alternating least squares [9]. This is consistent with the Rank-2 and Rank-3 examples of Section 6, where thousands of sweeps are needed to reach a relative error of \( 10^{-4} \) to \( 10^{-5} \), whereas Rank-1 converges to machine precision in effectively one sweep.

The dominant computational cost of each sweep is the product of a matricized tensor with a Khatri-Rao product (MTTKRP), with complexity about \( O(IJK \cdot R) \), growing linearly with the total number of tensor elements and with the rank \( R \). The cost of the \( R \times R \) Gram-matrix products is negligible by comparison. For faster convergence, the multiplicative updates in this note can be replaced with HALS or a projected-gradient or active-set NNLS subproblem solver, while the overall algorithmic framework of matricization plus alternating sweeps remains unchanged.

8. Conclusions

1. CP decomposition of an order-3 tensor represents the data as a sum of rank-1 tensors. Its matricized forms such as \( X_{(1)} = A(C \odot B)^{\mathsf{T}} \), together 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 optimization problem into a sequence of efficiently solvable non-negative least squares subproblems.
2. Under the non-negativity constraint, Lee-Seung-type multiplicative updates monotonically decrease the reconstruction error while guaranteeing non-negativity of the factor matrices, and their fixed points satisfy exactly the KKT complementary slackness conditions of the NNLS problem. This update rule has a unified form for any rank \( R \), with Rank-1, Rank-2, and Rank-3 being special cases that differ only in the size of the \( R \times R \) Gram matrix.
3. The numerical examples show that an exact rank-1 tensor is recovered in effectively one convergence step to floating-point precision. As the rank increases, off-diagonal Gram terms couple the components and the number of iterations needed for convergence grows significantly. The recovered component order may differ from the construction order, a direct manifestation of the theoretical property that CP decomposition is unique only up to permutation and scaling.
4. A minimal, dependency-free NN-CP implementation covers decomposition needs from Rank-1 to any Rank-R through a single rank parameter, with no separate code path per rank.

References

1. Bro, R. (1997). PARAFAC. Tutorial and applications. Chemometrics and Intelligent Laboratory Systems, 38(2), 149-171.
2. 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.
3. 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.
4. Kolda, T. G., & Bader, B. W. (2009). Tensor decompositions and applications. SIAM Review, 51(3), 455-500.
5. Lee, D. D., & Seung, H. S. (1999). Learning the parts of objects by non-negative matrix factorization. Nature, 401(6755), 788-791.
6. Welling, M., & Weber, M. (2001). Positive tensor factorization. Pattern Recognition Letters, 22(12), 1255-1261.
7. Lee, D. D., & Seung, H. S. (2001). Algorithms for non-negative matrix factorization. Advances in Neural Information Processing Systems, 13, 556-562.
8. 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.
9. Cichocki, A., & Phan, A. H. (2009). Fast local algorithms for large scale nonnegative matrix and tensor factorizations. IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences, E92-A(3), 708-721.
10. Lawson, C. L., & Hanson, R. J. (1974). Solving Least Squares Problems. Prentice-Hall.

Back to Technical Notes