Alternating least squares for Canonical Polyadic (CP) Decomposition
Copyright 2025 National Technology & Engineering Solutions of Sandia,
LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the
U.S. Government retains certain rights in this software.
The function cp_als computes an estimate of the best rank-\(R\) CP model of a tensor \(\mathcal{X}\) using the well-known alternating least-squares algorithm (see, e.g., Kolda and Bader, SIAM Review, 2009, for more information). The input \(\mathcal{X}\) can be almost any type of tensor including a tensor, sptensor, ktensor, or ttensor. The output CP model is a ktensor.
from __future__ import annotations
import os
import sys
import numpy as np
import pyttb as ttb
Generate low-rank data tensor
# Choose the rank and shape
R = 2
tensor_shape = (3, 4, 5)
# Set the random seed for reproducibility
np.random.seed(0)
# Create a low-rank dense tensor from a Kruskal tensor (i.e., ktensor)
factor_matrices = [np.random.randn(s, R) for s in tensor_shape]
M_true = ttb.ktensor(factor_matrices) # true solution
X = M_true.to_tensor()
Run cp_als using default parameters
The cp_als method has two required inputs: a data tensor (X) and the rank of the CP model (R) to compute.
By default, cp_als uses a random initial guess. At each iteration, it reports the fit f which is defined for a data tensor X and CP model M as
f = 1 - ( X.norm()**2 + M.norm()**2 - 2*<X,M> ) / X.norm()
and is loosely the proportion of the data described by the CP model, i.e., a fit of 1 is perfect.
# Compute a solution with final ktensor stored in M1
np.random.seed(1) # Set the random seed for reproducibility
M1, M1_init, M1_info = ttb.cp_als(X, R)
CP_ALS:
Iter 0: f = 5.629369e-01 f-delta = 5.6e-01
Iter 1: f = 6.994515e-01 f-delta = 1.4e-01
Iter 2: f = 7.186123e-01 f-delta = 1.9e-02
Iter 3: f = 7.270663e-01 f-delta = 8.5e-03
Iter 4: f = 7.305063e-01 f-delta = 3.4e-03
Iter 5: f = 7.317567e-01 f-delta = 1.3e-03
Iter 6: f = 7.321614e-01 f-delta = 4.0e-04
Iter 7: f = 7.323097e-01 f-delta = 1.5e-04
Iter 8: f = 7.323887e-01 f-delta = 7.9e-05
Final f = 7.323887e-01
The cp_als method returns the following:
M1: the solution as aktensor.M1_init: the initial guess as aktensorthat was used in computing the solution.M1_info: adictcontaining runtime information with keys:params: parameters used bycp_alsiters: number of iterations performednormresidual: the norm of the residualX.norm()**2 + M.norm()**2 - 2*<X,M>fit: the fitfdescribed above
print("M1_info:")
for k, v in M1_info.items():
print(f"\t{k}: {v}")
M1_info:
params: {'stoptol': 0.0001, 'maxiters': 1000, 'dimorder': array([0, 1, 2]), 'optdims': array([0, 1, 2]), 'printitn': 1, 'fixsigns': True}
iters: 8
normresidual: 3.1416394034929502
fit: 0.7323887077713218
Run cp_als using different initial guesses
Different random initial guesses can be generated and used by setting different random seeds (via numpy). You can also explicitly set the init parameter to "random" to make it clear that the default random initialization is being used, although this is not necessary as illustrated above.
np.random.seed(2) # Set seed for reproducibility
M2, M2_init, _ = ttb.cp_als(X, R, init="random") # leaving third output unassigned
CP_ALS:
Iter 0: f = 7.483199e-01 f-delta = 7.5e-01
Iter 1: f = 8.417330e-01 f-delta = 9.3e-02
Iter 2: f = 8.622039e-01 f-delta = 2.0e-02
Iter 3: f = 8.851075e-01 f-delta = 2.3e-02
Iter 4: f = 9.126355e-01 f-delta = 2.8e-02
Iter 5: f = 9.418493e-01 f-delta = 2.9e-02
Iter 6: f = 9.666331e-01 f-delta = 2.5e-02
Iter 7: f = 9.830955e-01 f-delta = 1.6e-02
Iter 8: f = 9.920977e-01 f-delta = 9.0e-03
Iter 9: f = 9.964589e-01 f-delta = 4.4e-03
Iter 10: f = 9.984440e-01 f-delta = 2.0e-03
Iter 11: f = 9.993222e-01 f-delta = 8.8e-04
Iter 12: f = 9.997058e-01 f-delta = 3.8e-04
Iter 13: f = 9.998725e-01 f-delta = 1.7e-04
Iter 14: f = 9.999448e-01 f-delta = 7.2e-05
Final f = 9.999448e-01
A specific ktensor can also be used as an initial guess. As an example, using the same initial guess (and all other parameters) as the previous run of cp_als gives the exact same solution.
M2_rerun, _, _ = ttb.cp_als(X, R, init=M2_init)
CP_ALS:
Iter 0: f = 7.483199e-01 f-delta = 7.5e-01
Iter 1: f = 8.417330e-01 f-delta = 9.3e-02
Iter 2: f = 8.622039e-01 f-delta = 2.0e-02
Iter 3: f = 8.851075e-01 f-delta = 2.3e-02
Iter 4: f = 9.126355e-01 f-delta = 2.8e-02
Iter 5: f = 9.418493e-01 f-delta = 2.9e-02
Iter 6: f = 9.666331e-01 f-delta = 2.5e-02
Iter 7: f = 9.830955e-01 f-delta = 1.6e-02
Iter 8: f = 9.920977e-01 f-delta = 9.0e-03
Iter 9: f = 9.964589e-01 f-delta = 4.4e-03
Iter 10: f = 9.984440e-01 f-delta = 2.0e-03
Iter 11: f = 9.993222e-01 f-delta = 8.8e-04
Iter 12: f = 9.997058e-01 f-delta = 3.8e-04
Iter 13: f = 9.998725e-01 f-delta = 1.7e-04
Iter 14: f = 9.999448e-01 f-delta = 7.2e-05
Final f = 9.999448e-01
The "nvecs" initialization option computes the initial guess using the eigenvectors of the outer product of the matricized data tensor. This initialization process will require more computation than the default random initialization, but it can often lead to better solutions in fewer iterations.
M2_nvecs, _, _ = ttb.cp_als(X, R, init="nvecs")
CP_ALS:
Iter 0: f = 9.399785e-01 f-delta = 9.4e-01
Iter 1: f = 9.735771e-01 f-delta = 3.4e-02
Iter 2: f = 9.889950e-01 f-delta = 1.5e-02
Iter 3: f = 9.955336e-01 f-delta = 6.5e-03
Iter 4: f = 9.981764e-01 f-delta = 2.6e-03
Iter 5: f = 9.992526e-01 f-delta = 1.1e-03
Iter 6: f = 9.996931e-01 f-delta = 4.4e-04
Iter 7: f = 9.998738e-01 f-delta = 1.8e-04
Iter 8: f = 9.999480e-01 f-delta = 7.4e-05
Final f = 9.999480e-01
Evaluate and compare the outputs of cp_als
Use the ktensor.score() method to compare outputs of cp_als to the true solution if known. In the examples above, the data tensor was generated from a ktensor, M_true, which can be used to evaluate solutions computed using cp_als. A score of 1 indicates a perfect match.
Note that the ktensor.score() method returns a tuple with the score as the first element and other information related to the score computation as the remaining elements. See the ktensor documentation for more information about the return values.
score_M1 = M1.score(M_true) # not a good solution
score_M2 = M2.score(M_true) # a better solution
score_M2_nvecs = M2_nvecs.score(M_true) # an even better solution
print(f"Score of M1 and M_true: {score_M1[0]}")
print(f"Score of M2 and M_true: {score_M2[0]}")
print(f"Score of M2_nvecs and M_true: {score_M2_nvecs[0]}")
Score of M1 and M_true: 0.16018319811919496
Score of M2 and M_true: 0.9999905577227326
Score of M2_nvecs and M_true: 0.999995876091137
When two solutions match, as is the case with M2 and M2_rerun, the score will be 1 (up to floating point error).
score_M2_rerun = M2.score(M2_rerun)
print(f"Score of M2 and M2_rerun: {score_M2_rerun[0]}")
Score of M2 and M2_rerun: 1.0000000000000007
Increase the maximum number of iterations
Note that the previous run kicked out at only 10 iterations, before reaching the specified convegence tolerance. Let’s increase the maximum number of iterations and try again, using the same initial guess.
more_iters = 100
M2_better, _, _ = ttb.cp_als(X, R, maxiters=more_iters, init=M_true)
CP_ALS:
Iter 0: f = 1.000000e+00 f-delta = 1.0e+00
Iter 1: f = 1.000000e+00 f-delta = 2.9e-08
Final f = 1.000000e+00
Changing the output frequency
Using the printitn option to change the output frequency.
M = ttb.cp_als(X, R, printitn=5)
CP_ALS:
Iter 0: f = 7.423397e-01 f-delta = 7.4e-01
Iter 5: f = 9.948635e-01 f-delta = 7.5e-03
Iter 10: f = 9.999370e-01 f-delta = 8.8e-05
Final f = 9.999370e-01
Suppress all output
Set printitn to zero to suppress all output.
M = ttb.cp_als(X, R, printitn=0) # No output
Use initial guess based
Use the "nvecs" option to initialize cp_als with the leading mode-\(n\) singular vectors of the input tensor. This initialization process will require more computation than the default "random initialization in general, but it can lead to better solutions.
Change the order of the dimensions in CP
Changing the order of the dimensions in which cp_als iterates over the input tensor can lead to a different solution.
M3, _, M3_info = ttb.cp_als(X, R, init="nvecs", printitn=0, dimorder=[2, 1, 0])
score_M3 = M3.score(M_true)
print(f"Score of M3 and M_true: {score_M2[0]}")
print(f"Score of M4 and M_true: {score_M3[0]}")
Score of M3 and M_true: 0.9999905577227326
Score of M4 and M_true: 0.9999943868161334
In the last example, we also collected the third output argument M4_info which has runtime information in it. The field M4_info["iters"] has the total number of iterations. The field M4_info["params"] has the information used to run the method. Unless the initialization method is "random", passing the parameters back to the method will yield the exact same results.
M3_rerun, _, M3_rerun_info = ttb.cp_als(X, R, init="nvecs", **M3_info["params"])
score_M3_rerun = M3.score(M3_rerun)
print(f"Score of M3 and M3_rerun: {score_M3_rerun[0]}")
print("M3_info:")
for k, v in M3_info.items():
print(f"\t{k}: {v}")
print("M3_rerun_info:")
for k, v in M3_rerun_info.items():
print(f"\t{k}: {v}")
Score of M3 and M3_rerun: 1.0
M3_info:
params: {'stoptol': 0.0001, 'maxiters': 1000, 'dimorder': array([2, 1, 0]), 'optdims': array([0, 1, 2]), 'printitn': 0, 'fixsigns': True}
iters: 8
normresidual: 0.000416248318775757
fit: 0.9999645431139068
M3_rerun_info:
params: {'stoptol': 0.0001, 'maxiters': 1000, 'dimorder': array([2, 1, 0]), 'optdims': array([0, 1, 2]), 'printitn': 0, 'fixsigns': True}
iters: 8
normresidual: 0.000416248318775757
fit: 0.9999645431139068
Change the stopping tolerance
It’s also possible to loosen or tighten the stopping tolerance on the change in the fit, using the stoptol parameter. Note that you may need to increase the number of iterations for it to converge. Note that you may need to increase the maimum number of iterations allowed, maxiters, in order to compute a decomposition whose fit is below the requested stopping tolerance.
M4 = ttb.cp_als(X, R, init="nvecs", maxiters=1000, stoptol=1e-12, printitn=100)
CP_ALS:
Iter 0: f = 9.399785e-01 f-delta = 9.4e-01
Iter 18: f = 1.000000e+00 f-delta = 0.0e+00
Final f = 1.000000e+00
Control sign ambiguity of factor matrices
The default behavior of cp_als is to make a call to fixsigns() to fix the sign ambiguity of the factor matrices. You can turn off this behavior by passing the fixsigns parameter value of False when calling cp_als.
# Create rank-2 tensor
X2 = ttb.ktensor(
factor_matrices=[
np.array([[1.0, 1.0], [-1.0, -10.0]]),
np.array([[1.0, 1.0], [-2.0, -10.0]]),
np.array([[1.0, 1.0], [-3.0, -10.0]]),
],
weights=np.array([1.0, 1.0]),
)
print(f"X2=\n{X2}\n")
M_fixsigns, _, _ = ttb.cp_als(X2, 2, printitn=0, init=ttb.ktensor(X2.factor_matrices))
print(f"M_fixsigns=\n{M_fixsigns}\n") # default behavior, fixsigns called
M_no_fixsigns, _, _ = ttb.cp_als(
X2, 2, printitn=0, init=ttb.ktensor(X2.factor_matrices), fixsigns=False
)
print(f"M_no_fixsigns=\n{M_no_fixsigns}\n") # fixsigns not called
X2=
ktensor of shape (2, 2, 2) with order F
weights=[1. 1.]
factor_matrices[0] =
[[ 1. 1.]
[ -1. -10.]]
factor_matrices[1] =
[[ 1. 1.]
[ -2. -10.]]
factor_matrices[2] =
[[ 1. 1.]
[ -3. -10.]]
M_fixsigns=
ktensor of shape (2, 2, 2) with order F
weights=[1015.03743773 10. ]
factor_matrices[0] =
[[-0.09950372 0.70710678]
[ 0.99503719 -0.70710678]]
factor_matrices[1] =
[[-0.09950372 -0.4472136 ]
[ 0.99503719 0.89442719]]
factor_matrices[2] =
[[ 0.09950372 -0.31622777]
[-0.99503719 0.9486833 ]]
M_no_fixsigns=
ktensor of shape (2, 2, 2) with order F
weights=[1015.03743773 10. ]
factor_matrices[0] =
[[ 0.09950372 0.70710678]
[-0.99503719 -0.70710678]]
factor_matrices[1] =
[[ 0.09950372 0.4472136 ]
[-0.99503719 -0.89442719]]
factor_matrices[2] =
[[ 0.09950372 0.31622777]
[-0.99503719 -0.9486833 ]]
Recommendations
Run multiple times with different guesses and select the solution with the best fit.
Try different ranks and choose the solution that is the best descriptor for your data based on the combination of the fit and the interpretation of the factors, e.g., by visualizing the results.