import numpy as np import warnings import tensorly as tl from ..random import check_random_state from ..base import unfold from ..kruskal_tensor import kruskal_to_tensor from ..tenalg import khatri_rao # Author: Jean Kossaifi # Author: Chris Swierczewski # Author: Sam Schneider # License: BSD 3 clause def normalize_factors(factors): """Normalizes factors to unit length and returns factor magnitudes 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 ---------- factors : ndarray list 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 ------- normalized_factors : list of ndarrays list of matrices with the same shape as `factors` weights : ndarray vector of length `R` holding normalizing constants """ # allocate variables for weights, and normalized factors rank = factors[0].shape[1] weights = tl.ones(rank, **tl.context(factors[0])) normalized_factors = [] # normalize columns of factor matrices for factor in factors: scales = tl.norm(factor, axis=0) weights *= scales scales_non_zero = tl.where(scales==0, tl.ones(tl.shape(scales), **tl.context(factors[0])), scales) normalized_factors.append(factor/scales_non_zero) return normalized_factors, weights def initialize_factors(tensor, rank, init='svd', svd='numpy_svd', random_state=None, non_negative=False): r"""Initialize factors used in `parafac`. The type of initialization is set using `init`. If `init == 'random'` then initialize factor matrices using `random_state`. If `init == 'svd'` then initialize the `m`th factor matrix using the `rank` left singular vectors of the `m`th unfolding of the input tensor. Parameters ---------- tensor : ndarray rank : int init : {'svd', 'random'}, optional svd : str, default is 'numpy_svd' function to use to compute the SVD, acceptable values in tensorly.SVD_FUNS non_negative : bool, default is False if True, non-negative factors are returned Returns ------- factors : ndarray list List of initialized factors of the CP decomposition where element `i` is of shape (tensor.shape[i], rank) """ rng = check_random_state(random_state) if init == 'random': factors = [tl.tensor(rng.random_sample((tensor.shape[i], rank)), **tl.context(tensor)) for i in range(tl.ndim(tensor))] if non_negative: return [tl.abs(f) for f in factors] else: return factors elif init == 'svd': try: svd_fun = tl.SVD_FUNS[svd] except KeyError: message = 'Got svd={}. However, for the current backend ({}), the possible choices are {}'.format( svd, tl.get_backend(), tl.SVD_FUNS) raise ValueError(message) factors = [] for mode in range(tl.ndim(tensor)): U, _, _ = svd_fun(unfold(tensor, mode), n_eigenvecs=rank) if tensor.shape[mode] < rank: # TODO: this is a hack but it seems to do the job for now # factor = tl.tensor(np.zeros((U.shape[0], rank)), **tl.context(tensor)) # factor[:, tensor.shape[mode]:] = tl.tensor(rng.random_sample((U.shape[0], rank - tl.shape(tensor)[mode])), **tl.context(tensor)) # factor[:, :tensor.shape[mode]] = U random_part = tl.tensor(rng.random_sample((U.shape[0], rank - tl.shape(tensor)[mode])), **tl.context(tensor)) U = tl.concatenate([U, random_part], axis=1) if non_negative: factors.append(tl.abs(U[:, :rank])) else: factors.append(U[:, :rank]) return factors raise ValueError('Initialization method "{}" not recognized'.format(init)) def parafac(tensor, rank, n_iter_max=100, init='svd', svd='numpy_svd', tol=1e-8, orthogonalise=False, random_state=None, verbose=False, return_errors=False): """CANDECOMP/PARAFAC decomposition via alternating least squares (ALS) Computes a rank-`rank` decomposition of `tensor` [1]_ such that, ``tensor = [| 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 'numpy_svd' function to use to compute the SVD, acceptable values in tensorly.SVD_FUNS 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 Returns ------- factors : ndarray list List of 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 algorithms. References ---------- .. [1] tl.G.Kolda and B.W.Bader, "Tensor Decompositions and Applications", SIAM REVIEW, vol. 51, n. 3, pp. 455-500, 2009. """ if orthogonalise and not isinstance(orthogonalise, int): orthogonalise = n_iter_max factors = initialize_factors(tensor, rank, init=init, svd=svd, random_state=random_state) rec_errors = [] norm_tensor = tl.norm(tensor, 2) for iteration in range(n_iter_max): if orthogonalise and iteration <= orthogonalise: factor = [tl.qr(factor)[0] for factor in factors] for mode in range(tl.ndim(tensor)): pseudo_inverse = tl.tensor(np.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) factor = tl.dot(unfold(tensor, mode), khatri_rao(factors, skip_matrix=mode)) factor = tl.transpose(tl.solve(tl.transpose(pseudo_inverse), tl.transpose(factor))) factors[mode] = factor if tol: rec_error = tl.norm(tensor - kruskal_to_tensor(factors), 2) / norm_tensor rec_errors.append(rec_error) if iteration > 1: if verbose: print('reconstruction error={}, variation={}.'.format( rec_errors[-1], rec_errors[-2] - rec_errors[-1])) if tol and abs(rec_errors[-2] - rec_errors[-1]) < tol: if verbose: print('converged in {} iterations.'.format(iteration)) break if return_errors: return factors, rec_errors else: return factors def non_negative_parafac(tensor, rank, n_iter_max=100, init='svd', svd='numpy_svd', tol=10e-7, random_state=None, verbose=0): """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 'numpy_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 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 = 10e-12 nn_factors = initialize_factors(tensor, rank, init=init, svd=svd, random_state=random_state, non_negative=True) n_factors = len(nn_factors) norm_tensor = tl.norm(tensor, 2) rec_errors = [] for iteration in range(n_iter_max): for mode in range(tl.ndim(tensor)): # khatri_rao(factors).tl.dot(khatri_rao(factors)) # simplifies to multiplications sub_indices = [i for i in range(n_factors) if i != mode] for i, e in enumerate(sub_indices): if i: accum = accum*tl.dot(tl.transpose(nn_factors[e]), nn_factors[e]) else: accum = tl.dot(tl.transpose(nn_factors[e]), nn_factors[e]) numerator = tl.dot(unfold(tensor, mode), khatri_rao(nn_factors, skip_matrix=mode)) numerator = tl.clip(numerator, a_min=epsilon, a_max=None) denominator = tl.dot(nn_factors[mode], accum) denominator = tl.clip(denominator, a_min=epsilon, a_max=None) nn_factors[mode] = nn_factors[mode]* numerator / denominator rec_error = tl.norm(tensor - kruskal_to_tensor(nn_factors), 2) / norm_tensor rec_errors.append(rec_error) if iteration > 1 and verbose: print('reconstruction error={}, variation={}.'.format( rec_errors[-1], rec_errors[-2] - rec_errors[-1])) if iteration > 1 and abs(rec_errors[-2] - rec_errors[-1]) < tol: if verbose: print('converged in {} iterations.'.format(iteration)) break return nn_factors def sample_khatri_rao(matrices, n_samples, skip_matrix=None, return_sampled_rows=False, random_state=None): """Random subsample of the Khatri-Rao product of the given list of matrices If one matrix only is given, that matrix is directly returned. Parameters ---------- matrices : ndarray list list of matrices with the same number of columns, i.e.:: for i in len(matrices): matrices[i].shape = (n_i, m) n_samples : int number of samples to be taken from the Khatri-Rao product skip_matrix : None or int, optional, default is None if not None, index of a matrix to skip random_state : None, int or numpy.random.RandomState if int, used to set the seed of the random number generator if numpy.random.RandomState, used to generate random_samples returned_sampled_rows : bool, default is False if True, also returns a list of the rows sampled from the full khatri-rao product Returns ------- sampled_Khatri_Rao : ndarray The sampled matricised tensor Khatri-Rao with `n_samples` rows indices : tuple list a list of indices sampled for each mode indices_kr : int list list of length `n_samples` containing the sampled row indices """ if random_state is None or not isinstance(random_state, np.random.RandomState): rng = check_random_state(random_state) warnings.warn('You are creating a new random number generator at each call.\n' 'If you are calling sample_khatri_rao inside a loop this will be slow:' ' best to create a rng outside and pass it as argument (random_state=rng).') else: rng = random_state if skip_matrix is not None: matrices = [matrices[i] for i in range(len(matrices)) if i != skip_matrix] n_factors = len(matrices) rank = tl.shape(matrices[0])[1] sizes = [tl.shape(m)[0] for m in matrices] # For each matrix, randomly choose n_samples indices for which to compute the khatri-rao product indices_list = [rng.randint(0, tl.shape(m)[0], size=n_samples, dtype=int) for m in matrices] if return_sampled_rows: # Compute corresponding rows of the full khatri-rao product indices_kr = np.zeros((n_samples), dtype=int) for size, indices in zip(sizes, indices_list): indices_kr = indices_kr*size + indices # Compute the Khatri-Rao product for the chosen indices sampled_kr = tl.ones((n_samples, rank), **tl.context(matrices[0])) for indices, matrix in zip(indices_list, matrices): sampled_kr = sampled_kr*matrix[indices, :] if return_sampled_rows: return sampled_kr, indices_list, indices_kr else: return sampled_kr, indices_list def randomised_parafac(tensor, rank, n_samples, n_iter_max=100, init='random', svd='numpy_svd', tol=10e-9, max_stagnation=20, random_state=None, verbose=1): """Randomised CP decomposition via sampled ALS Parameters ---------- tensor : ndarray rank : int number of components n_samples : int number of samples per ALS step n_iter_max : int maximum number of iteration init : {'svd', 'random'}, optional svd : str, default is 'numpy_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 max_stagnation: int, optional, default is 0 if not zero, the maximum allowed number of iterations with no decrease in fit random_state : {None, int, np.random.RandomState}, default is None verbose : int, optional level of verbosity Returns ------- factors : ndarray list list of positive factors of the CP decomposition element `i` is of shape ``(tensor.shape[i], rank)`` References ---------- .. [3] Casey Battaglino, Grey Ballard and Tamara G. Kolda, "A Practical Randomized CP Tensor Decomposition", """ rng = check_random_state(random_state) factors = initialize_factors(tensor, rank, init=init, svd=svd, random_state=random_state) rec_errors = [] n_dims = tl.ndim(tensor) norm_tensor = tl.norm(tensor, 2) min_error = 0 for iteration in range(n_iter_max): for mode in range(n_dims): kr_prod, indices_list = sample_khatri_rao(factors, n_samples, skip_matrix=mode, random_state=rng) indices_list = [i.tolist() for i in indices_list] # Keep all the elements of the currently considered mode indices_list.insert(mode, slice(None, None, None)) # MXNet will not be happy if this is a list insteaf of a tuple indices_list = tuple(indices_list) if mode: sampled_unfolding = tensor[indices_list] else: sampled_unfolding = tl.transpose(tensor[indices_list]) pseudo_inverse = tl.dot(tl.transpose(kr_prod), kr_prod) factor = tl.dot(tl.transpose(kr_prod), sampled_unfolding) factor = tl.transpose(tl.solve(pseudo_inverse, factor)) factors[mode] = factor if max_stagnation or tol: rec_error = tl.norm(tensor - kruskal_to_tensor(factors), 2) / norm_tensor if not min_error or rec_error < min_error: min_error = rec_error stagnation = -1 stagnation += 1 rec_errors.append(rec_error) if iteration > 1: if verbose: print('reconstruction error={}, variation={}.'.format( rec_errors[-1], rec_errors[-2] - rec_errors[-1])) if (tol and abs(rec_errors[-2] - rec_errors[-1]) < tol) or \ (stagnation and (stagnation > max_stagnation)): if verbose: print('converged in {} iterations.'.format(iteration)) break return factors