PaperlayerView exact revision on arXiv

Source-linked AI summary

Tensor Decomposition for Signal Processing and Machine Learning

Nicholas D. Sidiropoulos, Lieven De Lathauwer, Xiao Fu, Kejun Huang, Evangelos E. Papalexakis, Christos Faloutsos

arXiv:1607.01668v2stat.MLcs.LGmath.NA

TL;DR

Tensor research lacks a single sufficiently comprehensive and deep entry point spanning its mathematics, algorithms, and applications. This overview supplies that entry point, covering tensor models, uniqueness, algorithms, statistical analysis, and signal-processing and machine-learning applications. It also highlights practical boundaries, including difficult rank determination, non-guaranteed ALS convergence, and NP-hard low-rank approximation.

  • Problem

    Researchers lacked a single overview deep enough to connect tensor mathematics, algorithms, and signal-processing and machine-learning applications.

  • Method

    The article synthesizes tensor rank, decomposition models, uniqueness, approximation, algorithms, statistical performance analysis, and applications across signal processing and machine learning.

  • Results

    The overview provides a technically deep starting point for readers with graduate-level matrix algebra and probability to begin tensor research or develop related algorithms and software.

  • Takeaways & Limitations

    Tensor decompositions offer useful uniqueness properties and broad applications, but their practical use requires attention to rank ambiguity, approximation difficulty, and algorithmic convergence.

Abstract

from arXiv · show

Tensors or {\em multi-way arrays} are functions of three or more indices $(i,j,k,\cdots)$ -- similar to matrices (two-way arrays), which are functions of two indices $(r,c)$ for (row,column). Tensors have a rich history, stretching over almost a century, and touching upon numerous disciplines; but they have only recently become ubiquitous in signal and data analytics at the confluence of signal processing, statistics, data mining and machine learning. This overview article aims to provide a good starting point for researchers and practitioners interested in learning about and working with tensors. As such, it focuses on fundamentals and motivation (using various application examples), aiming to strike an appropriate balance of breadth {\em and depth} that will enable someone having taken first graduate courses in matrix algebra and probability to get started doing research and/or developing tensor algorithms and software. Some background in applied optimization is useful but not strictly required. The material covered includes tensor rank and rank decomposition; basic tensor factorization models and their relationships and properties (including fairly good coverage of identifiability); broad coverage of algorithms ranging from alternating optimization to stochastic gradient; statistical performance analysis; and applications ranging from source separation to collaborative filtering, mixture and topic modeling, classification, and multilinear subspace learning.

I. INTRODUCTION

This article provides a comprehensive, technically deep entry point to tensor decompositions, connecting signal-processing and machine-learning perspectives. It develops foundational concepts, products, models, algorithms, statistical analysis, and applications.

  • Foundations: Tensors are multi-way arrays indexed by three or more indices, generalizing matrices beyond two-way structure.The article focuses mainly on third-order tensors while noting that the definitions generalize naturally to higher orders.
  • Challenges: Tensor methods combine advantages and difficulties: low-rank factorization can be essentially unique under mild conditions, while rank determination and approximation can be hard or ill-posed.Observed tensors commonly contain noise, so applications often seek a useful signal rank rather than the actual tensor rank.
  • Motivation: The article addresses the lack of a single sufficiently comprehensive and deep starting point for tensor-decomposition research.It targets readers with graduate-level matrix algebra and probability, while applied optimization background is useful but not required.
  • Motivation: It bridges signal-processing and machine-learning viewpoints, which emphasize factor columns and rows, respectively, in CPD models.For ratings tensors, machine-learning researchers interpret rows as parsimonious latent-space representations, whereas signal-processing researchers focus on rank-1 factors for separation.
  • Scope: The roadmap covers matrix preliminaries, tensor rank and decompositions, tensor products, factorization models, uniqueness, approximation, and algorithmic methods.Products include Kronecker and Khatri–Rao products; algorithms include alternating optimization, gradient descent, Gauss–Newton, ADMM, and stochastic gradient.
  • Scope: The article also treats statistical performance and applications in communications, speech separation, harmonic retrieval, collaborative filtering, mixture and topic modeling, classification, and subspace learning.Its statistical analysis focuses on the Cramér–Rao bound and efficient numerical computation.

A. Low-rank tensor approximation

Low-rank tensor models are estimated by minimizing a nonconvex least-squares objective, commonly with Alternating Least Squares (ALS). Tensor rank and approximation behave differently from matrices: rank can be difficult to determine, and best low-rank approximations may fail to exist.

  • Model fitting: The CP model estimates factor matrices A, B, and C from noisy tensor data using a least-squares criterion.The objective uses the squared Frobenius norm of the residual.
  • Model fitting: Although the model is trilinear and nonconvex jointly, fixing two factor matrices makes the problem linear in the remaining factor.This conditional linearity permits alternating updates across the factor matrices.
  • Optimization: Alternating Least Squares repeatedly updates the factor matrices through these conditional subproblems until convergence.ALS is widely used because it is easy to program for approximate low-rank tensor models.
  • Rank bounds: A third-order tensor of size I × J × K can have rank as high as min(IJ, JK, IK), and this upper bound is attainable.The construction decomposes each frontal slab into at most min(I, J) rank-1 factors and combines the factors across slabs.
  • Rank behavior: Tensor rank depends on the field: a random 2 × 2 × 2 tensor has rank 2 almost surely over C, while complex-conjugate eigenvalues can imply real rank 3.Thus, the same tensor may have different ranks over the real and complex fields.
  • Approximation pathology: A rank-3 tensor can have border rank 2, allowing arbitrarily accurate rank-two approximations whose components diverge and approximately cancel.This phenomenon makes low-rank tensor approximation potentially ill-posed and can occur in practical tensor fitting.

IV. UNIQUENESS, DEMYSTIFIED

Tensor CPD uniqueness is explained through eigenvalue-based constructions and Kruskal-rank conditions, with uniqueness understood up to permutation and scaling. The section also identifies rank-estimation boundaries and increasingly relaxed generic conditions.

  • Essential uniqueness: CPD essential uniqueness means the rank-1 terms are determined up to permutation, with factor matrices additionally ambiguous under reciprocal diagonal scalings.The decomposition equivalence is represented by a common permutation and diagonal scaling matrices whose product is the identity.
  • Rank specification: Overestimating rank can destroy uniqueness unless additional conditions hold, including full column rank of the pairwise Khatri–Rao products; underestimating rank prevents full decomposition.An overestimated decomposition may admit an alternative representation with fewer rank-1 terms.
  • Eigenvalue-based uniqueness: Distinct generalized eigenvalues enable recovery of factor information from a matrix pencil, linking constructive CPD uniqueness proofs to EVD and ESPRIT.Under the stated assumptions, generalized eigenvectors recover corresponding columns of the factor matrices; repeated eigenvalues destroy this uniqueness mechanism.
  • Kruskal conditions: Kruskal rank kC ≥2 is necessary for uniqueness in the two-slice setting, and full-column-rank A and B then ensure rank(X)=F and essential uniqueness.The same condition extends to K-slice tensors, where kC ≥2 and rA=rB=F provide the stated sufficient result.
  • Generic conditions: Generic uniqueness conditions can be substantially relaxed by retaining structural information, while ignoring that structure can lose about a factor of 2 in the bound.The section reports numerical verification of the necessary-and-sufficient generic bound up to F = 24 and notes that Kruskal’s bound is sharp but not exhaustive.

A. Tucker and CPD

The Tucker model generalizes matrix SVD using a core tensor and factor matrices, allowing mixed interactions among factor columns. Its orthogonal compact form provides structured approximation and compression, but tensor-specific approximation and optimization remain challenging.

  • MLSVD and approximation: Core diagonality is generally impossible without restricting tensor rank, whereas all-orthogonality is available without loss of generality.This preserves a useful SVD analogy, although simply truncating core slabs does not generally produce the best low multilinear-rank approximation.
  • Tucker model: Tucker represents a tensor through factor matrices and a core tensor whose entries encode interactions among all factor-column combinations.Unlike CPD, Tucker permits mixed products involving noncorresponding columns.
  • Tucker model: The Tucker model is not unique because linear transformations can be transferred between factor matrices and the core.Choosing orthonormal bases for the mode spans yields a reduced core of dimensions r1 × r2 × r3.
  • MLSVD and approximation: MLSVD uses right singular vectors of the tensor unfoldings as orthonormal factor matrices.The resulting core has orthogonal one-lower-order slabs whose Frobenius norms correspond to singular values of tensor matrix views.
  • MLSVD and approximation: Discarding small-magnitude core elements provides an error bound without reconstructing the approximated tensor.For mode ranks (r′1, r′2, r′3), the approximation error can be at most 3 times the minimal error.
  • Optimization: Tucker-ALS alternates factor updates, with its cost sequence non-increasing and reward sequence non-decreasing, guaranteeing convergence of those sequences.The method is conceptually similar to ALS for CPD.

C. Compression as preprocessing

Orthogonal Tucker compression can reduce a tensor before CPD computation when exact-model assumptions hold. The compressed core preserves the CPD rank and can be expanded back to factors of the original tensor, but constraints may not survive compression.

  • Compression as preprocessing: The Tucker core admits a CPD with rank(G) ≤ rank(X), providing a smaller representation for subsequent decomposition.The core is obtained through orthogonal factor transformations.
  • Compression as preprocessing: If the CPD of the core is essentially unique, then rank(G) = rank(X), preserving the decomposition rank under compression.The original factors are recovered from the compressed factors through the Tucker matrices, up to permutation and scaling.
  • Compression as preprocessing: A practical strategy is to compress X, compute the CPD of the smaller core, and then expand the resulting factors.This procedure assumes exact CPD and orthogonal Tucker models when relating the compressed and original factors.
  • Limitations: Orthogonal compression does not preserve constrained CPD models such as nonnegative, monotonic, or sparse factors.Constrained ALS can instead update Tucker factors and the core through constrained or regularized least squares.

VI. OTHER DECOMPOSITIONS

Higher-order decompositions address scalability, structured components, and joint analysis beyond basic CPD and Tucker models. The section also reviews computational strategies and convergence issues for ALS and related fitting methods.

  • Other decompositions: For order N with equal mode rank r, the Tucker core has r^N entries, creating a curse of dimensionality as N grows.Tensor-train format avoids index cycles and permits matrix-like computations using established numerical linear algebra techniques.
  • Other decompositions: Block Term Decompositions represent tensors as sums of low-multilinear-rank terms that can be unique even when components are not rank-1.Expanding the blocks into repeated rank-1 vectors yields the PARALIND form.
  • Other decompositions: Coupled decompositions jointly analyze multiple matrices or tensors, and coupling can transfer uniqueness to an otherwise non-unique matrix decomposition.PARAFAC2 is an early coupled variant based on a shared covariance matrix.
  • Computation: Efficient CPD updates exploit Hadamard products and avoid explicitly storing the Khatri–Rao product.For sparse data, the resulting computation can use (2F +1)NNZ(X) operations without additional intermediate memory.
  • Computation: Tucker-ALS similarly uses column generation to avoid intermediate memory growth and exploit sparsity when computing unfolded tensor products.Oblique Tucker updates require pseudo-inverses of relatively small matrices, subject to rank conditions.
  • Optimization: ALS monotonically reduces the cost but is not guaranteed to converge to a stationary point.Maximum block improvement and proximal regularization offer convergence-oriented alternatives, with additional computational or modeling considerations.

C. Quasi-Newton and Nonlinear Least Squares

The paper surveys derivative-based and stochastic optimization methods for tensor models, including line search, missing-data handling, and constrained fitting. These methods trade computational structure, convergence behavior, and memory efficiency across settings.

  • Quasi-Newton and Nonlinear Least Squares: Newton, quasi-Newton, and nonlinear least-squares methods use gradient, Hessian, or Jacobian structure to update tensor-model parameters.Inexact nonlinear least squares can use a parallel ALS iteration as a preconditioner.
  • Step-size selection: Exact line search for CPD can require 11 loss evaluations because the step-size objective is a sixth-degree polynomial.In practice, the paper typically favors a small or good-enough step size because exact search costs more than gradient computation.
  • Missing data: Missing-data fitting can use derivative-based methods, stochastic gradients over observed entries, or expectation-maximization-style imputation.Imputation decreases the cost monotonically but is inefficient in memory for very large sparse data.
  • Stochastic optimization: For sparse CPD data, stochastic gradient updates affect only the factor rows associated with the sampled tensor entry.The approach can have limited parallelization at high order and sublinear convergence.
  • Constraints: Constraints can restore identifiability, improve estimation accuracy, ensure interpretability, and mitigate ill-posedness in challenging CPD problems.Examples include non-negativity, orthogonality, probability-simplex, and linear constraints.
  • Constraints: Direct ADMM requires sophisticated parameter tuning and is not guaranteed to converge, whereas an AO-ADMM hybrid was developed as a convergence-oriented alternative.Variable splitting for coupled constraints can also slow convergence by dividing variables into smaller blocks.

VIII. CRAM´ER-RAO BOUND

This section presents the Cramér-Rao bound for CPD models, derives the Fisher information matrix, and develops an efficient pseudo-inverse computation for its structured rank deficiency.

  • Performance benchmarking: The Cramér-Rao bound lower-bounds the variance of unbiased estimators using the pseudo-inverse of the Fisher information matrix.It remains a benchmark when empirical bias is small, even if unbiasedness is difficult to prove.
  • FIM derivation: The CPD Fisher information matrix is derived from the Jacobian of the nonlinear model and partitioned into blocks for the factor matrices.The derivation uses commutation-matrix properties to obtain diagonal and off-diagonal blocks.
  • Rank deficiency: The CPD FIM has rank at most (I + J + K)F − 2F, so it is necessarily rank deficient.The deficiency reflects the model structure and motivates using a pseudo-inverse for the CRB.
  • Efficient computation: When the deficiency equals 2F, the pseudo-inverse can be computed efficiently by exploiting the null space and structured block matrices.The approach uses a null-space construction and a matrix inversion lemma rather than direct inversion.
  • Efficient computation: O(F^6) flops are required for the most expensive step, versus O((I + J + K)^3F^3) flops for direct pseudo-inversion.The improvement applies when F ≪ min(I, J, K).

IX. APPLICATIONS

The applications section shows how tensor decompositions address blind signal separation and source-separation problems by exploiting multiway structure and diversity.

  • Blind Multiuser CDMA: In ideal single-user DS-CDMA, the received chip-rate matrix is rank one and its spreading code and bit sequence can be recovered up to sign ambiguity.With noise and imperfections, the principal component is extracted using SVD.
  • Blind Multiuser CDMA: Known linearly independent spreading codes permit zero-forcing recovery, but unknown codes make rank-two-or-higher matrix factorization non-unique in general.This non-uniqueness motivates using CPD for blind multiuser communications.
  • Blind Multiuser CDMA: CPD can recover spreading codes, information bits, and path losses up to permutation and scaling ambiguities, even when transmitters outnumber spreading-code length under an identifiability condition.The discussion notes practical concerns including multipath, noise, and imperfect conditions.
  • Source separation: Source separation can use time-varying powers or different power spectra to provide diversity for estimating the mixing matrix.These assumptions yield partially symmetric CPD models suited to speech separation and related signal-processing applications.
  • Source separation: Tensor factorization also supports source separation in the power-spectrum domain, including applications in radio astronomy and dynamic spectrum access.

C. Harmonics

This section connects tensor models to harmonic retrieval, recommender systems, mixture modeling, and multilinear classification and subspace learning.

  • Harmonics: Harmonic mixtures can be rewritten as two- or three-slab CPD models with Vandermonde structure in one or more modes.Vandermonde structure across all three modes leads to multidimensional harmonic retrieval.
  • Collaborative filtering: Collaborative filtering models sparse user–movie ratings with low-dimensional latent descriptions whose inner products predict missing entries.Only a very small percentage of ratings may be available in practice.
  • Collaborative filtering: Context such as timestamps can be represented as additional tensor modes, with smoothness constraints modeling slowly changing preferences.Coupled matrix-tensor decomposition can incorporate side information that is not naturally an extra tensor slice.
  • Classification: Multilinear classifiers use low-rank tensor weights to keep parameter counts low relative to training samples and support higher-order feature interactions.
  • Mixture modeling: Mixture-model parameters can be estimated from symmetric CPD structure in third-order statistics, but third-order moments alone retain a scaling ambiguity.Jointly fitting third- and second-order, or third- and first-order, statistics removes that ambiguity.

F. Topic modeling

The section develops tensor-based topic modeling and multilinear subspace learning, emphasizing moment tensors, identifiability, and class-discriminative reduced representations.

  • Topic modeling: For single-topic documents, word co-occurrence probabilities form a symmetric non-negative CPD whose factors encode topic-specific conditional word distributions.The factors and topic probabilities can be estimated from second- and third-order co-occurrence statistics.
  • Topic modeling: For documents mixing topics, the third-order word-statistics tensor admits a symmetric Tucker decomposition with topic-word factors and a core determined by topic proportions.The topic proportions are drawn from a distribution over the probability simplex.
  • Topic modeling: Tucker models are generally non-identifiable because of linear transformation freedom, which sparsity in the factors or core can alleviate.
  • Multilinear subspace learning: Discriminative subspace learning seeks a dimensionality-reducing transformation that places same-class points close and different-class points far apart.The solution uses the eigenvectors associated with the smallest eigenvalues of a constructed symmetric matrix.
  • Multilinear subspace learning: For vectorized tensors, common mode bases define a Kronecker-structured dimensionality reduction whose low-dimensional cores are optimized for class separation.

A. Detailed proof of uniqueness via eigendecomposition

The proof establishes essential uniqueness of a CP decomposition by reducing tensor-slab relations to eigendecomposition and then resolving the resulting permutation and scaling ambiguities under rank conditions.

  • Eigendecomposition setup: Slab transformations produce square nonsingular correlation equations, enabling recovery through an eigendecomposition rather than directly handling possibly tall factors.The transformed matrices R1 and R2 are F × F and nonsingular.
  • Ambiguities and recovery: The eigendecomposition recovers the latent factor matrix only up to a common column permutation and nonsingular diagonal scaling.The same ambiguity propagates to eA, while B and C can then be recovered from the original slab equations.
  • Eigendecomposition setup: Distinct eigenvalues are necessary for uniqueness because repeated eigenvalues permit arbitrary linear combinations of the corresponding eigenvectors.The eigendecomposition therefore identifies components only when the eigenvalues in D are distinct.
  • Uniqueness argument: Under the stated rank and k-rank conditions, any alternative factorization must map each rank-1 column of A ⊙B to a unique corresponding column of the alternative product.The proof uses slab combinations, rank inequalities, and the permutation lemma to constrain the alternative factors.
  • Uniqueness conclusion: The resulting relation between alternative and original factors is a common permutation and scaling, and the remaining truncated rows of a tall C are consequently unique.This completes the proof of essential uniqueness under the theorem’s conditions.

C. Rank and k-rank of the Khatri–Rao product

The section characterizes rank and k-rank of Khatri–Rao products, showing both restrictive failure conditions and generic full-rank behavior under continuous random factors.

  • Rank bounds: If the Khatri–Rao product has deficient k-rank, its deficiency is bounded by the component k-ranks through kB⊙A ≥ min(F, kA + kB − 1).The proof excludes dependent-column sets below the relevant k-rank thresholds using Sylvester’s inequality.
  • Rank bounds: Full rank of A and B does not necessarily guarantee full rank of A ⊙B, as a Vandermonde example yields rank 5 for a 6-column product.With I = 3 and J = 2 the corresponding 6 × 6 product is full rank, whereas swapping dimensions produces rank 5.
  • Generic rank: For continuously distributed random A and B, the Khatri–Rao product has rank and k-rank min(IJ, F) almost surely.The proof establishes this by showing suitable determinants are nonzero analytic functions.
  • Product structure: The order of Khatri–Rao multiplications changes only the row order, not the resulting rank or k-rank.This invariance applies to the final product under reordered multiplications.

F. CRB for matrix and CP tensor factorization

The article reviews the Cramér–Rao bound for matrix and CP factorizations, emphasizing constrained models, Fisher information computation, and efficient handling of large matrices.

  • Motivation: The Cramér–Rao bound is a central estimation benchmark in signal processing and is asymptotically achievable by maximum-likelihood estimators at high SNR in many cases.The section uses the CRB as the reference point for statistical performance analysis.
  • Scope and contribution: The matrix and CP factorization Fisher information matrices are unaffected by constraints on latent factors, although constraints such as non-negativity can be crucial for identifiability.The section also discusses efficient pseudo-inversion to reduce the computational and memory burden of FIM calculations.

1) The Cram´er-Rao Bound:

This section derives the Cramér–Rao bound through Fisher information, then treats singular information matrices, equality constraints, Gaussian noise, and selected non-Gaussian noise models.

  • General CRB: The Fisher information matrix is defined from the parameterized measurement density and supplies the covariance lower bound under the regularity conditions.The derivation relies on interchanging differentiation and integration and on unbiased estimation.
  • Singular FIM: When the Fisher information matrix is singular, a generalized Schur-complement expression using Φ† remains a valid, though looser and theoretically unattainable, bound.The stated conditions include Φ ⪰0 and (I −ΦΦ†)G = 0.
  • Identifiability connection: Fisher-information singularity is generally neither implied by nor does it imply non-identifiability, although rank deficiency often reflects trivial ambiguities.Permutation and scaling ambiguities should be fixed consistently when comparing MSE with the generalized CRB.
  • Constraints: Equality constraints modify the CRB through an orthonormal basis spanning the null space of the constraint Jacobian, whereas inequality constraints roughly do not affect it.The constrained bound is formed using the matrix Q associated with the equality-constraint tangent space.
  • Noise models: For additive i.i.d. Gaussian noise, the Fisher information matrix has a simple model-specific form, and related non-Gaussian CRBs can be scaled versions of Gaussian counterparts.The discussion includes Laplacian and Cauchy noise as examples of alternative models.

2) Cram´er-Rao Bound for Matrix Factorization Models:

The matrix-factorization CRB is derived under Gaussian noise, accounting for a rank-deficient Fisher information matrix and exploiting its null-space structure for efficient pseudo-inversion.

  • Under Gaussian noise, the matrix-factorization model parameterizes W and H, whose stacked unknown vector has size (m+n)k.
  • The null space is explicitly characterized through linearly independent vectors formed from perturbations of the factor matrices.
  • Brute-force SVD pseudo-inversion is approximately cubic in the FIM dimension, while direct Woodbury or blockwise inversion is obstructed by singular components and complicated structure.
  • The proposed computation completes the range using a null-space basis, yielding a block-diagonal matrix whose diagonal blocks can be inverted efficiently.
  • The resulting CRB quantities βW and βH bound the average squared estimation errors for unbiased estimators of W and H.

3) Cram´er-Rao Bound for CP Factorization Models:

The CP-factorization CRB is derived separately from the matrix case because its Fisher information matrix has different null-space structure. The derivation identifies that structure and uses it to compute the pseudo-inverse more efficiently.

  • The CP model under Gaussian noise parameterizes factor matrices H1,...,HN, whose stacked unknown vector has size (n1+...+nN)k.
  • The CP Fisher information matrix is expressed in block form and alternatively as a block-diagonal-plus-low-rank structure using Jacobian and Khatri–Rao-product terms.
  • The CP Fisher information matrix has rank at most (n1+...+nN)k−(N−1)k when all factor matrices have full column rank.
  • The (N−1)k-dimensional null space arises from scaling ambiguities across CP-factor columns, unlike the k^2-dimensional deficiency in matrix factorization.
  • A null-space basis is constructed and combined with an efficient inversion identity to obtain the pseudo-inverse needed for the CP Cramér–Rao bound.
  • Unlike the matrix case, the completed CP matrix lacks a convenient block-diagonal structure, although its diagonal-plus-low-rank form still supports efficient inversion when K+EET is invertible.
  • The computational improvement applies when Nk<n1+...+nN; otherwise, directly pseudo-inverting the FIM is not considered a bad option.
Loading 1607.01668v2…