Source code for tensorly.decomposition._nn_cp

import warnings
import tensorly as tl
from ._base_decomposition import DecompositionMixin
from ._cp import initialize_cp
from ..tenalg.proximal import hals_nnls
from ..cp_tensor import (
    CPTensor,
    unfolding_dot_khatri_rao,
    cp_norm,
    cp_normalize,
    validate_cp_rank,
)
from ..tenalg.svd import svd_interface

# Authors: Jean Kossaifi <jean.kossaifi+tensors@gmail.com>
#          Chris Swierczewski <csw@amazon.com>
#          Sam Schneider <samjohnschneider@gmail.com>
#          Aaron Meurer <asmeurer@gmail.com>
#          Aaron Meyer <tensorly@ameyer.me>
#          Jeremy Cohen <jeremy.cohen@irisa.fr>
#          Axel Marmoret <axel.marmoret@inria.fr>
#          Caglayan TUna <caglayantun@gmail.com>

# License: BSD 3 clause


[docs] def non_negative_parafac( tensor, rank, n_iter_max=100, init="svd", svd="truncated_svd", tol=10e-7, random_state=None, verbose=0, normalize_factors=False, return_errors=False, mask=None, cvg_criterion="abs_rec_error", fixed_modes=None, ): """ Non-negative CP decomposition Uses multiplicative updates, see [2]_ Parameters ---------- tensor : ndarray rank : int number of components n_iter_max : int maximum number of iteration init : {'svd', 'random'}, optional svd : str, default is 'truncated_svd' function to use to compute the SVD, acceptable values in tensorly.SVD_FUNS tol : float, optional tolerance: the algorithm stops when the variation in the reconstruction error is less than the tolerance random_state : {None, int, np.random.RandomState} verbose : int, optional level of verbosity normalize_factors : if True, aggregate the weights of each factor in a 1D-tensor of shape (rank, ), which will contain the norms of the factors fixed_modes : list, default is None A list of modes for which the initial value is not modified. The last mode cannot be fixed due to error computation. Returns ------- factors : ndarray list list of positive factors of the CP decomposition element `i` is of shape ``(tensor.shape[i], rank)`` References ---------- .. [2] Amnon Shashua and Tamir Hazan, "Non-negative tensor factorization with applications to statistics and computer vision", In Proceedings of the International Conference on Machine Learning (ICML), pp 792-799, ICML, 2005 """ epsilon = tl.eps(tensor.dtype) rank = validate_cp_rank(tl.shape(tensor), rank=rank) weights, factors = initialize_cp( tensor, rank, init=init, svd=svd, non_negative=True, mask=mask, random_state=random_state, normalize_factors=normalize_factors, ) rec_errors = [] norm_tensor = tl.norm(tensor, 2) if fixed_modes is None: fixed_modes = [] if tl.ndim(tensor) - 1 in fixed_modes: warnings.warn( "You asked for fixing the last mode, which is not supported while tol is fixed.\n The last mode will not be fixed. Consider using tl.moveaxis()" ) fixed_modes.remove(tl.ndim(tensor) - 1) modes_list = [mode for mode in range(tl.ndim(tensor)) if mode not in fixed_modes] for iteration in range(n_iter_max): if verbose > 1: print("Starting iteration", iteration + 1) for mode in modes_list: if verbose > 1: print("Mode", mode, "of", tl.ndim(tensor)) accum = 1 # khatri_rao(factors).tl.dot(khatri_rao(factors)) # simplifies to multiplications sub_indices = [i for i in range(len(factors)) if i != mode] for i, e in enumerate(sub_indices): if i: accum *= tl.dot(tl.transpose(factors[e]), factors[e]) else: accum = tl.dot(tl.transpose(factors[e]), factors[e]) accum = tl.reshape(weights, (-1, 1)) * accum * tl.reshape(weights, (1, -1)) if mask is not None: tensor = tensor * mask + tl.cp_to_tensor( (weights, factors), mask=1 - mask ) mttkrp = unfolding_dot_khatri_rao(tensor, (weights, factors), mode) numerator = tl.clip(mttkrp, a_min=epsilon, a_max=None) denominator = tl.dot(factors[mode], accum) denominator = tl.clip(denominator, a_min=epsilon, a_max=None) factor = factors[mode] * numerator / denominator factors[mode] = factor if normalize_factors and mode != modes_list[-1]: weights, factors = cp_normalize((weights, factors)) if tol: # ||tensor - rec||^2 = ||tensor||^2 + ||rec||^2 - 2*<tensor, rec> factors_norm = cp_norm((weights, factors)) # mttkrp and factor for the last mode. This is equivalent to the # inner product <tensor, factorization> iprod = tl.sum(tl.sum(mttkrp * factor, axis=0)) rec_error = ( tl.sqrt(tl.abs(norm_tensor**2 + factors_norm**2 - 2 * iprod)) / norm_tensor ) rec_errors.append(rec_error) if iteration >= 1: rec_error_decrease = rec_errors[-2] - rec_errors[-1] if verbose: print( f"iteration {iteration}, reconstraction error: {rec_error}, decrease = {rec_error_decrease}" ) if cvg_criterion == "abs_rec_error": stop_flag = abs(rec_error_decrease) < tol elif cvg_criterion == "rec_error": stop_flag = rec_error_decrease < tol else: raise TypeError("Unknown convergence criterion") if stop_flag: if verbose: print(f"PARAFAC converged after {iteration} iterations") break else: if verbose: print(f"reconstruction error={rec_errors[-1]}") if normalize_factors: weights, factors = cp_normalize((weights, factors)) cp_tensor = CPTensor((weights, factors)) if return_errors: return cp_tensor, rec_errors else: return cp_tensor
[docs] def non_negative_parafac_hals( tensor, rank, n_iter_max=100, init="svd", svd="truncated_svd", tol=10e-8, random_state=None, sparsity_coefficients=None, fixed_modes=None, nn_modes="all", exact=False, normalize_factors=False, verbose=False, return_errors=False, cvg_criterion="abs_rec_error", ): """ Non-negative CP decomposition via HALS Uses Hierarchical ALS (Alternating Least Squares) which updates each factor column-wise (one column at a time while keeping all other columns fixed), see [1]_ Parameters ---------- tensor : ndarray rank : int number of components n_iter_max : int maximum number of iteration init : {'svd', 'random'}, optional svd : str, default is 'truncated_svd' function to use to compute the SVD, acceptable values in tensorly.SVD_FUNS tol : float, optional tolerance: the algorithm stops when the variation in the reconstruction error is less than the tolerance Default: 1e-8 random_state : {None, int, np.random.RandomState} sparsity_coefficients: array of float (of length the number of modes) The sparsity coefficients on each factor. If set to None, the algorithm is computed without sparsity Default: None, fixed_modes: array of integers (between 0 and the number of modes) Has to be set not to update a factor, 0 and 1 for U and V respectively Default: None nn_modes: None, 'all' or array of integers (between 0 and the number of modes) Used to specify which modes to impose non-negativity constraints on. If 'all', then non-negativity is imposed on all modes. Default: 'all' exact: If it is True, the algorithm gives a results with high precision but it needs high computational cost. If it is False, the algorithm gives an approximate solution Default: False normalize_factors : if True, aggregate the weights of each factor in a 1D-tensor of shape (rank, ), which will contain the norms of the factors verbose: boolean Indicates whether the algorithm prints the successive reconstruction errors or not Default: False return_errors: boolean Indicates whether the algorithm should return all reconstruction errors and computation time of each iteration or not Default: False cvg_criterion : {'abs_rec_error', 'rec_error'}, optional Stopping criterion for ALS, works if `tol` is not None. If 'rec_error', ALS stops at current iteration if ``(previous rec_error - current rec_error) < tol``. If 'abs_rec_error', ALS terminates when `|previous rec_error - current rec_error| < tol`. sparsity : float or int random_state : {None, int, np.random.RandomState} Returns ------- factors : ndarray list list of positive factors of the CP decomposition element `i` is of shape ``(tensor.shape[i], rank)`` errors: list A list of reconstruction errors at each iteration of the algorithm. References ---------- .. [1] N. Gillis and F. Glineur, Accelerated Multiplicative Updates and Hierarchical ALS Algorithms for Nonnegative Matrix Factorization, Neural Computation 24 (4): 1085-1105, 2012. """ weights, factors = initialize_cp( tensor, rank, init=init, svd=svd, non_negative=True, random_state=random_state, normalize_factors=normalize_factors, ) norm_tensor = tl.norm(tensor, 2) n_modes = tl.ndim(tensor) if sparsity_coefficients is None or isinstance(sparsity_coefficients, float): sparsity_coefficients = [sparsity_coefficients] * n_modes if fixed_modes is None: fixed_modes = [] if nn_modes == "all": nn_modes = set(range(n_modes)) elif nn_modes is None: nn_modes = set() # Avoiding errors for fixed_value in fixed_modes: sparsity_coefficients[fixed_value] = None for mode in range(n_modes): if sparsity_coefficients[mode] is not None: warnings.warn("Sparsity coefficient is ignored in unconstrained modes.") # Generating the mode update sequence modes = [mode for mode in range(n_modes) if mode not in fixed_modes] # initialisation - declare local varaibles rec_errors = [] # Iteratation for iteration in range(n_iter_max): # One pass of least squares on each updated mode for mode in modes: # Computing Hadamard of cross-products pseudo_inverse = tl.ones((rank, rank), **tl.context(tensor)) for i, factor in enumerate(factors): if i != mode: pseudo_inverse = pseudo_inverse * tl.dot( tl.transpose(factor), factor ) pseudo_inverse = ( tl.reshape(weights, (-1, 1)) * pseudo_inverse * tl.reshape(weights, (1, -1)) ) mttkrp = unfolding_dot_khatri_rao(tensor, (weights, factors), mode) if mode in nn_modes: # Call the hals resolution with nnls, optimizing the current mode nn_factor, _, _, _ = hals_nnls( tl.transpose(mttkrp), pseudo_inverse, tl.transpose(factors[mode]), n_iter_max=100, sparsity_coefficient=sparsity_coefficients[mode], exact=exact, ) factors[mode] = tl.transpose(nn_factor) else: factor = tl.solve(tl.transpose(pseudo_inverse), tl.transpose(mttkrp)) factors[mode] = tl.transpose(factor) if normalize_factors and mode != modes[-1]: weights, factors = cp_normalize((weights, factors)) if tol: factors_norm = cp_norm((weights, factors)) iprod = tl.sum(tl.sum(mttkrp * factors[-1], axis=0)) rec_error = ( tl.sqrt(tl.abs(norm_tensor**2 + factors_norm**2 - 2 * iprod)) / norm_tensor ) rec_errors.append(rec_error) if iteration >= 1: rec_error_decrease = rec_errors[-2] - rec_errors[-1] if verbose: print( f"iteration {iteration}, reconstruction error: {rec_error}, decrease = {rec_error_decrease}" ) if cvg_criterion == "abs_rec_error": stop_flag = abs(rec_error_decrease) < tol elif cvg_criterion == "rec_error": stop_flag = rec_error_decrease < tol else: raise TypeError("Unknown convergence criterion") if stop_flag: if verbose: print(f"PARAFAC converged after {iteration} iterations") break else: if verbose: print(f"reconstruction error={rec_errors[-1]}") if normalize_factors: weights, factors = cp_normalize((weights, factors)) cp_tensor = CPTensor((weights, factors)) if return_errors: return cp_tensor, rec_errors else: return cp_tensor
class CP_NN(DecompositionMixin): """ Non-Negative Candecomp-Parafac decomposition via Alternating-Least Square Computes a rank-`rank` decomposition of `tensor` [1]_ such that, ``tensor = [|weights; factors[0], ..., factors[-1] |]``. Parameters ---------- tensor : ndarray rank : int Number of components. n_iter_max : int Maximum number of iteration init : {'svd', 'random'}, optional Type of factor matrix initialization. See `initialize_factors`. svd : str, default is 'truncated_svd' function to use to compute the SVD, acceptable values in tensorly.SVD_FUNS normalize_factors : if True, aggregate the weights of each factor in a 1D-tensor of shape (rank, ), which will contain the norms of the factors tol : float, optional (Default: 1e-6) Relative reconstruction error tolerance. The algorithm is considered to have found the global minimum when the reconstruction error is less than `tol`. random_state : {None, int, np.random.RandomState} verbose : int, optional Level of verbosity return_errors : bool, optional Activate return of iteration errors mask : ndarray array of booleans with the same shape as ``tensor`` should be 0 where the values are missing and 1 everywhere else. Note: if tensor is sparse, then mask should also be sparse with a fill value of 1 (or True). Allows for missing values [2]_ cvg_criterion : {'abs_rec_error', 'rec_error'}, optional Stopping criterion for ALS, works if `tol` is not None. If 'rec_error', ALS stops at current iteration if (previous rec_error - current rec_error) < tol. If 'abs_rec_error', ALS terminates when |previous rec_error - current rec_error| < tol. sparsity : float or int If `sparsity` is not None, we approximate tensor as a sum of low_rank_component and sparse_component, where low_rank_component = cp_to_tensor((weights, factors)). `sparsity` denotes desired fraction or number of non-zero elements in the sparse_component of the `tensor`. fixed_modes : list, default is None A list of modes for which the initial value is not modified. The last mode cannot be fixed due to error computation. svd_mask_repeats: int If using a tensor with masked values, this initializes using SVD multiple times to remove the effect of these missing values on the initialization. Returns ------- CPTensor : (weight, factors) * weights : 1D array of shape (rank, ) all ones if normalize_factors is False (default), weights of the (normalized) factors otherwise * factors : List of factors of the CP decomposition element `i` is of shape (tensor.shape[i], rank) * sparse_component : nD array of shape tensor.shape. Returns only if `sparsity` is not None. errors : list A list of reconstruction errors at each iteration of the algorithms. References ---------- .. [1] T.G.Kolda and B.W.Bader, "Tensor Decompositions and Applications", SIAM REVIEW, vol. 51, n. 3, pp. 455-500, 2009. .. [2] Tomasi, Giorgio, and Rasmus Bro. "PARAFAC and missing values." Chemometrics and Intelligent Laboratory Systems 75.2 (2005): 163-180. .. [3] R. Bro, "Multi-Way Analysis in the Food Industry: Models, Algorithms, and Applications", PhD., University of Amsterdam, 1998 """ def __init__( self, rank, n_iter_max=100, init="svd", svd="truncated_svd", tol=10e-7, random_state=None, verbose=0, normalize_factors=False, mask=None, cvg_criterion="abs_rec_error", fixed_modes=None, ): self.n_iter_max = n_iter_max self.init = init self.svd = svd self.tol = tol self.random_state = random_state self.verbose = verbose self.normalize_factors = normalize_factors self.mask = mask self.cvg_criterion = cvg_criterion self.fixed_modes = fixed_modes def fit_transform(self, tensor): """Decompose an input tensor Parameters ---------- tensor : tensorly tensor input tensor to decompose Returns ------- CPTensor decomposed tensor """ cp_tensor, errors = non_negative_parafac( tensor, n_iter_max=self.n_iter_max, init=self.init, svd=self.svd, tol=self.tol, random_state=self.random_state, verbose=self.verbose, normalize_factors=self.normalize_factors, mask=self.mask, cvg_criterion=self.cvg_criterion, fixed_modes=self.fixed_modes, return_errors=True, ) self.decomposition_ = cp_tensor self.errors_ = errors return self.decomposition_ def __repr__(self): return f"Rank-{self.rank} Non-Negative CP decomposition."
[docs] class CP_NN_HALS(DecompositionMixin): """ Non-Negative Candecomp-Parafac decomposition via Alternating-Least Square Computes a rank-`rank` decomposition of `tensor` [1]_ such that:: ``tensor = [|weights; factors[0], ..., factors[-1] |]``. Parameters ---------- tensor : ndarray rank : int Number of components. n_iter_max : int Maximum number of iteration init : {'svd', 'random'}, optional Type of factor matrix initialization. See `initialize_factors`. svd : str, default is 'truncated_svd' function to use to compute the SVD, acceptable values in tensorly.SVD_FUNS normalize_factors : if True, aggregate the weights of each factor in a 1D-tensor of shape (rank, ), which will contain the norms of the factors tol : float, optional (Default: 1e-6) Relative reconstruction error tolerance. The algorithm is considered to have found the global minimum when the reconstruction error is less than `tol`. random_state : {None, int, np.random.RandomState} verbose : int, optional Level of verbosity return_errors : bool, optional Activate return of iteration errors mask : ndarray array of booleans with the same shape as ``tensor`` should be 0 where the values are missing and 1 everywhere else. Note: if tensor is sparse, then mask should also be sparse with a fill value of 1 (or True). Allows for missing values [2]_ cvg_criterion : {'abs_rec_error', 'rec_error'}, optional Stopping criterion for ALS, works if `tol` is not None. If 'rec_error', ALS stops at current iteration if (previous rec_error - current rec_error) < tol. If 'abs_rec_error', ALS terminates when ``|previous rec_error - current rec_error| < tol``. sparsity : float or int If `sparsity` is not None, we approximate tensor as a sum of low_rank_component and sparse_component, where low_rank_component = cp_to_tensor((weights, factors)). `sparsity` denotes desired fraction or number of non-zero elements in the sparse_component of the `tensor`. fixed_modes : list, default is None A list of modes for which the initial value is not modified. The last mode cannot be fixed due to error computation. svd_mask_repeats: int If using a tensor with masked values, this initializes using SVD multiple times to remove the effect of these missing values on the initialization. Returns ------- CPTensor : (weight, factors) * weights : 1D array of shape (rank, ) all ones if normalize_factors is False (default), weights of the (normalized) factors otherwise * factors : List of factors of the CP decomposition element `i` is of shape (tensor.shape[i], rank) * sparse_component : nD array of shape tensor.shape. Returns only if `sparsity` is not None. errors : list A list of reconstruction errors at each iteration of the algorithms. References ---------- .. [1] T.G.Kolda and B.W.Bader, "Tensor Decompositions and Applications", SIAM REVIEW, vol. 51, n. 3, pp. 455-500, 2009. .. [2] Tomasi, Giorgio, and Rasmus Bro. "PARAFAC and missing values." Chemometrics and Intelligent Laboratory Systems 75.2 (2005): 163-180. .. [3] R. Bro, "Multi-Way Analysis in the Food Industry: Models, Algorithms, and Applications", PhD., University of Amsterdam, 1998 """ def __init__( self, rank, n_iter_max=100, init="svd", svd="truncated_svd", tol=10e-8, sparsity_coefficients=None, fixed_modes=None, nn_modes="all", exact=False, verbose=False, normalize_factors=False, cvg_criterion="abs_rec_error", random_state=None, ): self.rank = rank self.n_iter_max = n_iter_max self.init = init self.svd = svd self.tol = tol self.sparsity_coefficients = sparsity_coefficients self.random_state = random_state self.fixed_modes = fixed_modes self.nn_modes = nn_modes self.exact = exact self.verbose = verbose self.normalize_factors = normalize_factors self.cvg_criterion = cvg_criterion self.random_state = random_state
[docs] def fit_transform(self, tensor): """Decompose an input tensor Parameters ---------- tensor : tensorly tensor input tensor to decompose Returns ------- CPTensor decomposed tensor """ cp_tensor, errors = non_negative_parafac_hals( tensor, rank=self.rank, n_iter_max=self.n_iter_max, init=self.init, svd=self.svd, tol=self.tol, random_state=self.random_state, sparsity_coefficients=self.sparsity_coefficients, fixed_modes=self.fixed_modes, nn_modes=self.nn_modes, exact=self.exact, verbose=self.verbose, normalize_factors=self.normalize_factors, return_errors=True, cvg_criterion=self.cvg_criterion, ) self.decomposition_ = cp_tensor self.errors_ = errors return self.decomposition_
def __repr__(self): return f"Rank-{self.rank} Non-Negative CP decomposition."