Sum of Structured Tensors

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.

When certain operations are performed on a tensor which is formed as a sum of tensors, it can be beneficial to avoid explicitly forming the sum. For example, if a tensor is formed as a sum of a low rank tensor and a sparse tensor, the structure of the summands can make storage, decomposition and operations with other tensors significantly more efficient. A sumtensor exploits this structure. Here we explain the basics of defining and using sumtensors.

import pyttb as ttb
import numpy as np
import pickle
def estimate_mem_usage(variable) -> int:
    """
    Python variables contain references to memory.
    Quickly estimate memory usage of custom types.
    """
    return len(pickle.dumps(variable))

Creating sumtensors

A sumtensor T can only be declared as a sum of same-shaped tensors T1, T2,…,TN. The summand tensors are stored internally, which define the “parts” of the sumtensor. The parts of a sumtensor can be (dense) tensors (tensor), sparse tensors ( sptensor), Kruskal tensors (ktensor), or Tucker tensors (ttensor). An example of the use of the sumtensor constructor follows.

T1 = ttb.tenones((3, 3, 3))
T2 = ttb.sptensor(
    subs=np.array([[0, 0, 0], [1, 1, 1], [2, 2, 1], [1, 0, 0]]),
    vals=np.ones((4, 1)),
    shape=(3, 3, 3),
)
T = ttb.sumtensor([T1, T2])
print(T)
sumtensor of shape (3, 3, 3) with 2 parts:
Part 0: 
	tensor of shape (3, 3, 3) with order F
	data[:, :, 0] =
	[[1. 1. 1.]
	 [1. 1. 1.]
	 [1. 1. 1.]]
	data[:, :, 1] =
	[[1. 1. 1.]
	 [1. 1. 1.]
	 [1. 1. 1.]]
	data[:, :, 2] =
	[[1. 1. 1.]
	 [1. 1. 1.]
	 [1. 1. 1.]]
Part 1: 
	sparse tensor of shape (3, 3, 3) with 4 nonzeros and order F
	[0, 0, 0] = 1.0
	[1, 1, 1] = 1.0
	[2, 2, 1] = 1.0
	[1, 0, 0] = 1.0

A Magnitude Example

For large-scale problems, the sumtensor class may make the difference as to whether or not a tensor can be stored in memory. Consider the following example, where \(\mathcal{T}\) is of size \(1000 x 1000 x 1000\), formed from the sum of a ktensor and an sptensor.

np.random.seed(0)
X1 = np.random.rand(50, 3)
X2 = np.random.rand(50, 3)
X3 = np.random.rand(50, 3)
K = ttb.ktensor([X1, X2, X3], np.ones((3,)), copy=False)
S = ttb.sptenrand((50, 50, 50), 1e-100)

ST = ttb.sumtensor([K, S])
TT = ST.full()
print(
    f"Size of sumtensor: {estimate_mem_usage(ST)}\n"
    f"Size of tensor: {estimate_mem_usage(TT)}"
)
WARNING:root:Selected no copy, but input factor matrices aren't F ordered so must copy.
Size of sumtensor: 4200
Size of tensor: 1000224

Further examples of the sumtensor constructor

We can declare an empty sumtensor, with no parts.

P = ttb.sumtensor()
print(P)
Empty sumtensor

sumtensor also supports a copy constructor.

S = P.copy()
print(S)
Empty sumtensor

Ndims and shape for dimensions of a sumtensor

For a given sumtensor, ndims returns the number of modes and shape returns the shape in each dimension of the sumtensor.

print(f"Ndims: {T.ndims}\n" f"Shape: {T.shape}")
Ndims: 3
Shape: (3, 3, 3)

Use full to convert sumtensor to a dense tensor

The full method can convert all the parts of a sumtensor to a dense tensor. Note that for large tensors, this can use a large amount of memory to expand then sum the parts.

print(T.full())
tensor of shape (3, 3, 3) with order F
data[:, :, 0] =
[[2. 1. 1.]
 [2. 1. 1.]
 [1. 1. 1.]]
data[:, :, 1] =
[[1. 1. 1.]
 [1. 2. 1.]
 [1. 1. 2.]]
data[:, :, 2] =
[[1. 1. 1.]
 [1. 1. 1.]
 [1. 1. 1.]]

Use double to convert to a numpy array

The double method can convert the parts of a sumtensor to a dense numpy array. Similar warnings for memory usages as full.

print(T.double())
[[[2. 1. 1.]
  [1. 1. 1.]
  [1. 1. 1.]]

 [[2. 1. 1.]
  [1. 2. 1.]
  [1. 1. 1.]]

 [[1. 1. 1.]
  [1. 1. 1.]
  [1. 2. 1.]]]

Matricized Khatri-Rao product of a sumtensor

The mttkrp method computes the Khatri-Rao product of a matricized tensor and sumtensor. The required arguments are:

  • A list of matrices (or a ktensor)

  • A mode n

The list of matrices must consist of m matrices, where m is the number of modes in the sumtensor. The number of columns in all matrices should be the same and the number of rows of matrix i should match the dimension of the sumtensor shape in mode i. For more details see the documentation of tensor.mttkrp.

matrices = [np.eye(3), np.ones((3, 3)), np.random.rand(3, 3)]
n = 1
T.mttkrp(matrices, n)
array([[2.23269848, 1.57660651, 1.94578466],
       [1.28174587, 2.07389535, 1.94578466],
       [1.28174587, 1.34318625, 2.82750487]])

Innerproducts of sumtensors

The innerprod method computes the inner product of a sumtensors parts with other tensor types.

S = ttb.sptensor(
    subs=np.array([[0, 0, 0], [1, 1, 1], [2, 2, 1], [1, 0, 0]]),
    vals=np.ones((4, 1)),
    shape=(3, 3, 3),
)
T.innerprod(S)
8.0

Norm compatibility interface

The norm method just returns 0 and issues a warning. Norm cannot be distributed, but some algorithms access the norm for verbose details.

T.norm()
/home/docs/checkouts/readthedocs.org/user_builds/pyttb/envs/374/lib/python3.9/site-packages/pyttb/sumtensor.py:451: UserWarning: Sumtensor doesn't actually support norm. Returning 0 for compatibility.
  warnings.warn(
0.0

Use CP-ALS with sumtensor

One of the primary motivations for defining the sumtensor class is for efficient decomposition. In particular, when trying to find a CP decomposition of a tensor using alternating least squares, the subproblems can be efficiently created and solved using mttkrp and innerprod. Both of these operations can be performed more efficiently by exploiting extra structure in the tensors which form the sum, so the performance of cp_als is also improved. Consider the following example, where a cp_als is run on a sumtensor.

result, _, _ = ttb.cp_als(T, 2, maxiters=10)
print(result)
CP_ALS:
 Iter 0: f = -3.565172e+01 f-delta = 3.6e+01
 Iter 1: f = -3.697690e+01 f-delta = 1.3e+00
 Iter 2: f = -3.709311e+01 f-delta = 1.2e-01
 Iter 3: f = -3.719637e+01 f-delta = 1.0e-01
 Iter 4: f = -3.729319e+01 f-delta = 9.7e-02
 Iter 5: f = -3.737617e+01 f-delta = 8.3e-02
 Iter 6: f = -3.744308e+01 f-delta = 6.7e-02
 Iter 7: f = -3.749770e+01 f-delta = 5.5e-02
 Iter 8: f = -3.754461e+01 f-delta = 4.7e-02
 Iter 9: f = -3.758560e+01 f-delta = 4.1e-02
 Final f = -3.758560e+01
ktensor of shape (3, 3, 3) with order F
weights=[4.91985657 2.34303787]
factor_matrices[0] =
[[0.61330126 0.34390829]
 [0.63326328 0.53070383]
 [0.47205846 0.77464865]]
factor_matrices[1] =
[[ 0.8094263  -0.09808623]
 [ 0.46429877  0.56423944]
 [ 0.3595215   0.81976396]]
factor_matrices[2] =
[[0.73488876 0.19632928]
 [0.48660539 0.85888068]
 [0.47239148 0.47305263]]

It follows that in cases where \(\mathcal{T}\) is too large for its full expansion to be stored in memory, we may still be able find a CP decomposition by exploiting the sumtensor structure.

Note that the fit returned by cp_als is not correct for sumtensor, because the norm operation is not supported.

Addition with sumtensors

Sumtensors can be added to any other type of tensor. The result is a new sumtensor with the tensor appended to the parts of the original sumtensor. Note that the tensor is always appended, despite the order of the operation.

# Equivalent additions despite the order
print(f"T+S:\n{T+S}\n")
print(f"S+T:\n{S+T}")
T+S:
sumtensor of shape (3, 3, 3) with 3 parts:
Part 0: 
	tensor of shape (3, 3, 3) with order F
	data[:, :, 0] =
	[[1. 1. 1.]
	 [1. 1. 1.]
	 [1. 1. 1.]]
	data[:, :, 1] =
	[[1. 1. 1.]
	 [1. 1. 1.]
	 [1. 1. 1.]]
	data[:, :, 2] =
	[[1. 1. 1.]
	 [1. 1. 1.]
	 [1. 1. 1.]]
Part 1: 
	sparse tensor of shape (3, 3, 3) with 4 nonzeros and order F
	[0, 0, 0] = 1.0
	[1, 1, 1] = 1.0
	[2, 2, 1] = 1.0
	[1, 0, 0] = 1.0
Part 2: 
	sparse tensor of shape (3, 3, 3) with 4 nonzeros and order F
	[0, 0, 0] = 1.0
	[1, 1, 1] = 1.0
	[2, 2, 1] = 1.0
	[1, 0, 0] = 1.0

S+T:
sumtensor of shape (3, 3, 3) with 3 parts:
Part 0: 
	tensor of shape (3, 3, 3) with order F
	data[:, :, 0] =
	[[1. 1. 1.]
	 [1. 1. 1.]
	 [1. 1. 1.]]
	data[:, :, 1] =
	[[1. 1. 1.]
	 [1. 1. 1.]
	 [1. 1. 1.]]
	data[:, :, 2] =
	[[1. 1. 1.]
	 [1. 1. 1.]
	 [1. 1. 1.]]
Part 1: 
	sparse tensor of shape (3, 3, 3) with 4 nonzeros and order F
	[0, 0, 0] = 1.0
	[1, 1, 1] = 1.0
	[2, 2, 1] = 1.0
	[1, 0, 0] = 1.0
Part 2: 
	sparse tensor of shape (3, 3, 3) with 4 nonzeros and order F
	[0, 0, 0] = 1.0
	[1, 1, 1] = 1.0
	[2, 2, 1] = 1.0
	[1, 0, 0] = 1.0

Accessing sumtensor parts

Subscripted reference can be used to access individual parts of the sumtensor.

print(f"Part 0:\n{T.parts[0]}\n\n" f"Part 1:\n{T.parts[1]}")
Part 0:
tensor of shape (3, 3, 3) with order F
data[:, :, 0] =
[[1. 1. 1.]
 [1. 1. 1.]
 [1. 1. 1.]]
data[:, :, 1] =
[[1. 1. 1.]
 [1. 1. 1.]
 [1. 1. 1.]]
data[:, :, 2] =
[[1. 1. 1.]
 [1. 1. 1.]
 [1. 1. 1.]]

Part 1:
sparse tensor of shape (3, 3, 3) with 4 nonzeros and order F
[0, 0, 0] = 1.0
[1, 1, 1] = 1.0
[2, 2, 1] = 1.0
[1, 0, 0] = 1.0