Recent Advances in Tensor Decomposition

A tensor is a multi-way array of numbers. An order-1 tensor is just a vector $u \in \mathbb{R}^n$. An order-2 tensor is a matrix $M \in \mathbb{R}^{n_1 \times n_2}$. An order-3 tensor is a 3-way array $T \in \mathbb{R}^{n_1 \times n_2 \times n_3}$, and so on.

Analogous to the rank-1 matrix $uv^\top$ for (column) vectors $u$ and $v$, a rank-1 tensor of shape $n_1 \times n_2 \times n_3$ is written $u \otimes v \otimes w$ for vectors $u \in \mathbb{R}^{n_1}$, $v \in \mathbb{R}^{n_2}$, $w \in \mathbb{R}^{n_3}$, and this notation means the tensor has entries $T_{ijk} = u_i v_j w_k$.

Tensors have a bunch of applications in statistics and data science [see 6, 8, 10, 11]. For instance, certain datasets might naturally be represented as a 3-way array encoding 3-way interactions between 3 different variables. Another key example in statistics is the method of moments: given many samples of an $n$-dimensional random vector, it may be useful to compute the moments. The first moment is the mean (expected value), which is an $n$-dimensional vector; the second moment is the $n \times n$ covariance matrix; the third moment is an $n \times n \times n$ tensor; and so on. We are used to performing various primitive computations on matrices: eigenvalues and eigenvectors, singular value decomposition (SVD), low-rank approximation, and so on. For tensors of order 3 and above, the analogous operations tend to be much more difficult to compute, or even ill-defined.

One of the most basic (and useful) computational tasks involving tensors is to decompose a given $n_1 \times n_2 \times n_3$ tensor $T$ as a sum of rank-1 terms:

\begin{equation}\label{eq:decomp}
T = \sum_{i=1}^r u^{(i)} \otimes v^{(i)} \otimes w^{(i)} \qquad \text{where } u^{(i)} \in \mathbb{R}^{n_1}, \, v^{(i)} \in \mathbb{R}^{n_2}, \, w^{(i)} \in \mathbb{R}^{n_3}.
\end{equation}

Here, addition of tensors is defined entry-wise as for matrices: $(T+S)_{ijk} = T_{ijk} + S_{ijk}$. The smallest $r$ for which $T$ admits a decomposition into $r$ rank-1 terms is called the rank of $T$, analogous to the usual notion of matrix rank. For a rank-$r$ matrix, a decomposition into $r$ rank-1 terms can be found efficiently on a computer using the SVD, but the solution is not unique due to the "rotation problem" [See 8]. In contrast, rank-$r$ tensors of order 3 and above tend to admit a unique decomposition into $r$ rank-1 terms, as long as $r$ is small enough. (Of course, one can re-order the $r$ terms and represent a rank-1 tensor in multiple ways: $u \otimes v \otimes w = (2u) \otimes (v/2) \otimes w$. So the meaning of "unique" here is that the collection of rank-1 summands is unique.) This uniqueness is key to why tensor decomposition is a useful subroutine in many applications!

For such a basic problem, the algorithmic question of decomposing order-3 tensors (with the minimal number of terms) is still not so well understood. In fact, the problem is known to be NP-hard, meaning we do not expect there to exist a fast (polynomial-time) algorithm that solves the problem for all input tensors [3, 4]. However, there are algorithms that succeed for "almost all" input tensors of low enough rank. Let's assume for simplicity that the input $T$ is $n \times n \times n$, and assume further that it is generated as the sum of $r$ rank-1 tensors as in [1], where the components $u^{(i)},v^{(i)},w^{(i)}$ are generically chosen. Here, "generic" has a formal meaning based in algebraic geometry, but essentially this means there can be a measure-zero set of 
"bad" components for which the algorithm fails, but otherwise it should succeed at finding the rank-$r$ decomposition.

What is known about the generic $n \times n \times n$ decomposition task defined above? As long as $r \le cn^2$ for a constant $c > 0$, the decomposition that generated $T$ will be the unique rank-$r$ decomposition, and the rank of $T$ will be exactly $r$ [1]. This means the solution is identifiable in principle, but does not guarantee a fast algorithm. (A related fact: no $n \times n \times n$ tensor can have rank larger than $n^2$, since each of the $n$ "slices" is an $n \times n$ matrix which has rank at most $n$.) As far as algorithms are concerned, a classical method (circa 1970) called simultaneous diagonalization or Jennrich's algorithm succeeds when $r \le n$ [see 8]. This barrier remained unbroken for many years, and the first progress happened here at the UC Davis Department of Mathematics: Ph.D. student Haolin Chen and his advisor Luis Rademacher gave an algorithm that succeeds when $r = n+k$ for any constant $k$ [2]. More recently, Pascal Koiran (ENS Lyon) improved this condition to $r \le 4n/3$ [5]. In forthcoming joint work with Pravesh Kothari (Princeton) and Ankur Moitra (MIT), we improve this further to $r \le (2-\epsilon)n$ for any constant $\epsilon > 0$. All of these algorithms involve only basic linear algebraic operations such as Gram-Schmidt, matrix diagonalization, and the like, but these are combined in somewhat subtle ways. The runtime is polynomial in $n$, but exponential in the "constants" $k$ or $\epsilon$.

To explain the key ideas in our new algorithm, it is instructive to think about a related but easier task known as rank detection, which is a necessary prerequisite for decomposition. We are again given a tensor $T$ that is constructed as the sum of $r$ rank-1 terms with generic components, but now the goal is simply to determine $r$. Building on the Ph.D. thesis of Madalina Persu (MIT) [9], we approach this by constructing a "flattening," that is, a map that takes a tensor and outputs a matrix, hopefully with the property that the matrix rank reflects the tensor rank in a predictable way. One basic example of a flattening is the so-called "trivial" flattening where the $n$ "slices" of the tensor are stacked to form an $n^2 \times n$ matrix as shown in the figure.

One can check that the trivial flattening transforms a rank-1 tensor into a rank-1 matrix. The flattening is a linear map, so a rank-$r$ tensor flattens to the sum of $r$ rank-1 matrices, which we might hope has rank exactly $r$. This indeed pans out if the tensor components are generic and $r \le n$, allowing rank detection to be solved for $r \le n$. When $r > n$, this approach breaks down because the $n^2 \times n$ matrix cannot have rank exceeding its smaller dimension $n$. The core component of our new algorithm is a certain non-trivial flattening, based on the so-called Koszul-Young flattening [7], which can solve rank detection up to rank $r \le (2-\epsilon)n$.

So now we have new algorithms for rank detection and decomposition that tolerate higher rank tensors, by roughly a factor of 2, than the classical methods. The big question remaining is whether this can be improved further, since, after all, this is still quite far from the identifiability threshold $r \approx n^2$. I suspect that rank $2n$ might actually be a fundamental breaking point for all fast algorithms, and we can prove that a certain class of flattening-based approaches cannot surpass this. On the other hand, maybe there is a completely different strategy waiting to be discovered, that could break this barrier!