Revision

**85fbc1e20d309b4d1d31353bc8ed0f74bfd37050**authored by Meraj Hashemizadeh on**30 August 2021, 20:41:50 UTC**, committed by Meraj Hashemizadeh on**30 August 2021, 20:41:50 UTC****1 parent**5a6992a

cp_tensor.py

```
"""
Core operations on CP tensors.
"""
from . import backend as T
from .base import fold, tensor_to_vec
from ._factorized_tensor import FactorizedTensor
from .tenalg import khatri_rao, multi_mode_dot
from .utils import DefineDeprecated
import numpy as np
# Author: Jean Kossaifi
# License: BSD 3 clause
class CPTensor(FactorizedTensor):
def __init__(self, cp_tensor):
super().__init__()
shape, rank = _validate_cp_tensor(cp_tensor)
weights, factors = cp_tensor
# Should we allow None weights?
if weights is None:
weights = T.ones(rank, **T.context(factors[0]))
self.shape = shape
self.rank = rank
self.factors = factors
self.weights = weights
def __getitem__(self, index):
if index == 0:
return self.weights
elif index == 1:
return self.factors
else:
raise IndexError('You tried to access index {} of a CP tensor.\n'
'You can only access index 0 and 1 of a CP tensor'
'(corresponding respectively to the weights and factors)'.format(index))
def __setitem__(self, index, value):
if index == 0:
self.weights = value
elif index == 1:
self.factors = value
else:
raise IndexError('You tried to set the value at index {} of a CP tensor.\n'
'You can only set index 0 and 1 of a CP tensor'
'(corresponding respectively to the weights and factors)'.format(index))
def __iter__(self):
yield self.weights
yield self.factors
def __len__(self):
return 2
def __repr__(self):
message = '(weights, factors) : rank-{} CPTensor of shape {} '.format(self.rank, self.shape)
return message
def to_tensor(self):
return cp_to_tensor(self)
def to_vec(self):
return cp_to_vec(self)
def to_unfolded(self, mode):
return cp_to_unfolded(self, mode)
def mode_dot(self, matrix_or_vector, mode, keep_dim=False, copy=True):
"""n-mode product of a CP tensor and a matrix or vector at the specified mode
Parameters
----------
cp_tensor : tl.CPTensor or (core, factors)
matrix_or_vector : ndarray
1D or 2D array of shape ``(J, i_k)`` or ``(i_k, )``
matrix or vectors to which to n-mode multiply the tensor
mode : int
Returns
-------
CPTensor = (core, factors)
`mode`-mode product of `tensor` by `matrix_or_vector`
* of shape :math:`(i_1, ..., i_{k-1}, J, i_{k+1}, ..., i_N)` if matrix_or_vector is a matrix
* of shape :math:`(i_1, ..., i_{k-1}, i_{k+1}, ..., i_N)` if matrix_or_vector is a vector
See also
--------
kruskal_multi_mode_dot : chaining several mode_dot in one call
"""
return cp_mode_dot(self, matrix_or_vector, mode, keep_dim=keep_dim, copy=copy)
def norm(self):
"""Returns the l2 norm of a CP tensor
Parameters
----------
cp_tensor : tl.CPTensor or (core, factors)
Returns
-------
l2-norm : int
Notes
-----
This is ||cp_to_tensor(factors)||^2
You can see this using the fact that
khatria-rao(A, B)^T x khatri-rao(A, B) = A^T x A * B^T x B
"""
return cp_norm(self)
def normalize(self, inplace=True):
"""Normalizes the factors to unit length
Turns ``factors = [|U_1, ... U_n|]`` into ``[weights; |V_1, ... V_n|]``,
where the columns of each `V_k` are normalized to unit Euclidean length
from the columns of `U_k` with the normalizing constants absorbed into
`weights`. In the special case of a symmetric tensor, `weights` holds the
eigenvalues of the tensor.
Parameters
----------
cp_tensor : CPTensor = (weight, factors)
factors is list of matrices, all with the same number of columns
i.e.::
for u in U:
u[i].shape == (s_i, R)
where `R` is fixed while `s_i` can vary with `i`
inplace : bool, default is True
if False, returns a normalized Copy
otherwise the tensor modifies itself and returns itself
Returns
-------
CPTensor = (normalisation_weights, normalised_factors)
returns itself if inplace is False, a normalized copy otherwise
"""
self.weights, self.factors = cp_normalize(self)
def _validate_cp_tensor(cp_tensor):
"""Validates a cp_tensor in the form (weights, factors)
Returns the rank and shape of the validated tensor
Parameters
----------
cp_tensor : CPTensor or (weights, factors)
Returns
-------
(shape, rank) : (int tuple, int)
size of the full tensor and rank of the CP tensor
"""
if isinstance(cp_tensor, CPTensor):
# it's already been validated at creation
return cp_tensor.shape, cp_tensor.rank
elif isinstance(cp_tensor, (float, int)): #0-order tensor
return 0, 0
weights, factors = cp_tensor
# if len(factors) < 2:
# raise ValueError('A CP tensor should be composed of at least two factors.'
# 'However, {} factor was given.'.format(len(factors)))
if T.ndim(factors[0]) == 2:
rank = int(T.shape(factors[0])[1])
else:
rank = 1
shape = []
for i, factor in enumerate(factors):
s = T.shape(factor)
if len(s) == 2:
current_mode_size, current_rank = s
else:
current_mode_size, current_rank = s, 1
if current_rank != rank:
raise ValueError('All the factors of a CP tensor should have the same number of column.'
'However, factors[0].shape[1]={} but factors[{}].shape[1]={}.'.format(
rank, i, T.shape(factor)[1]))
shape.append(current_mode_size)
if weights is not None and T.shape(weights) != (rank, ):
raise ValueError('Given factors for a rank-{} CP tensor but len(weights)={}.'.format(
rank, T.shape(weights)))
return tuple(shape), rank
def _cp_n_param(tensor_shape, rank, weights=False):
"""Number of parameters of a CP decomposition for a given `rank` and full `tensor_shape`.
Parameters
----------
tensor_shape : int tuple
shape of the full tensor to decompose (or approximate)
rank : tuple
rank of the CP decomposition
Returns
-------
n_params : int
Number of parameters of a CP decomposition of rank `rank` of a full tensor of shape `tensor_shape`
"""
factors_params = rank*np.sum(tensor_shape)
if weights:
return factors_params + rank
else:
return factors_params
def validate_cp_rank(tensor_shape, rank='same', rounding='round'):
"""Returns the rank of a CP Decomposition
Parameters
----------
tensor_shape : tupe
shape of the tensor to decompose
rank : {'same', float, int}, default is same
way to determine the rank, by default 'same'
if 'same': rank is computed to keep the number of parameters (at most) the same
if float, computes a rank so as to keep rank percent of the original number of parameters
if int, just returns rank
rounding = {'round', 'floor', 'ceil'}
Returns
-------
rank : int
rank of the decomposition
"""
if rounding == 'ceil':
rounding_fun = np.ceil
elif rounding == 'floor':
rounding_fun = np.floor
elif rounding == 'round':
rounding_fun = np.round
else:
raise ValueError(f'Rounding should be of round, floor or ceil, but got {rounding}')
if rank == 'same':
rank = float(1)
if isinstance(rank, float):
rank = int(rounding_fun(np.prod(tensor_shape)*rank/np.sum(tensor_shape)))
return rank
def cp_normalize(cp_tensor):
"""Returns cp_tensor with factors normalised to unit length
Turns ``factors = [|U_1, ... U_n|]`` into ``[weights; |V_1, ... V_n|]``,
where the columns of each `V_k` are normalized to unit Euclidean length
from the columns of `U_k` with the normalizing constants absorbed into
`weights`. In the special case of a symmetric tensor, `weights` holds the
eigenvalues of the tensor.
Parameters
----------
cp_tensor : CPTensor = (weight, factors)
factors is list of matrices, all with the same number of columns
i.e.::
for u in U:
u[i].shape == (s_i, R)
where `R` is fixed while `s_i` can vary with `i`
Returns
-------
CPTensor = (normalisation_weights, normalised_factors)
"""
_, rank = _validate_cp_tensor(cp_tensor)
weights, factors = cp_tensor
if weights is None:
weights = T.ones(rank, **T.context(factors[0]))
normalized_factors = []
for i, factor in enumerate(factors):
if i == 0:
factor = factor*weights
weights = T.ones(rank, **T.context(factor))
scales = T.norm(factor, axis=0)
scales_non_zero = T.where(scales==0, T.ones(T.shape(scales), **T.context(factor)), scales)
weights = weights*scales
normalized_factors.append(factor / T.reshape(scales_non_zero, (1, -1)))
return CPTensor((weights, normalized_factors))
def cp_flip_sign(cp_tensor, mode=0, func=None):
"""Returns cp_tensor with factors flipped to have positive signs.
The sign of a given column is determined by `func`, which is the mean
by default. Any negative signs are assigned to the mode indicated by `mode`.
Parameters
----------
cp_tensor : CPTensor = (weight, factors)
factors is list of matrices, all with the same number of columns
i.e.::
for u in U:
u[i].shape == (s_i, R)
where `R` is fixed while `s_i` can vary with `i`
mode: int
mode that should receive negative signs
func: tensorly function
a function that should summarize the sign of a column
it must be able to take an axis argument
Returns
-------
CPTensor = (normalisation_weights, normalised_factors)
"""
_validate_cp_tensor(cp_tensor)
weights, factors = cp_tensor
if func is None:
func = T.mean
for jj in range(0, len(factors)):
# Skip the target mode
if jj == mode:
continue
# Calculate the sign of the current factor in each component
column_signs = T.sign(func(factors[jj], axis=0))
# Update both the current and receiving factor
factors[mode] = factors[mode]*column_signs[np.newaxis, :]
factors[jj] = factors[jj]*column_signs[np.newaxis, :]
# Check the weight signs
weight_signs = T.sign(weights)
factors[mode] = factors[mode]*weight_signs[np.newaxis, :]
weights = T.abs(weights)
return CPTensor((weights, factors))
def cp_lstsq_grad(cp_tensor, tensor, return_loss=False, mask=None):
"""This function computes (for a third-order tensor)
.. math::
\nabla 0.5 ||\\mathcal{X} - [\\mathbf{w}; \\mathbf{A}, \\mathbf{B}, \\mathbf{C}]||^2
where :math:`[\\mathbf{w}; \\mathbf{A}, \\mathbf{B}, \\mathbf{C}]` is the CP decomposition with weights
:math:`\\mathbf{w}` and factor matrices :math:`\\mathbf{A}`, :math:`\\mathbf{B}` and :math:`\\mathbf{C}`.
Note that this does not return the gradient with respect to the weights even if CP is normalized.
Parameters
----------
cp_tensor : CPTensor = (weight, factors)
factors is a list of factor matrices, all with the same number of columns
i.e. for all matrix U in factor_matrices:
U has shape ``(s_i, R)``, where R is fixed and s_i varies with i
mask : ndarray
A mask to be applied to the final tensor. It should be
broadcastable to the shape of the final tensor, that is
``(U[1].shape[0], ... U[-1].shape[0])``.
return_loss : bool
Optionally return the scalar loss function along with the gradient.
Returns
-------
cp_gradient : CPTensor = (None, factors)
factors is a list of factor matrix gradients, all with the same number of columns
i.e. for all matrix U in factor_matrices:
U has shape ``(s_i, R)``, where R is fixed and s_i varies with i
loss : float
Scalar quantity of the loss function corresponding to cp_gradient. Only returned
if return_loss = True.
"""
_validate_cp_tensor(cp_tensor)
_, factors = cp_tensor
diff = tensor - cp_to_tensor(cp_tensor)
if mask is not None:
diff = diff * mask
grad_fac = [-unfolding_dot_khatri_rao(diff, cp_tensor, ii) for ii in range(len(factors))]
if return_loss:
return CPTensor((None, grad_fac)), 0.5*T.sum(diff**2)
return CPTensor((None, grad_fac))
def cp_to_tensor(cp_tensor, mask=None):
"""Turns the Khatri-product of matrices into a full tensor
``factor_matrices = [|U_1, ... U_n|]`` becomes
a tensor shape ``(U[1].shape[0], U[2].shape[0], ... U[-1].shape[0])``
Parameters
----------
cp_tensor : CPTensor = (weight, factors)
factors is a list of factor matrices, all with the same number of columns
i.e. for all matrix U in factor_matrices:
U has shape ``(s_i, R)``, where R is fixed and s_i varies with i
mask : ndarray a mask to be applied to the final tensor. It should be
broadcastable to the shape of the final tensor, that is
``(U[1].shape[0], ... U[-1].shape[0])``.
Returns
-------
ndarray
full tensor of shape ``(U[1].shape[0], ... U[-1].shape[0])``
Notes
-----
This version works by first computing the mode-0 unfolding of the tensor
and then refolding it.
There are other possible and equivalent alternate implementation, e.g.
summing over r and updating an outer product of vectors.
"""
shape, _ = _validate_cp_tensor(cp_tensor)
if not shape: # 0-order tensor
return cp_tensor
weights, factors = cp_tensor
if len(shape) == 1: # just a vector
return T.sum(weights*factors[0], axis=1)
if weights is None:
weights = 1
if mask is None:
full_tensor = T.dot(factors[0]*weights,
T.transpose(khatri_rao(factors, skip_matrix=0)))
else:
full_tensor = T.sum(khatri_rao([factors[0]*weights]+factors[1:], mask=mask), axis=1)
return fold(full_tensor, 0, shape)
def cp_to_unfolded(cp_tensor, mode):
"""Turns the khatri-product of matrices into an unfolded tensor
turns ``factors = [|U_1, ... U_n|]`` into a mode-`mode`
unfolding of the tensor
Parameters
----------
cp_tensor : CPTensor = (weight, factors)
factors is a list of matrices, all with the same number of columns
ie for all u in factor_matrices:
u[i] has shape (s_u_i, R), where R is fixed
mode: int
mode of the desired unfolding
Returns
-------
ndarray
unfolded tensor of shape (tensor_shape[mode], -1)
Notes
-----
Writing factors = [U_1, ..., U_n], we exploit the fact that
``U_k = U[k].dot(khatri_rao(U_1, ..., U_k-1, U_k+1, ..., U_n))``
"""
_validate_cp_tensor(cp_tensor)
weights, factors = cp_tensor
if weights is not None:
return T.dot(factors[mode]*weights, T.transpose(khatri_rao(factors, skip_matrix=mode)))
else:
return T.dot(factors[mode], T.transpose(khatri_rao(factors, skip_matrix=mode)))
def cp_to_vec(cp_tensor):
"""Turns the khatri-product of matrices into a vector
(the tensor ``factors = [|U_1, ... U_n|]``
is converted into a raveled mode-0 unfolding)
Parameters
----------
cp_tensor : CPTensor = (weight, factors)
factors is a list of matrices, all with the same number of columns
i.e.::
for u in U:
u[i].shape == (s_i, R)
where `R` is fixed while `s_i` can vary with `i`
Returns
-------
ndarray
vectorised tensor
"""
return tensor_to_vec(cp_to_tensor(cp_tensor))
def cp_mode_dot(cp_tensor, matrix_or_vector, mode, keep_dim=False, copy=False):
"""n-mode product of a CP tensor and a matrix or vector at the specified mode
Parameters
----------
cp_tensor : tl.CPTensor or (core, factors)
matrix_or_vector : ndarray
1D or 2D array of shape ``(J, i_k)`` or ``(i_k, )``
matrix or vectors to which to n-mode multiply the tensor
mode : int
Returns
-------
CPTensor = (core, factors)
`mode`-mode product of `tensor` by `matrix_or_vector`
* of shape :math:`(i_1, ..., i_{k-1}, J, i_{k+1}, ..., i_N)` if matrix_or_vector is a matrix
* of shape :math:`(i_1, ..., i_{k-1}, i_{k+1}, ..., i_N)` if matrix_or_vector is a vector
See also
--------
kruskal_multi_mode_dot : chaining several mode_dot in one call
"""
shape, _ = _validate_cp_tensor(cp_tensor)
weights, factors = cp_tensor
contract = False
if T.ndim(matrix_or_vector) == 2: # Tensor times matrix
# Test for the validity of the operation
if matrix_or_vector.shape[1] != shape[mode]:
raise ValueError(
'shapes {0} and {1} not aligned in mode-{2} multiplication: {3} (mode {2}) != {4} (dim 1 of matrix)'.format(
shape, matrix_or_vector.shape, mode, shape[mode], matrix_or_vector.shape[1]
))
elif T.ndim(matrix_or_vector) == 1: # Tensor times vector
if matrix_or_vector.shape[0] != shape[mode]:
raise ValueError(
'shapes {0} and {1} not aligned for mode-{2} multiplication: {3} (mode {2}) != {4} (vector size)'.format(
shape, matrix_or_vector.shape, mode, shape[mode], matrix_or_vector.shape[0]
))
if not keep_dim:
contract = True # Contract over that mode
else:
raise ValueError('Can only take n_mode_product with a vector or a matrix.')
if copy:
factors = [T.copy(f) for f in factors]
weights = T.copy(weights)
if contract:
factor = factors.pop(mode)
factor = T.dot(matrix_or_vector, factor)
mode = max(mode - 1, 0)
factors[mode] *= factor
else:
factors[mode] = T.dot(matrix_or_vector, factors[mode])
if copy:
return CPTensor((weights, factors))
else:
cp_tensor.shape = tuple(f.shape[0] for f in factors)
return cp_tensor
def unfolding_dot_khatri_rao(tensor, cp_tensor, mode):
"""mode-n unfolding times khatri-rao product of factors
Parameters
----------
tensor : tl.tensor
tensor to unfold
factors : tl.tensor list
list of matrices of which to the khatri-rao product
mode : int
mode on which to unfold `tensor`
Returns
-------
mttkrp
dot(unfold(tensor, mode), khatri-rao(factors))
Notes
-----
This is a variant of::
unfolded = unfold(tensor, mode)
kr_factors = khatri_rao(factors, skip_matrix=mode)
mttkrp2 = tl.dot(unfolded, kr_factors)
Multiplying with the Khatri-Rao product is equivalent to multiplying,
for each rank, with the kronecker product of each factor.
In code::
mttkrp_parts = []
for r in range(rank):
component = tl.tenalg.multi_mode_dot(tensor, [f[:, r] for f in factors], skip=mode)
mttkrp_parts.append(component)
mttkrp = tl.stack(mttkrp_parts, axis=1)
return mttkrp
This can be done by taking n-mode-product with the full factors
(faster but more memory consuming)::
projected = multi_mode_dot(tensor, factors, skip=mode, transpose=True)
ndims = T.ndim(tensor)
res = []
for i in range(factors[0].shape[1]):
index = tuple([slice(None) if k == mode else i for k in range(ndims)])
res.append(projected[index])
return T.stack(res, axis=-1)
The same idea could be expressed using einsum::
ndims = tl.ndim(tensor)
tensor_idx = ''.join(chr(ord('a') + i) for i in range(ndims))
rank = chr(ord('a') + ndims + 1)
op = tensor_idx
for i in range(ndims):
if i != mode:
op += ',' + ''.join([tensor_idx[i], rank])
else:
result = ''.join([tensor_idx[i], rank])
op += '->' + result
factors = [f for (i, f) in enumerate(factors) if i != mode]
return tl_einsum(op, tensor, *factors)
"""
mttkrp_parts = []
_, rank = _validate_cp_tensor(cp_tensor)
weights, factors = cp_tensor
for r in range(rank):
component = multi_mode_dot(tensor, [f[:, r] for f in factors], skip=mode)
mttkrp_parts.append(component)
if weights is None:
return T.stack(mttkrp_parts, axis=1)
else:
return T.stack(mttkrp_parts, axis=1)*T.reshape(weights, (1, -1))
def cp_norm(cp_tensor):
"""Returns the l2 norm of a CP tensor
Parameters
----------
cp_tensor : tl.CPTensor or (core, factors)
Returns
-------
l2-norm : int
Notes
-----
This is ||cp_to_tensor(factors)||^2
You can see this using the fact that
khatria-rao(A, B)^T x khatri-rao(A, B) = A^T x A * B^T x B
"""
_ = _validate_cp_tensor(cp_tensor)
weights, factors = cp_tensor
norm = 1
for factor in factors:
norm *= T.dot(T.transpose(factor), factor)
if weights is not None:
#norm = T.dot(T.dot(weights, norm), weights)
norm = norm * (T.reshape(weights, (-1, 1))*T.reshape(weights, (1, -1)))
# We sum even if weigths is not None
# as e.g. MXNet would return a 1D tensor, not a 0D tensor
return T.sqrt(T.sum(norm))
# Deprecated classes and functions
KruskalTensor = DefineDeprecated(deprecated_name='KruskalTensor', use_instead=CPTensor)
kruskal_norm = DefineDeprecated(deprecated_name='kruskal_norm', use_instead=cp_norm)
kruskal_mode_dot = DefineDeprecated(deprecated_name='kruskal_mode_dot', use_instead=cp_mode_dot)
kruskal_to_tensor = DefineDeprecated(deprecated_name='kruskal_to_tensor', use_instead=cp_to_tensor)
kruskal_to_unfolded = DefineDeprecated(deprecated_name='kruskal_to_unfolded', use_instead=cp_to_unfolded)
kruskal_to_vec = DefineDeprecated(deprecated_name='kruskal_to_vec', use_instead=cp_to_vec)
```

Computing file changes ...