Alternating least squares for Canonical Polyadic (CP) Decomposition
Copyright 2022 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.
import os
import sys
import pyttb as ttb
import numpy as np
Generate data
# Pick the shape and rank
R = 3
np.random.seed(0) # Set seed for reproducibility
X = ttb.tenrand(shape=(6, 8, 10))
Basic call to the method, specifying the data tensor and its rank
This uses a random initial guess. At each iteration, it reports the fit f which is defined 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(0) # Set seed for reproducibility
short_tutorial = 10 # Cut off solve early for demo
M1 = ttb.cp_als(X, R, maxiters=short_tutorial)
CP_ALS:
Iter 0: f = 5.251960e-01 f-delta = 5.3e-01
Iter 1: f = 5.332788e-01 f-delta = 8.1e-03
Iter 2: f = 5.356720e-01 f-delta = 2.4e-03
Iter 3: f = 5.378495e-01 f-delta = 2.2e-03
Iter 4: f = 5.396460e-01 f-delta = 1.8e-03
Iter 5: f = 5.408062e-01 f-delta = 1.2e-03
Iter 6: f = 5.415884e-01 f-delta = 7.8e-04
Iter 7: f = 5.422437e-01 f-delta = 6.6e-04
Iter 8: f = 5.428984e-01 f-delta = 6.5e-04
Iter 9: f = 5.436212e-01 f-delta = 7.2e-04
Final f = 5.436212e-01
Since we set only a single output, M1 is actually a tuple containing:
M1[0]: the solution as aktensor.M1[1]: the initial guess as aktensorthat was generated at runtime since no initial guess was provided.M1[2]: a dictionary containing 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(f"M1[2]['params']: {M1[2]['params']}")
print(f"M1[2]['iters']: {M1[2]['iters']}")
print(f"M1[2]['normresidual']: {M1[2]['normresidual']}")
print(f"M1[2]['fit']: {M1[2]['fit']}")
M1[2]['params']: {'stoptol': 0.0001, 'maxiters': 10, 'dimorder': array([0, 1, 2]), 'optdims': array([0, 1, 2]), 'printitn': 1, 'fixsigns': True}
M1[2]['iters']: 9
M1[2]['normresidual']: 5.781169892859293
M1[2]['fit']: 0.5436212193118705
Run again with a different initial guess, output the initial guess.
np.random.seed(1) # Set seed for reproducibility
M2bad, Minit, _ = ttb.cp_als(X, R, maxiters=short_tutorial)
CP_ALS:
Iter 0: f = 5.171905e-01 f-delta = 5.2e-01
Iter 1: f = 5.310255e-01 f-delta = 1.4e-02
Iter 2: f = 5.354620e-01 f-delta = 4.4e-03
Iter 3: f = 5.378540e-01 f-delta = 2.4e-03
Iter 4: f = 5.389924e-01 f-delta = 1.1e-03
Iter 5: f = 5.396324e-01 f-delta = 6.4e-04
Iter 6: f = 5.400953e-01 f-delta = 4.6e-04
Iter 7: f = 5.404937e-01 f-delta = 4.0e-04
Iter 8: f = 5.408788e-01 f-delta = 3.9e-04
Iter 9: f = 5.412842e-01 f-delta = 4.1e-04
Final f = 5.412842e-01
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.
less_short_tutorial = 10 * short_tutorial
M2 = ttb.cp_als(X, R, maxiters=less_short_tutorial, init=Minit)
CP_ALS:
Iter 0: f = 5.171905e-01 f-delta = 5.2e-01
Iter 1: f = 5.310255e-01 f-delta = 1.4e-02
Iter 2: f = 5.354620e-01 f-delta = 4.4e-03
Iter 3: f = 5.378540e-01 f-delta = 2.4e-03
Iter 4: f = 5.389924e-01 f-delta = 1.1e-03
Iter 5: f = 5.396324e-01 f-delta = 6.4e-04
Iter 6: f = 5.400953e-01 f-delta = 4.6e-04
Iter 7: f = 5.404937e-01 f-delta = 4.0e-04
Iter 8: f = 5.408788e-01 f-delta = 3.9e-04
Iter 9: f = 5.412842e-01 f-delta = 4.1e-04
Iter 10: f = 5.417389e-01 f-delta = 4.5e-04
Iter 11: f = 5.422686e-01 f-delta = 5.3e-04
Iter 12: f = 5.428873e-01 f-delta = 6.2e-04
Iter 13: f = 5.435793e-01 f-delta = 6.9e-04
Iter 14: f = 5.442831e-01 f-delta = 7.0e-04
Iter 15: f = 5.449112e-01 f-delta = 6.3e-04
Iter 16: f = 5.454076e-01 f-delta = 5.0e-04
Iter 17: f = 5.457786e-01 f-delta = 3.7e-04
Iter 18: f = 5.460630e-01 f-delta = 2.8e-04
Iter 19: f = 5.462960e-01 f-delta = 2.3e-04
Iter 20: f = 5.464972e-01 f-delta = 2.0e-04
Iter 21: f = 5.466759e-01 f-delta = 1.8e-04
Iter 22: f = 5.468362e-01 f-delta = 1.6e-04
Iter 23: f = 5.469800e-01 f-delta = 1.4e-04
Iter 24: f = 5.471087e-01 f-delta = 1.3e-04
Iter 25: f = 5.472236e-01 f-delta = 1.1e-04
Iter 26: f = 5.473257e-01 f-delta = 1.0e-04
Iter 27: f = 5.474162e-01 f-delta = 9.0e-05
Final f = 5.474162e-01
Compare the two solutions
Use the ktensor score() member function to compare the two solutions. A score of 1 indicates a perfect match.
M1_ktns = M1[0]
M2_ktns = M2[0]
score = M1_ktns.score(M2_ktns)
Here, score() returned a tuple score with the score as the first element:
score[0]
np.float64(0.5476368270211815)
See the ktensor documentation for more information about the return values of score().
Rerun with same initial guess
Using the same initial guess (and all other parameters) gives the exact same solution.
M2alt = ttb.cp_als(X, R, maxiters=less_short_tutorial, init=Minit)
M2alt_ktns = M2alt[0]
score = M2_ktns.score(M2alt_ktns) # Score of 1 indicates the same solution
print(f"Score: {score[0]}.")
CP_ALS:
Iter 0: f = 5.171905e-01 f-delta = 5.2e-01
Iter 1: f = 5.310255e-01 f-delta = 1.4e-02
Iter 2: f = 5.354620e-01 f-delta = 4.4e-03
Iter 3: f = 5.378540e-01 f-delta = 2.4e-03
Iter 4: f = 5.389924e-01 f-delta = 1.1e-03
Iter 5: f = 5.396324e-01 f-delta = 6.4e-04
Iter 6: f = 5.400953e-01 f-delta = 4.6e-04
Iter 7: f = 5.404937e-01 f-delta = 4.0e-04
Iter 8: f = 5.408788e-01 f-delta = 3.9e-04
Iter 9: f = 5.412842e-01 f-delta = 4.1e-04
Iter 10: f = 5.417389e-01 f-delta = 4.5e-04
Iter 11: f = 5.422686e-01 f-delta = 5.3e-04
Iter 12: f = 5.428873e-01 f-delta = 6.2e-04
Iter 13: f = 5.435793e-01 f-delta = 6.9e-04
Iter 14: f = 5.442831e-01 f-delta = 7.0e-04
Iter 15: f = 5.449112e-01 f-delta = 6.3e-04
Iter 16: f = 5.454076e-01 f-delta = 5.0e-04
Iter 17: f = 5.457786e-01 f-delta = 3.7e-04
Iter 18: f = 5.460630e-01 f-delta = 2.8e-04
Iter 19: f = 5.462960e-01 f-delta = 2.3e-04
Iter 20: f = 5.464972e-01 f-delta = 2.0e-04
Iter 21: f = 5.466759e-01 f-delta = 1.8e-04
Iter 22: f = 5.468362e-01 f-delta = 1.6e-04
Iter 23: f = 5.469800e-01 f-delta = 1.4e-04
Iter 24: f = 5.471087e-01 f-delta = 1.3e-04
Iter 25: f = 5.472236e-01 f-delta = 1.1e-04
Iter 26: f = 5.473257e-01 f-delta = 1.0e-04
Iter 27: f = 5.474162e-01 f-delta = 9.0e-05
Final f = 5.474162e-01
Score: 1.0000000000000007.
Changing the output frequency
Using the printitn option to change the output frequency.
M2alt2 = ttb.cp_als(X, R, maxiters=less_short_tutorial, init=Minit, printitn=20)
CP_ALS:
Iter 0: f = 5.171905e-01 f-delta = 5.2e-01
Iter 20: f = 5.464972e-01 f-delta = 2.0e-04
Iter 27: f = 5.474162e-01 f-delta = 9.0e-05
Final f = 5.474162e-01
Suppress all output
Set printitn to zero to suppress all output.
M2alt2 = ttb.cp_als(X, R, printitn=0) # No output
Use HOSVD initial guess
Use the "nvecs" option to use the leading mode-\(n\) singular vectors as the initial guess.
M3 = ttb.cp_als(X, R, init="nvecs", printitn=20)
s = M2[0].score(M3[0])
print(f"score(M2,M3) = {s[0]}")
CP_ALS:
Iter 0: f = 5.404929e-01 f-delta = 5.4e-01
Iter 14: f = 5.544095e-01 f-delta = 9.3e-05
Final f = 5.544095e-01
score(M2,M3) = 0.26711137527148554
Change the order of the dimensions in CP
M4, _, info = ttb.cp_als(X, 3, dimorder=[1, 2, 0], init="nvecs", printitn=20)
s = M2[0].score(M4)
print(f"score(M2,M4) = {s[0]}")
CP_ALS:
Iter 0: f = 5.359701e-01 f-delta = 5.4e-01
Iter 16: f = 5.536798e-01 f-delta = 1.0e-04
Final f = 5.536798e-01
score(M2,M4) = 0.32904706039891574
In the last example, we also collected the third output argument info which has runtime information in it. The field info["iters"] has the total number of iterations. The field 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.
M4alt, _, info = ttb.cp_als(X, 3, **info["params"])
s = M4alt.score(M4)
print(f"score(M4alt,M4) = {s[0]}")
CP_ALS:
Iter 0: f = 5.277447e-01 f-delta = 5.3e-01
Iter 20: f = 5.445464e-01 f-delta = 2.5e-04
Iter 26: f = 5.454859e-01 f-delta = 9.1e-05
Final f = 5.454859e-01
score(M4alt,M4) = 0.3501120666752804
Change the tolerance
It’s also possible to loosen or tighten the tolerance on the change in the fit. You may need to increase the number of iterations for it to converge.
M5 = ttb.cp_als(X, 3, init="nvecs", stoptol=1e-12, printitn=100)
CP_ALS:
Iter 0: f = 5.404929e-01 f-delta = 5.4e-01
Iter 100: f = 5.551901e-01 f-delta = 9.7e-07
Iter 200: f = 5.552383e-01 f-delta = 2.4e-07
Iter 300: f = 5.552543e-01 f-delta = 1.1e-07
Iter 400: f = 5.552623e-01 f-delta = 6.0e-08
Iter 500: f = 5.552671e-01 f-delta = 3.9e-08
Iter 600: f = 5.552704e-01 f-delta = 2.7e-08
Iter 700: f = 5.552727e-01 f-delta = 2.0e-08
Iter 800: f = 5.552744e-01 f-delta = 1.5e-08
Iter 900: f = 5.552757e-01 f-delta = 1.2e-08
Final f = 5.552768e-01
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.
X = ttb.ktensor(
factor_matrices=[
np.array([[1.0, 1.0], [1.0, -10.0]]),
np.array([[1.0, 1.0], [1.0, -10.0]]),
],
weights=np.array([1.0, 1.0]),
)
M1 = ttb.cp_als(X, 2, printitn=1, init=ttb.ktensor(X.factor_matrices))
print(M1[0]) # default behavior, fixsigns called
M2 = ttb.cp_als(X, 2, printitn=1, init=ttb.ktensor(X.factor_matrices), fixsigns=False)
print(M2[0]) # fixsigns not called
CP_ALS:
Iter 0: f = 1.000000e+00 f-delta = 1.0e+00
Iter 1: f = 1.000000e+00 f-delta = 0.0e+00
Final f = 1.000000e+00
ktensor of shape (2, 2) with order F
weights=[101. 2.]
factor_matrices[0] =
[[-0.09950372 0.70710678]
[ 0.99503719 0.70710678]]
factor_matrices[1] =
[[-0.09950372 0.70710678]
[ 0.99503719 0.70710678]]
CP_ALS:
Iter 0: f = 1.000000e+00 f-delta = 1.0e+00
Iter 1: f = 1.000000e+00 f-delta = 0.0e+00
Final f = 1.000000e+00
ktensor of shape (2, 2) with order F
weights=[101. 2.]
factor_matrices[0] =
[[ 0.09950372 0.70710678]
[-0.99503719 0.70710678]]
factor_matrices[1] =
[[ 0.09950372 0.70710678]
[-0.99503719 0.70710678]]
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.