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 rank and shape
R = 3
tensor_shape = (6, 8, 10)

# Set seed for reproducibility
np.random.seed(0)

# Create a low-rank dense tensor from a Kruskal tensor
factor_matrices = [np.random.rand(s, R) for s in tensor_shape]
M_true = ttb.ktensor(factor_matrices)
X = M_true.to_tensor()

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(1)  # Set seed for reproducibility
few_iters = 10  # Cut off solve early for demo
M1, M1_init, M1_info = ttb.cp_als(X, R, maxiters=few_iters)
CP_ALS:
 Iter 0: f = 8.426960e-01 f-delta = 8.4e-01
 Iter 1: f = 9.014378e-01 f-delta = 5.9e-02
 Iter 2: f = 9.128262e-01 f-delta = 1.1e-02
 Iter 3: f = 9.209368e-01 f-delta = 8.1e-03
 Iter 4: f = 9.286166e-01 f-delta = 7.7e-03
 Iter 5: f = 9.365676e-01 f-delta = 8.0e-03
 Iter 6: f = 9.444309e-01 f-delta = 7.9e-03
 Iter 7: f = 9.514539e-01 f-delta = 7.0e-03
 Iter 8: f = 9.570292e-01 f-delta = 5.6e-03
 Iter 9: f = 9.610205e-01 f-delta = 4.0e-03
 Final f = 9.610205e-01

The cp_als method returns a the following objects:

  1. M1: the solution as a ktensor.

  2. M1_init: the initial guess as a ktensor that was generated at runtime since no initial guess was provided.

  3. M1_info: a dictionary containing runtime information with keys:

    • params: parameters used by cp_als

    • iters: number of iterations performed

    • normresidual: the norm of the residual X.norm()**2 + M.norm()**2 - 2*<X,M>

    • fit: the fit f described above

print("M1_info:")
for k, v in M1_info.items():
    print(f"\t{k}: {v}")
M1_info:
	params: {'stoptol': 0.0001, 'maxiters': 10, 'dimorder': array([0, 1, 2]), 'optdims': array([0, 1, 2]), 'printitn': 1, 'fixsigns': True}
	iters: 9
	normresidual: 0.3792127669517397
	fit: 0.9610204980859988

Run again with a different initial guess, output the initial guess

# Compute a solution with final ktensor stored in M2
np.random.seed(2)  # Set seed for reproducibility
M2, M2_init, _ = ttb.cp_als(
    X, R, maxiters=few_iters
)  # Will not use third output, so leaving unassigned
CP_ALS:
 Iter 0: f = 9.232027e-01 f-delta = 9.2e-01
 Iter 1: f = 9.415585e-01 f-delta = 1.8e-02
 Iter 2: f = 9.460143e-01 f-delta = 4.5e-03
 Iter 3: f = 9.485320e-01 f-delta = 2.5e-03
 Iter 4: f = 9.502871e-01 f-delta = 1.8e-03
 Iter 5: f = 9.516351e-01 f-delta = 1.3e-03
 Iter 6: f = 9.527532e-01 f-delta = 1.1e-03
 Iter 7: f = 9.537623e-01 f-delta = 1.0e-03
 Iter 8: f = 9.547515e-01 f-delta = 9.9e-04
 Iter 9: f = 9.557875e-01 f-delta = 1.0e-03
 Final f = 9.557875e-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.

more_iters = 10 * few_iters
M2_better, _, _ = ttb.cp_als(X, R, maxiters=more_iters, init=M2_init)
CP_ALS:
 Iter 0: f = 9.232027e-01 f-delta = 9.2e-01
 Iter 1: f = 9.415585e-01 f-delta = 1.8e-02
 Iter 2: f = 9.460143e-01 f-delta = 4.5e-03
 Iter 3: f = 9.485320e-01 f-delta = 2.5e-03
 Iter 4: f = 9.502871e-01 f-delta = 1.8e-03
 Iter 5: f = 9.516351e-01 f-delta = 1.3e-03
 Iter 6: f = 9.527532e-01 f-delta = 1.1e-03
 Iter 7: f = 9.537623e-01 f-delta = 1.0e-03
 Iter 8: f = 9.547515e-01 f-delta = 9.9e-04
 Iter 9: f = 9.557875e-01 f-delta = 1.0e-03
 Iter 10: f = 9.569224e-01 f-delta = 1.1e-03
 Iter 11: f = 9.581995e-01 f-delta = 1.3e-03
 Iter 12: f = 9.596531e-01 f-delta = 1.5e-03
 Iter 13: f = 9.613014e-01 f-delta = 1.6e-03
 Iter 14: f = 9.631317e-01 f-delta = 1.8e-03
 Iter 15: f = 9.650854e-01 f-delta = 2.0e-03
 Iter 16: f = 9.670592e-01 f-delta = 2.0e-03
 Iter 17: f = 9.689361e-01 f-delta = 1.9e-03
 Iter 18: f = 9.706308e-01 f-delta = 1.7e-03
 Iter 19: f = 9.721134e-01 f-delta = 1.5e-03
 Iter 20: f = 9.733984e-01 f-delta = 1.3e-03
 Iter 21: f = 9.745196e-01 f-delta = 1.1e-03
 Iter 22: f = 9.755127e-01 f-delta = 9.9e-04
 Iter 23: f = 9.764073e-01 f-delta = 8.9e-04
 Iter 24: f = 9.772262e-01 f-delta = 8.2e-04
 Iter 25: f = 9.779860e-01 f-delta = 7.6e-04
 Iter 26: f = 9.786985e-01 f-delta = 7.1e-04
 Iter 27: f = 9.793724e-01 f-delta = 6.7e-04
 Iter 28: f = 9.800139e-01 f-delta = 6.4e-04
 Iter 29: f = 9.806278e-01 f-delta = 6.1e-04
 Iter 30: f = 9.812174e-01 f-delta = 5.9e-04
 Iter 31: f = 9.817854e-01 f-delta = 5.7e-04
 Iter 32: f = 9.823337e-01 f-delta = 5.5e-04
 Iter 33: f = 9.828638e-01 f-delta = 5.3e-04
 Iter 34: f = 9.833771e-01 f-delta = 5.1e-04
 Iter 35: f = 9.838744e-01 f-delta = 5.0e-04
 Iter 36: f = 9.843567e-01 f-delta = 4.8e-04
 Iter 37: f = 9.848247e-01 f-delta = 4.7e-04
 Iter 38: f = 9.852789e-01 f-delta = 4.5e-04
 Iter 39: f = 9.857198e-01 f-delta = 4.4e-04
 Iter 40: f = 9.861479e-01 f-delta = 4.3e-04
 Iter 41: f = 9.865637e-01 f-delta = 4.2e-04
 Iter 42: f = 9.869675e-01 f-delta = 4.0e-04
 Iter 43: f = 9.873597e-01 f-delta = 3.9e-04
 Iter 44: f = 9.877407e-01 f-delta = 3.8e-04
 Iter 45: f = 9.881107e-01 f-delta = 3.7e-04
 Iter 46: f = 9.884700e-01 f-delta = 3.6e-04
 Iter 47: f = 9.888190e-01 f-delta = 3.5e-04
 Iter 48: f = 9.891579e-01 f-delta = 3.4e-04
 Iter 49: f = 9.894870e-01 f-delta = 3.3e-04
 Iter 50: f = 9.898065e-01 f-delta = 3.2e-04
 Iter 51: f = 9.901168e-01 f-delta = 3.1e-04
 Iter 52: f = 9.904180e-01 f-delta = 3.0e-04
 Iter 53: f = 9.907104e-01 f-delta = 2.9e-04
 Iter 54: f = 9.909942e-01 f-delta = 2.8e-04
 Iter 55: f = 9.912696e-01 f-delta = 2.8e-04
 Iter 56: f = 9.915369e-01 f-delta = 2.7e-04
 Iter 57: f = 9.917963e-01 f-delta = 2.6e-04
 Iter 58: f = 9.920480e-01 f-delta = 2.5e-04
 Iter 59: f = 9.922922e-01 f-delta = 2.4e-04
 Iter 60: f = 9.925290e-01 f-delta = 2.4e-04
 Iter 61: f = 9.927588e-01 f-delta = 2.3e-04
 Iter 62: f = 9.929816e-01 f-delta = 2.2e-04
 Iter 63: f = 9.931977e-01 f-delta = 2.2e-04
 Iter 64: f = 9.934073e-01 f-delta = 2.1e-04
 Iter 65: f = 9.936105e-01 f-delta = 2.0e-04
 Iter 66: f = 9.938075e-01 f-delta = 2.0e-04
 Iter 67: f = 9.939984e-01 f-delta = 1.9e-04
 Iter 68: f = 9.941835e-01 f-delta = 1.9e-04
 Iter 69: f = 9.943629e-01 f-delta = 1.8e-04
 Iter 70: f = 9.945368e-01 f-delta = 1.7e-04
 Iter 71: f = 9.947053e-01 f-delta = 1.7e-04
 Iter 72: f = 9.948686e-01 f-delta = 1.6e-04
 Iter 73: f = 9.950268e-01 f-delta = 1.6e-04
 Iter 74: f = 9.951801e-01 f-delta = 1.5e-04
 Iter 75: f = 9.953287e-01 f-delta = 1.5e-04
 Iter 76: f = 9.954725e-01 f-delta = 1.4e-04
 Iter 77: f = 9.956119e-01 f-delta = 1.4e-04
 Iter 78: f = 9.957469e-01 f-delta = 1.4e-04
 Iter 79: f = 9.958777e-01 f-delta = 1.3e-04
 Iter 80: f = 9.960044e-01 f-delta = 1.3e-04
 Iter 81: f = 9.961270e-01 f-delta = 1.2e-04
 Iter 82: f = 9.962458e-01 f-delta = 1.2e-04
 Iter 83: f = 9.963609e-01 f-delta = 1.2e-04
 Iter 84: f = 9.964723e-01 f-delta = 1.1e-04
 Iter 85: f = 9.965802e-01 f-delta = 1.1e-04
 Iter 86: f = 9.966847e-01 f-delta = 1.0e-04
 Iter 87: f = 9.967859e-01 f-delta = 1.0e-04
 Iter 88: f = 9.968838e-01 f-delta = 9.8e-05
 Final f = 9.968838e-01

Compare the two solutions

Use the ktensor score() member function to compare the two solutions. A score of 1 indicates a perfect match.

score_M2 = M2.score(M_true)
score_M2_better = M2_better.score(M_true)

Here, score() returns a tuple with the score as the first element:

print(f"Score of M2 and M_true: {score_M2[0]}")
print(f"Score of M2_better and M_true: {score_M2_better[0]}")
Score of M2 and M_true: 0.5167723537387411
Score of M2_better and M_true: 0.9519267496183947

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.

M2_rerun, _, _ = ttb.cp_als(X, R, maxiters=few_iters, init=M2_init)
score_M2_rerun = M2.score(M2_rerun)  # Score of 1 indicates the same solution
print(f"Score: {score_M2_rerun[0]}")
CP_ALS:
 Iter 0: f = 9.232027e-01 f-delta = 9.2e-01
 Iter 1: f = 9.415585e-01 f-delta = 1.8e-02
 Iter 2: f = 9.460143e-01 f-delta = 4.5e-03
 Iter 3: f = 9.485320e-01 f-delta = 2.5e-03
 Iter 4: f = 9.502871e-01 f-delta = 1.8e-03
 Iter 5: f = 9.516351e-01 f-delta = 1.3e-03
 Iter 6: f = 9.527532e-01 f-delta = 1.1e-03
 Iter 7: f = 9.537623e-01 f-delta = 1.0e-03
 Iter 8: f = 9.547515e-01 f-delta = 9.9e-04
 Iter 9: f = 9.557875e-01 f-delta = 1.0e-03
 Final f = 9.557875e-01
Score: 1.0000000000000004

Changing the output frequency

Using the printitn option to change the output frequency.

M = ttb.cp_als(X, R, maxiters=few_iters, printitn=2)
CP_ALS:
 Iter 0: f = 8.805093e-01 f-delta = 8.8e-01
 Iter 2: f = 9.379232e-01 f-delta = 1.5e-02
 Iter 4: f = 9.445392e-01 f-delta = 3.0e-03
 Iter 6: f = 9.492378e-01 f-delta = 2.2e-03
 Iter 8: f = 9.532334e-01 f-delta = 2.0e-03
 Final f = 9.552078e-01

Suppress all output

Set printitn to zero to suppress all output.

M = ttb.cp_als(X, R, maxiters=few_iters, printitn=0)  # No output

Use HOSVD initial guess

Use the "nvecs" option to use the leading mode-\(n\) singular vectors as the initial guess. The initialization process will require more computation in general, but it may lead to better solutions.

M3, _, _ = ttb.cp_als(X, R, maxiters=few_iters, init="nvecs", printitn=0)
score_M3 = M3.score(M_true)
print(f"Score of M2 and M_true: {score_M2[0]}")
print(f"Score of M3 and M_true: {score_M3[0]}")
Score of M2 and M_true: 0.5167723537387411
Score of M3 and M_true: 0.5406177000260893

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.

M4, _, M4_info = ttb.cp_als(
    X, R, maxiters=few_iters, init="nvecs", printitn=0, dimorder=[2, 1, 0]
)
score_M4 = M4.score(M_true)
print(f"Score of M3 and M_true: {score_M3[0]}")
print(f"Score of M4 and M_true: {score_M4[0]}")
Score of M3 and M_true: 0.5406177000260893
Score of M4 and M_true: 0.6314317590547762

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.

M4_rerun, _, M4_rerun_info = ttb.cp_als(X, R, init="nvecs", **M4_info["params"])
score_M4_rerun = M4.score(M4_rerun)
print(f"Score of M4 and M4_rerun: {score_M4_rerun[0]}")

print("M4_info:")
for k, v in M4_info.items():
    print(f"\t{k}: {v}")

print("M4_rerun_info:")
for k, v in M4_rerun_info.items():
    print(f"\t{k}: {v}")
Score of M4 and M4_rerun: 0.999999999999997
M4_info:
	params: {'stoptol': 0.0001, 'maxiters': 10, 'dimorder': array([2, 1, 0]), 'optdims': array([0, 1, 2]), 'printitn': 0, 'fixsigns': True}
	iters: 9
	normresidual: 0.45173622073803993
	fit: 0.9535657698910674
M4_rerun_info:
	params: {'stoptol': 0.0001, 'maxiters': 10, 'dimorder': array([2, 1, 0]), 'optdims': array([0, 1, 2]), 'printitn': 0, 'fixsigns': True}
	iters: 9
	normresidual: 0.45173622073800845
	fit: 0.9535657698910707

Change the stopping tolerance

It’s also possible to loosen or tighten the stopping tolerance on the change in the fit. Note that you may need to increase the number of iterations for it to converge.

M5 = ttb.cp_als(X, R, init="nvecs", maxiters=1000, stoptol=1e-12, printitn=100)
CP_ALS:
 Iter 0: f = 8.214438e-01 f-delta = 8.2e-01
 Iter 100: f = 9.983903e-01 f-delta = 5.1e-05
 Iter 200: f = 9.999079e-01 f-delta = 2.5e-06
 Iter 300: f = 9.999932e-01 f-delta = 1.8e-07
 Iter 400: f = 9.999995e-01 f-delta = 1.4e-08
 Iter 467: f = 9.999999e-01 f-delta = 0.0e+00
 Final f = 9.999999e-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.

# 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]]),
    ],
    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) with order F
weights=[1. 1.]
factor_matrices[0] =
[[  1.   1.]
 [ -1. -10.]]
factor_matrices[1] =
[[  1.   1.]
 [ -2. -10.]]

M_fixsigns=
ktensor of shape (2, 2) with order F
weights=[101.           3.16227766]
factor_matrices[0] =
[[-0.09950372 -0.70710678]
 [ 0.99503719  0.70710678]]
factor_matrices[1] =
[[-0.09950372 -0.4472136 ]
 [ 0.99503719  0.89442719]]

M_no_fixsigns=
ktensor of shape (2, 2) with order F
weights=[101.           3.16227766]
factor_matrices[0] =
[[ 0.09950372  0.70710678]
 [-0.99503719 -0.70710678]]
factor_matrices[1] =
[[ 0.09950372  0.4472136 ]
 [-0.99503719 -0.89442719]]

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.