SIAM Journal on Matrix Analysis and Applications

This issue of SIAM Journal on Matrix Analysis and Applications was motivated by the Workshop on Tensor Decompositions and Applications, held in Luminy, France from August 29 to September 2, 2005. The issue was announced through both the SIMAX and the workshop web sites.

Though decompositions of higher-order tensors have been around for many years, the door is now opening on greater mathematical understanding and new applications. Tensor decompositions have been applied by researchers in psychometrics and chemometrics since the seventies. More recently, tensors have found their way to signal processing via the use of high-order statistics, joint matrix techniques, and applications in telecommunications. Other applications include complexity theory, (blind) system identification, biomedical engineering, numerical analysis, and data mining, among others.

This volume contains 15 papers that together form a nice cross-section of current research on tensor decompositions. The papers present new algorithms, basic algebraic results including the introduction of new decompositions, and applications in signal processing, scientific computing, and large-scale problems.

Many thanks are due to all the authors for their valuable contributions. We would also like to thank Henk van der Vorst, Mitch Chernoff, and the SIAM staff for their support.

Since the initial submission of the papers, two of the authors have passed away. Richard Harshman and Gene Golub were both authorities in their field.

We consider Tucker-like approximations with an $r \times r \times r$ core tensor for three-dimensional $n \times n \times n$ arrays in the case of $r \ll n$ and possibly very large n (up to $10^4$–$10^6$). As the approximation contains only $\mathcal{O}(rn + r^3)$ parameters, it is natural to ask if it can be computed using only a small amount of entries of the given array. A similar question for matrices (two-dimensional tensors) was asked and positively answered in [S. A. Goreinov, E. E. Tyrtyshnikov, and N. L. Zamarashkin, A theory of pseudo-skeleton approximations, Linear Algebra Appl., 261 (1997), pp. 1–21]. In the present paper we extend the positive answer to the case of three-dimensional tensors. More specifically, it is shown that if the tensor admits a good Tucker approximation for some (small) rank r, then this approximation can be computed using only $\mathcal{O}(nr)$ entries with $\mathcal{O}(nr^{3})$ complexity.

Motivated by numerous applications in which the data may be modeled by a variable subscripted by three or more indices, we develop a tensor-based extension of the matrix CUR decomposition. The tensor-CUR decomposition is most relevant as a data analysis tool when the data consist of one mode that is qualitatively different from the others. In this case, the tensor-CUR decomposition approximately expresses the original data tensor in terms of a basis consisting of underlying subtensors that are actual data elements and thus that have a natural interpretation in terms of the processes generating the data. Assume the data may be modeled as a $(2+1)$-tensor, i.e., an $m \times n \times p$ tensor $\mathcal{A}$ in which the first two modes are similar and the third is qualitatively different. We refer to each of the p different $m \times n$ matrices as “slabs” and each of the $mn$ different p-vectors as “fibers.” In this case, the tensor-CUR algorithm computes an approximation to the data tensor $\mathcal{A}$ that is of the form $\mathcal{CUR}$, where $\mathcal{C}$ is an $m \times n \times c$ tensor consisting of a small number c of the slabs, $\mathcal{R}$ is an $r \times p$ matrix consisting of a small number r of the fibers, and $\mathcal{U}$ is an appropriately defined and easily computed $c \times r$ encoding matrix. Both $\mathcal{C}$ and $\mathcal{R}$ may be chosen by randomly sampling either slabs or fibers according to a judiciously chosen and data-dependent probability distribution, and both c and r depend on a rank parameter k, an error parameter $\epsilon$, and a failure probability $\delta$. Under appropriate assumptions, provable bounds on the Frobenius norm of the error tensor $\mathcal{A} - \mathcal{CUR}$ are obtained. In order to demonstrate the general applicability of this tensor decomposition, we apply it to problems in two diverse domains of data analysis: hyperspectral medical image analysis and consumer recommendation system analysis. In the hyperspectral data application, the tensor-CUR decomposition is used to compress the data, and we show that classification quality is not substantially reduced even after substantial data compression. In the recommendation system application, the tensor-CUR decomposition is used to reconstruct missing entries in a user-product-product preference tensor, and we show that high quality recommendations can be made on the basis of a small number of basis users and a small number of product-product comparisons from a new user.

We consider the low-rank approximation over the real field of generic $p \times q \times 2$ arrays. For all possible combinations of p, q, and R, we present conjectures on the existence of a best rank-R approximation. Our conjectures are motivated by a detailed analysis of the boundary of the set of arrays with at most rank R. We link these results to the Candecomp/Parafac (CP) model for three-way component analysis. Essentially, CP tries to find a best rank-R approximation to a given three-way array. In the case of $p \times q \times 2$ arrays, we show (under some regularity condition) that if a best rank-R approximation does not exist, then any sequence of CP updates will exhibit diverging CP components, which implies that several components are highly correlated in all three modes and their component weights become arbitrarily large. This extends Stegeman [Psychometrika, 71 (2006), pp. 483–501], who considers $p \times p \times 2$ arrays of rank $p+1$ or higher. We illustrate our results by means of simulations.

The higher-order singular value decomposition (HOSVD) is a generalization of the singular value decomposition (SVD) to higher-order tensors (i.e., arrays with more than two indices) and plays an important role in various domains. Unfortunately, this decomposition is computationally demanding. Indeed, the HOSVD of a third-order tensor involves the computation of the SVD of three matrices, which are referred to as “modes” or “matrix unfoldings.” In this paper, we present fast algorithms for computing the full and the rank-truncated HOSVD of third-order structured (symmetric, Toeplitz, and Hankel) tensors. These algorithms are derived by considering two specific ways to unfold a structured tensor, leading to structured matrix unfoldings whose SVD can be efficiently computed.

In this paper we study a generalization of Kruskal's permutation lemma to partitioned matrices. We define the k'-rank of partitioned matrices as a generalization of the k-rank of matrices. We derive a lower-bound on the k'-rank of Khatri–Rao products of partitioned matrices. We prove that Khatri–Rao products of partitioned matrices are generically full column rank.

In this paper we introduce a new class of tensor decompositions. Intuitively, we decompose a given tensor block into blocks of smaller size, where the size is characterized by a set of mode-n ranks. We study different types of such decompositions. For each type we derive conditions under which essential uniqueness is guaranteed. The parallel factor decomposition and Tucker's decomposition can be considered as special cases in the new framework. The paper sheds new light on fundamental aspects of tensor algebra.

In this paper we derive alternating least squares algorithms for the computation of the block term decompositions introduced in Part II. We show that degeneracy can also occur for block term decompositions.

There has been continued interest in seeking a theorem describing optimal low-rank approximations to tensors of order 3 or higher that parallels the Eckart–Young theorem for matrices. In this paper, we argue that the naive approach to this problem is doomed to failure because, unlike matrices, tensors of order 3 or higher can fail to have best rank-r approximations. The phenomenon is much more widespread than one might suspect: examples of this failure can be constructed over a wide range of dimensions, orders, and ranks, regardless of the choice of norm (or even Brègman divergence). Moreover, we show that in many instances these counterexamples have positive volume: they cannot be regarded as isolated phenomena.  In one extreme case, we exhibit a tensor space in which no rank-3 tensor has an optimal rank-2 approximation. The notable exceptions to this misbehavior are rank-1 tensors and order-2 tensors (i.e., matrices). In a more positive spirit, we propose a natural way of overcoming the ill-posedness of the low-rank approximation problem, by using weak solutions when true solutions do not exist. For this to work, it is necessary to characterize the set of weak solutions, and we do this  in the case of rank 2, order 3 (in arbitrary dimensions). In our work we emphasize the importance of closely studying concrete low-dimensional examples as a first step toward more general results. To this end, we present a detailed analysis of equivalence classes of $2 \times 2 \times 2$ tensors, and we develop methods for extending results upward to higher orders and dimensions. Finally, we link our work to existing studies of tensors from an algebraic geometric point of view. The rank of a tensor can in theory be given a semialgebraic description; in other words, it can be determined by a system of polynomial inequalities. We study some of these polynomials in cases of interest to us; in particular, we make extensive use of the hyperdeterminant $\Delta$ on $\mathbb{R}^{2\times 2 \times 2}$.

Several modifications have been proposed to speed up the alternating least squares (ALS) method of fitting the PARAFAC model. The most widely used is line search, which extrapolates from linear trends in the parameter changes over prior iterations to estimate the parameter values that would be obtained after many additional ALS iterations. We propose some extensions of this approach that incorporate a more sophisticated extrapolation, using information on nonlinear trends in the parameters and changing all the parameter sets simultaneously. The new method, called “enhanced line search (ELS),” can be implemented at different levels of complexity, depending on how many different extrapolation parameters (for different modes) are jointly optimized during each iteration. We report some tests of the simplest parameter version, using simulated data. The performance of this lowest-level of ELS depends on the nature of the convergence difficulty. It significantly outperforms standard LS when there is a “convergence bottleneck,” a situation where some modes have almost collinear factors but others do not, but is somewhat less effective in classic “swamp” situations where factors are highly collinear in all modes. This is illustrated by examples. To demonstrate how ELS can be adapted to different N-way decompositions, we also apply it to a four-way array to perform a blind identification of an under-determined mixture (UDM). Since analysis of this dataset happens to involve a serious convergence “bottleneck” (collinear factors in two of the four modes), it provides another example of a situation in which ELS dramatically outperforms standard line search.

We investigate the sensitivity of the problem of nonorthogonal (matrix) joint diagonalization (NOJD). First, we consider the uniqueness conditions for the problem of exact joint diagonalization (EJD), which is closely related to the issue of uniqueness in tensor decompositions. As a byproduct, we derive the well-known identifiability conditions for independent component analysis (ICA) based on an EJD formulation of ICA. We next introduce some known cost functions for NOJD and derive flows based on these cost functions for NOJD. Then we define and investigate the noise sensitivity of the stationary points of these flows. We show that the condition number of the joint diagonalizer and uniqueness of the joint diagonalizer as measured by modulus of uniqueness (as defined in this paper) affect the sensitivity. We also investigate the effect of the number of matrices on the sensitivity. Our numerical experiments confirm the theoretical results.While this paper was under review a few of its results were presented in ICASSP07 conference in Honolulu, HI .

This paper presents some recent filtering methods based on the lower-rank tensor approximation approach for denoising tensor signals. In this approach, multicomponent data are represented by tensors, that is, multiway arrays, and the presented tensor filtering methods rely on multilinear algebra. First, the classical channel-by-channel SVD-based filtering method is overviewed. Then, an extension of the classical matrix filtering method is presented. It is based on the lower rank-$(K_1,\dots,K_N)$ truncation of the higher order SVD which performs a multimode principal component analysis (PCA) and is implicitly developed for an additive white Gaussian noise. Two tensor filtering methods recently developed by the authors are also overviewed. The first method consists of an improvement of the multimode PCA-based tensor filtering in the case of an additive correlated Gaussian noise. This improvement is specially done thanks to the fourth order cumulant slice matrix. The second method consists of an extension of Wiener filtering for data tensors. The performances and comparative results between all these tensor filtering methods are presented for the cases of noise reduction in color images, multispectral images, and multicomponent seismic data.

In this paper we present a Newton method to jointly approximately diagonalize a set of positive definite Hermitian matrices. To this end, we derive the local gradient and Hessian of the underlying cost function in closed form. The algorithm is derived for the complex case and can also update a nonsquare diagonalization matrix. We analyze the cost function at the critical points and show its relation to a different cost function that is commonly studied.

Suppose $\mathcal{A}=(a_{ijk})\in\mathbb{R}^{n\times n\times n}$ is a three-way array or third-order tensor. Many of the powerful tools of linear algebra such as the singular value decomposition (SVD) do not, unfortunately, extend in a straightforward way to tensors of order three or higher. In the two-dimensional case, the SVD is particularly illuminating, since it reduces a matrix to diagonal form. Although it is not possible in general to diagonalize a tensor (i.e., $a_{ijk}=0$ unless $i=j=k$), our goal is to “condense” a tensor in fewer nonzero entries using orthogonal transformations. We propose an algorithm for tensors of the form $\mathcal{A}\in\mathbb{R}^{n\times n\times n}$ that is an extension of the Jacobi SVD algorithm for matrices. The resulting tensor decomposition reduces $\mathcal{A}$ to a form such that the quantity $\sum_{i=1}^n a_{iii}^2$ or $\sum_{i=1}^n a_{iii}$ is maximized.

The Kronecker tensor-product approximation combined with the $\mathcal{H}$-matrix techniques provides an efficient tool to represent integral operators as well as a discrete elliptic operator inverse $A^{-1}\in\mathbb{R}^{N\times N}$ in $\mathbb{R}^{d}$ (the discrete Green's function) with a high spatial dimension d. In the present paper we give a survey on modern methods of the structured tensor-product approximation to multidimensional integral operators and Green's functions and present some new results on the existence of low tensor-rank decompositions to a class of function-related operators. The memory space of the considered data-sparse representations is estimated by $\mathcal{O}(dn\log^{q}n)$ with q independent of d, retaining the approximation accuracy of order $O(n^{-\delta})$, where $n=N^{1/d}$ is the dimension of the discrete problem in one space direction. In particular, we apply the results to the Newton, Yukawa, and Helmholtz kernels $\frac{1}{|x-y|}$, $\frac{e^{-\lambda|x-y|}}{|x-y|}$, and $\frac{\cos(\lambda|x-y|)}{|x-y|}$, respectively, with $x,y \in \mathbb{R}^{d}$.

A symmetric tensor is a higher order generalization of a symmetric matrix. In this paper, we study various properties of symmetric tensors in relation to a decomposition into a symmetric sum of outer product of vectors. A rank-1 order-k tensor is the outer product of k nonzero vectors. Any symmetric tensor can be decomposed into a linear combination of rank-1 tensors, each of which is symmetric or not. The rank of a symmetric tensor is the minimal number of rank-1 tensors that is necessary to reconstruct it. The symmetric rank is obtained when the constituting rank-1 tensors are imposed to be themselves symmetric. It is shown that rank and symmetric rank are equal in a number of cases and that they always exist in an algebraically closed field. We will discuss the notion of the generic symmetric rank, which, due to the work of Alexander and Hirschowitz [J. Algebraic Geom., 4 (1995), pp. 201–222], is now known for any values of dimension and order. We will also show that the set of symmetric tensors of symmetric rank at most r is not closed unless $r=1$.