Least-Squares CP Fitting by Alternating Least Squares

Why Fit a CP Model Alternately?

A CP model represents a tensor as a sum of rank-one tensors. Fitting all factor matrices simultaneously is a nonconvex optimization problem: products of unknown factor entries appear in the reconstruction. Yet if all but one factor matrix are held fixed, the remaining factor enters linearly. Alternating least squares (ALS) uses this fact, solving a sequence of ordinary least-squares problems, one factor matrix at a time.

Each block update can be derived exactly, but that does not make the joint problem convex or guarantee a globally best fit. Understanding the update, its rank conditions, and the model’s indeterminacies is essential both for correct computation and for interpreting the fitted factors.

Definition: Let \(\mathcal{X}\in\mathbb{R}^{I_1\times\cdots\times I_N}\) be a complete data tensor, and choose a target CP rank \(R\). A CP approximation has factor matrices \(A^{(n)}\in\mathbb{R}^{I_n\times R}\) and entries \[ \hat{x}_{i_1,\ldots,i_N} =\sum_{r=1}^{R}\prod_{n=1}^{N}a^{(n)}_{i_n r}. \] The least-squares CP objective is \(F(A^{(1)},\ldots,A^{(N)})=\frac12\|\mathcal{X}-\hat{\mathcal{X}}\|_F^2\), where \(\|\cdot\|_F\) is the Frobenius norm.

Unfolding the Model and Deriving an Update

For mode \(n\), let \(\mathcal{X}_{(n)}\) denote the mode-\(n\) unfolding: its rows index \(i_n\), and its columns index the remaining index tuple in the chosen unfolding order. For matrices with the same number of columns, their Khatri–Rao product is the columnwise Kronecker product. Define

Z^{(n)} = A^{(N)}\odot\cdots\odot A^{(n+1)} \odot A^{(n-1)}\odot\cdots\odot A^{(1)} \in\mathbb{R}^{J_n\times R}, \qquad J_n=\prod_{m\ne n} I_m ,

with the row ordering chosen to agree with the columns of \(\mathcal{X}_{(n)}\). Thus \(Z^{(n)}\)'s \(r\)th column contains the products of entries from component \(r\) in all modes other than \(n\). The unfolded CP approximation is

\hat{\mathcal{X}}_{(n)}=A^{(n)}(Z^{(n)})^T.

When the other factor matrices are fixed, write \(X_n=\mathcal{X}_{(n)}\), \(A=A^{(n)}\), and \(Z=Z^{(n)}\). The update minimizes \(\frac12\|X_n-AZ^T\|_F^2\) over \(A\). The dimensions are \(X_n:I_n\times J_n\), \(A:I_n\times R\), and \(Z:J_n\times R\), so \(AZ^T\) has the same dimensions as \(X_n\).

Theorem (mode-wise least-squares update): For fixed \(Z\), a minimizer of \(\|X_n-AZ^T\|_F^2\) satisfies \(A(Z^TZ)=X_nZ\). If \(Z\) has full column rank, the unique minimizer is \(A=X_nZ(Z^TZ)^{-1}\). Without full column rank, a minimum-Frobenius-norm minimizer is \(A=X_n((Z^T)^\dagger)=X_n(Z^\dagger)^T\), where \(\dagger\) is the Moore–Penrose pseudoinverse.

To derive the normal equations, differentiate the squared Frobenius norm with respect to \(A\). Its gradient is \(2(AZ^TZ-X_nZ)\); setting this equal to zero gives the stated equation. If \(Z\) has full column rank, \(Z^TZ\) is invertible, and right multiplication by its inverse gives the unique solution.

For the rank-deficient case, set \(B=Z^T\), so the problem is \(\min_A\|AB-X_n\|_F^2\). The standard minimum-norm solution to this matrix least-squares problem is \(A=X_nB^\dagger\). Here \(B^\dagger=(Z^\dagger)^T\), giving the formula above. In particular, \(Z^\dagger\) has dimensions \(R\times J_n\), so \((Z^\dagger)^T\) has dimensions \(J_n\times R\); the product \(X_n(Z^\dagger)^T\) is well-defined and has dimensions \(I_n\times R\). Rank deficiency can make the minimizer nonunique, even though the minimum-norm choice is specified. \(\square\)

In the full-rank case, \(Z^TZ\) has a useful structure. Since the \(r\)th column of \(Z\) is a Kronecker product of the other mode-\(r\) factor columns, inner products factor across modes:

Z^TZ=\underset{m\ne n}{\ast}\bigl((A^{(m)})^T A^{(m)}\bigr),

where \(\ast\) denotes entrywise (Hadamard) multiplication. Its \((r,s)\) entry is the product of the corresponding inner products across modes \(m\ne n\). In computation, solve the normal equations as a linear system rather than explicitly forming a matrix inverse; if the system is ill-conditioned, the solve can still be numerically unstable.

Worked Factor Updates and Indeterminacies

Worked Example: A Rank-One Update for a Three-Way Tensor

Consider a \(2\times2\times2\) tensor with mode-1 unfolding \[ X_1= \begin{bmatrix} 1&2&3&4\\ 2&0&1&3 \end{bmatrix}. \] For a rank-one fit, hold the mode-2 factor \(b=(1,2)^T\) and mode-3 factor \(c=(1,1)^T\) fixed. Using columns ordered by the mode-2, then mode-3 index, the Khatri–Rao design column is \[ z=b\otimes c=(1,1,2,2)^T. \] The least-squares update is \(a=X_1z/(z^Tz)\). The denominator is \(1^2+1^2+2^2+2^2=10\), and the two rowwise numerators are \[ 1(1)+2(1)+3(2)+4(2)=17,\qquad 2(1)+0(1)+1(2)+3(2)=10. \] Therefore \(a=(17/10,1)^T=(1.7,1)^T\). Each entry is the least-squares coefficient for one row of the unfolding, using the same fixed design vector \(z\).

Worked Example: A Rank-Deficient Design and the Minimum-Norm Update

Take \(X_1=(2,0,2)\), a \(1\times3\) matrix, and let \[ Z= \begin{bmatrix} 1&1\\ 0&0\\ 1&1 \end{bmatrix} \in\mathbb{R}^{3\times2}. \] Its two columns are identical, so the two coefficients in \(A=(a_1,a_2)\) cannot be distinguished: \(AZ^T=(a_1+a_2,0,a_1+a_2)\). The least-squares fit equals \(X_1\) whenever \(a_1+a_2=2\), and among these solutions the smallest Euclidean norm occurs at \(a_1=a_2=1\).

The pseudoinverse confirms this result. Write \(Z=uv^T\), where \(u=(1,0,1)^T\) and \(v=(1,1)^T\). For a rank-one matrix, \(Z^\dagger=vu^T/(\|u\|^2\|v\|^2)\), so \[ Z^\dagger= \frac14 \begin{bmatrix} 1&0&1\\ 1&0&1 \end{bmatrix}, \qquad (Z^\dagger)^T= \frac14 \begin{bmatrix} 1&1\\ 0&0\\ 1&1 \end{bmatrix}. \] The dimensions are \(Z^\dagger:2\times3\) and \((Z^\dagger)^T:3\times2\). Hence \[ X_1(Z^\dagger)^T =(2,0,2)\frac14 \begin{bmatrix} 1&1\\ 0&0\\ 1&1 \end{bmatrix} =(1,1). \] The minimum-norm update is therefore well-defined and gives one of the many exact minimizers.

Definition (scaling indeterminacy): In a CP component, multiplying its mode-\(n\) factor column by a nonzero scalar \(c_n\) in every mode leaves the component unchanged whenever \(\prod_{n=1}^N c_n=1\). A common permutation of the \(R\) columns across all factor matrices also leaves the fitted tensor unchanged.

Worked Example: The Same Component With Different Factor Scales

For three modes, let the factor columns be \(a=(2,4)^T\), \(b=(1,-1)^T\), and \(c=(2,1)^T\). The entry at index \((1,1,1)\) of their outer product is \(2\cdot1\cdot2=4\). Replace \(a\) by \(3a\) and \(b\) by \(b/3\), leaving \(c\) unchanged. The same entry becomes \((3\cdot2)(1/3)\cdot2=4\), and every other entry is unchanged because the two rescaling factors cancel in every component product. The factor columns have changed, but the represented tensor has not.

This freedom means raw factor magnitudes are not uniquely interpretable. A common convention is to normalize factor columns and store the component magnitudes separately, while remembering that normalization fixes a representation convention rather than changing the fit.

What ALS Guarantees—and What It Does Not

ALS initializes factor matrices and cycles through the modes. At each mode it minimizes the objective over that factor while holding all other factors fixed. The resulting sequence of objective values cannot increase, provided each update is an exact least-squares minimizer.

Theorem (monotone descent of exact CP-ALS): Suppose every mode update is an exact minimizer of its conditional least-squares problem, and the objective is evaluated on the same complete tensor. After any individual update, the objective is no greater than immediately before that update.

To prove the claim, consider an update of mode \(n\), with all other factors fixed. Before the update, the current \(A^{(n)}\) is one feasible choice for the conditional problem. The updated factor is a minimizer over all choices, so its residual norm is no greater than the residual norm from the previous choice. Squaring and multiplying by \(1/2\) preserves this inequality. The other factors do not change during this block update, so this is also the comparison of the full objective before and after it. Applying the same argument at each step proves monotone descent through the entire cycle. \(\square\)

This result is a descent guarantee, not a global optimality guarantee. The joint objective is nonconvex: a factor update is a least-squares problem, but the full objective depends on products across factor matrices. Monotone objective values therefore do not rule out convergence to a stationary point that is not a global minimizer. Nor does a small change in the objective prove that the factors are uniquely or stably determined.

Diagnosing Slow Convergence and Unstable Fits

Three issues deserve particular attention in practice:

Scaling indeterminacy. Factor columns can become very large in some modes and very small in others without changing the fitted component. This can produce poor numerical scaling. Normalizing columns and recording their magnitudes separately can improve interpretability, but it does not resolve every instability of the fitted decomposition.
Slow convergence. If the conditional design matrix is nearly rank deficient, its Gram matrix \(Z^TZ\) is ill-conditioned. The least-squares update can then be sensitive to small numerical errors. Even when updates are computed accurately, successive factor changes may become small while the objective decreases only slowly. Monitor both the residual objective and the factor or fitted-tensor changes, and use a stated stopping rule rather than treating a few small improvements as proof of optimality.
Local minima and initialization. Different initial factor matrices can lead to different final objective values. Multiple starts, followed by comparison of their attained objectives, can reveal sensitivity to initialization. They do not prove that the best observed solution is globally optimal.

A common mistake is to assume that because every ALS block update is an exact least-squares solution, the final CP fit must be the best rank-\(R\) approximation. The blockwise calculation only optimizes one factor at a time with the others fixed; it does not minimize over all factor matrices jointly in one step. Another mistake is to interpret the scale of an individual factor column as an intrinsic component strength. The rescaling example shows why that quantity depends on the chosen representation. Report the fitted tensor or residual as well as any factor-based interpretation.

Key takeaway: CP-ALS turns a joint nonlinear fit into a sequence of linear least-squares updates. The mode-\(n\) design is the Khatri–Rao product of the other factors, and the update solves \(A^{(n)}(Z^{(n)})^T\approx\mathcal{X}_{(n)}\). Exact updates ensure nonincreasing residual error, but scaling indeterminacy, ill-conditioning, and nonconvexity still require careful diagnosis.

Check Your Understanding

Use the definitions and update derivation above to answer the following.

For mode \(n\), what are the dimensions of \(\mathcal{X}_{(n)}\), \(A^{(n)}\), and \(Z^{(n)}\), and how do they make the unfolded approximation conformable?
State the full-column-rank least-squares update and the condition that permits the inverse in that formula.
Why is \(X_n(Z^\dagger)^T\), rather than \(X_nZ^\dagger\), dimensionally appropriate as the minimum-norm solution when \(Z\) is rank deficient?
Explain why exact ALS updates imply nonincreasing objective values but do not imply global optimality.
Give a nonzero rescaling of two factor columns that leaves a three-way CP component unchanged, and state the constraint the rescaling factors must satisfy.