Higher-order singular value decomposition

Short description: Tensor decomposition

In multilinear algebra, the higher-order singular value decomposition (HOSVD) is a misnomer. There does not exist a single tensor decomposition that retains all the defining properties of the matrix SVD. The matrix SVD simultaneously yields a

• rank-𝑅 decomposition and
• orthonormal subspaces for the row and column spaces.

These properties are not realized within a single algorithm for higher-order tensors, but are instead realized by two distinct algorithmic developments and represent two distinct research directions. Harshman, as well as, the team of Carol and Chang proposed Canonical polyadic decomposition (CPD), which is a variant of the tensor rank decomposition, in which a tensor is approximated as a sum of K rank-1 tensors for a user-specified K. L. R. Tucker proposed a strategy for computing orthonormal subspaces for third order tensors. Aspects of these algorithms can be traced as far back as F. L. Hitchcock in 1928.

De Lathauwer et al. introduced clarity to the Tucker concepts, while Vasilescu and Terzopoulos introduced algorithmic clarity. Vasilescu and Terzopoulos introduced the M-mode SVD, which is the classic algorithm that is currently referred in the literature as the Tucker or the HOSVD. The Tucker approach and De Lathauwer's implementation are both sequential and rely on iterative procedures such as gradient descent or the power method. By contrast, the M-mode SVD provides a closed-form solution that can be executed sequentially and is well-suited for parallel computation.

This misattribution has had lasting impact on the scholarly record, obscuring the original source of a widely adopted algorithm, and complicating efforts to trace its development, reproduce results, and recognizing the respective contributions of different research efforts.

The term M-mode SVD accurately reflects the algorithm employed. It captures the actual computation, a set of SVDs on mode-flattenings without making assumptions about the structure of the core tensor or implying a rank decomposition.

Robust and L1-norm-based variants of this decomposition framework have since been proposed.

Definition

For the purpose of this article, the abstract tensor ${𝒜}$ is assumed to be given in coordinates with respect to some basis as a M-way array, also denoted by ${𝒜}\in {\mathbb{ℂ}}^{I_1\times I_2\cdots \times \cdots I_m\cdots \times I_M}$, where M is the number of modes and the order of the tensor. $\mathbb{ℂ}$ is the complex numbers and it includes both the real numbers $\mathbb{ℝ}$ and the pure imaginary numbers.

Let ${{𝒜}}_{[m]}\in {\mathbb{ℂ}}^{I_m\times (I_1{I}_2\cdots I_{m-1}{I}_{m+1}\cdots I_M)}$ denote the mode-m flattening of ${𝒜}$, so that the left index of ${{𝒜}}_{[m]}$ corresponds to the $m$'th index ${𝒜}$ and the right index of ${{𝒜}}_{[m]}$ corresponds to all other indices of ${𝒜}$ combined. Let ${\mathbf{𝐔}}_m\in {\mathbb{ℂ}}^{I_m\times I_m}$be a unitary matrix containing a basis of the left singular vectors of the ${{𝒜}}_{[m]}$ such that the jth column ${\mathbf{𝐮}}_j$ of ${\mathbf{𝐔}}_m$ corresponds to the jth largest singular value of ${{𝒜}}_{[m]}$. Observe that the mode/factor matrix ${\mathbf{𝐔}}_m$ does not depend on the particular on the specific definition of the mode m flattening. By the properties of the multilinear multiplication, we have$$\begin{array}{rcl}{𝒜}& =& {𝒜}\times (\mathbf{𝐈},\mathbf{𝐈},\dots ,\mathbf{𝐈})\\ & =& {𝒜}\times ({\mathbf{𝐔}}_1{\mathbf{𝐔}}_1^H,{\mathbf{𝐔}}_2{\mathbf{𝐔}}_2^H,\dots ,{\mathbf{𝐔}}_M{\mathbf{𝐔}}_M^H)\\ & =& ({𝒜}\times ({\mathbf{𝐔}}_1^H,{\mathbf{𝐔}}_2^H,\dots ,{\mathbf{𝐔}}_M^H))\times ({\mathbf{𝐔}}_1,{\mathbf{𝐔}}_2,\dots ,{\mathbf{𝐔}}_M),\end{array}$$where ${\cdot }^H$ denotes the conjugate transpose. The second equality is because the ${\mathbf{𝐔}}_m$'s are unitary matrices. Define now the core tensor$${𝒮}:={𝒜}\times ({\mathbf{𝐔}}_1^H,{\mathbf{𝐔}}_2^H,\dots ,{\mathbf{𝐔}}_M^H).$$Then, the M-mode SVD(HOSVD) of ${𝒜}$ is the decomposition$${𝒜}={𝒮}\times ({\mathbf{𝐔}}_1,{\mathbf{𝐔}}_2,\dots ,{\mathbf{𝐔}}_M).$$ The above construction shows that every tensor has a M-mode SVD(HOSVD).

Compact M-mode SVD (HOSVD)

As in the case of the compact singular value decomposition of a matrix, where the rows and columns corresponding to vanishing singular values are dropped, it is also possible to consider a compact M-mode SVD(HOSVD), which is very useful in applications.

Assume that ${\mathbf{𝐔}}_m\in {\mathbb{ℂ}}^{I_m\times R_m}$ is a matrix with unitary columns containing a basis of the left singular vectors corresponding to the nonzero singular values of the standard factor-m flattening ${{𝒜}}_{[m]}$ of ${𝒜}$. Let the columns of ${\mathbf{𝐔}}_m$ be sorted such that the $r_m$ th column ${\mathbf{𝐮}}_{r_m}$ of ${\mathbf{𝐔}}_m$ corresponds to the $r_m$th largest nonzero singular value of ${{𝒜}}_{[m]}$. Since the columns of ${\mathbf{𝐔}}_m$ form a basis for the image of ${{𝒜}}_{[m]}$, we have$${{𝒜}}_{[m]}={\mathbf{𝐔}}_m{\mathbf{𝐔}}_m^H{{𝒜}}_{[m]}=({𝒜}{\times }_m({\mathbf{𝐔}}_m{\mathbf{𝐔}}_m^H){)}_{[m]},$$where the first equality is due to the properties of orthogonal projections (in the Hermitian inner product) and the last equality is due to the properties of multilinear multiplication. As flattenings are bijective maps and the above formula is valid for all $m=1,2,\dots ,m,\dots ,M$, we find as before that$$\begin{array}{rcl}{𝒜}& =& {𝒜}\times ({\mathbf{𝐔}}_1{\mathbf{𝐔}}_1^H,{\mathbf{𝐔}}_2{\mathbf{𝐔}}_2^H,\dots ,{\mathbf{𝐔}}_M{\mathbf{𝐔}}_M^H)\\ & =& ({𝒜}\times ({\mathbf{𝐔}}_1^H,{\mathbf{𝐔}}_2^H,\dots ,{\mathbf{𝐔}}_M^H))\times ({\mathbf{𝐔}}_1,{\mathbf{𝐔}}_2,\dots ,{\mathbf{𝐔}}_M)\\ & =& {𝒮}\times ({\mathbf{𝐔}}_1,{\mathbf{𝐔}}_2,\dots ,{\mathbf{𝐔}}_M),\end{array}$$where the core tensor ${𝒮}$ is now of size $R_1\times R_2\times \cdots \times R_M$.

Multilinear rank

The multilinear rank of ${𝒜}$ is denoted with rank-$(R_1,R_2,\dots ,R_M)$. The multilinear rank is a tuple in ${\mathbb{\mathbb{N}}}^M$ where $R_m:=\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}({{𝒜}}_{[m]})$. Not all tuples in ${\mathbb{\mathbb{N}}}^M$ are multilinear ranks. The multilinear ranks are bounded by $1\le R_m\le I_m$ and it satisfy the constraint $R_m\le \prod _{i\ne m}{R}_i$ must hold.

The compact M-mode SVD(HOSVD) is a rank-revealing decomposition in the sense that the dimensions of its core tensor correspond with the components of the multilinear rank of the tensor.

Interpretation

The following geometric interpretation is valid for both the full and compact M-mode SVD(HOSVD). Let $(R_1,R_2,\dots ,R_M)$ be the multilinear rank of the tensor ${𝒜}$. Since ${𝒮}\in {\mathbb{ℂ}}^{R_1\times R_2\times \cdots \times R_M}$ is a multidimensional array, we can expand it as follows$${𝒮}=\sum _{r_1=1}^{R_1}\sum _{r_2=1}^{R_2}\cdots \sum _{r_M=1}^{R_M}{s}_{r_1,r_2,\dots ,r_M}{\mathbf{𝐞}}_{r_1}\otimes {\mathbf{𝐞}}_{r_2}\otimes \cdots \otimes {\mathbf{𝐞}}_{r_M},$$where ${\mathbf{𝐞}}_{r_m}$ is the $r_m$th standard basis vector of ${\mathbb{ℂ}}^{I_m}$. By definition of the multilinear multiplication, it holds that$${𝒜}=\sum _{r_1=1}^{R_1}\sum _{r_2=1}^{R_2}\cdots \sum _{r_M=1}^{R_M}{s}_{r_1,r_2,\dots ,r_M}{\mathbf{𝐮}}_{r_1}\otimes {\mathbf{𝐮}}_{r_2}\otimes \cdots \otimes {\mathbf{𝐮}}_{r_M},$$where the ${\mathbf{𝐮}}_{r_m}$ are the columns of ${\mathbf{𝐔}}_m\in {\mathbb{ℂ}}^{I_m\times R_m}$. It is easy to verify that $B=\{{\mathbf{𝐮}}_{r_1}\otimes {\mathbf{𝐮}}_{r_2}\otimes \cdots \otimes {\mathbf{𝐮}}_{r_M}{\}}_{r_1,r_2,\dots ,r_M}$ is an orthonormal set of tensors. This means that the M-mode SVD(HOSVD) can be interpreted as a way to express the tensor ${𝒜}$ with respect to a specifically chosen orthonormal basis $B$ with the coefficients given as the multidimensional array ${𝒮}$.

Computation

Let ${𝒜}\in {\mathbb{ℂ}}^{I_1\times I_2\times \cdots \times I_M}$ be a tensor with a rank-$(R_1,R_2,\dots ,R_M)$, where $\mathbb{ℂ}$ contains the reals $\mathbb{ℝ}$ as a subset.

Classic computation

While De Lathauwer et al. clarified Tucker's concepts through two influential papers, Vasilescu and Terzopoulos provided algoritmic clarity. The Tucker algorithm and De Lathauwer et al. companion algorithm are sequential, relying on iterative methods such as gradient descent or the power method. In contrast, the M-mode SVD is computes the othonormal subspaces in closed-form that can be executed sequentially, but it is also well-suited for parallel computation.

M-mode SVD (also referred to as HOSVD or Tucker)

What is commonly referred to as the HOSVD or Tucker was developed by Vasilescu and Terzopoulos under the name M-mode SVD.

• For $m=1,\dots ,M$, do the following:

1. Construct the mode-m flattening ${{𝒜}}_{[m]}$;
2. Compute the (compact) singular value decomposition
   ${{𝒜}}_{[m]}={\mathbf{𝐔}}_m{\mathrm{\Sigma }}_m{\mathbf{𝐕}}_m^T$, and store the left singular vectors $\mathbf{𝐔}\in {\mathbb{ℂ}}^{I_m\times R_m}$;

• Compute the core tensor ${𝒮}$ via the mode-m product
  ${𝒮}={𝒜}{\times }_1{\mathbf{𝐔}}_1^H{\times }_2{\mathbf{𝐔}}_2^H\dots {\times }_m{\mathbf{𝐔}}_m^H\dots {\times }_M{\mathbf{𝐔}}_M^H$

Interlacing computation

A strategy that is significantly faster when some or all $R_m\ll I_m$ consists of interlacing the computation of the core tensor and the factor matrices, as follows:

• Set ${{𝒜}}^0={𝒜}$;
• For $m=1,2\dots ,M$ perform the following:
  1. Construct the standard mode-m flattening ${{𝒜}}_{[m]}^{m-1}$;
  2. Compute the (compact) singular value decomposition ${{𝒜}}_{[m]}^{m-1}=U_m{\mathrm{\Sigma }}_m{V}_m^T$, and store the left singular vectors $U_m\in F^{I_m\times R_m}$;
  3. Set ${{𝒜}}^m=U_m^H{\times }_m{{𝒜}}^{m-1}$, or, equivalently, ${{𝒜}}_{[m]}^m={\mathrm{\Sigma }}_m{V}_m^T$.

In-place computation

The M-mode SVD (HOSVD) can be computed in-place via the Fused In-place Sequentially Truncated Higher Order Singular Value Decomposition (FIST-HOSVD) algorithm by overwriting the original tensor by the M-mode SVD (HOSVD) core tensor, significantly reducing the memory consumption of computing HOSVD.

Approximation

In applications, such as those mentioned below, a common problem consists of approximating a given tensor ${𝒜}\in {\mathbb{ℂ}}^{I_1\times I_2\times \cdots \times I_m\cdots \times I_M}$ by one with a reduced multilinear rank. Formally, if the multilinear rank of ${𝒜}$ is denoted by $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\mathrm{-}(R_1,R_2,\dots ,R_m,\dots ,R_M)$, then computing the optimal $\overline{{A}}$ that approximates ${𝒜}$ for a given reduced $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\mathrm{-}({\overline{R}}_1,{\overline{R}}_2,\dots ,{\overline{R}}_m,\dots ,{\overline{R}}_M)$ is a nonlinear non-convex ${\ell }_2$-optimization problem $$\underset{\overline{{A}}\in {\mathbb{ℂ}}^{I_1\times I_2\times \cdots \times I_M}}{\min }\frac{1}{2}\Vert {𝒜}-\overline{{A}}{\Vert }_F^2\text{s.t.}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}\mathrm{-}({\overline{R}}_1,{\overline{R}}_2,\dots ,{\overline{R}}_M),$$where $({\overline{R}}_1,{\overline{R}}_2,\dots ,{\overline{R}}_M)\in {\mathbb{\mathbb{N}}}^M$ is the reduced multilinear rank with $1\le {\overline{R}}_m<R_m\le I_m$, and the norm $\Vert .{\Vert }_F$ is the Frobenius norm.

A simple idea for trying to solve this optimization problem is to truncate the (compact) SVD in step 2 of either the classic or the interlaced computation. A classically truncated M-mode SVD/HOSVD is obtained by replacing step 2 in the classic computation by

• Compute a rank-${\overline{R}}_m$ truncated SVD ${{𝒜}}_{[m]}\approx U_m{\mathrm{\Sigma }}_m{V}_m^T$, and store the top ${\overline{R}}_m$ left singular vectors $U_m\in F^{I_m\times {\overline{R}}_m}$;

while a sequentially truncated M-mode SVD (HOSVD) (or successively truncated M-mode SVD(HOSVD)) is obtained by replacing step 2 in the interlaced computation by

• Compute a rank-${\overline{R}}_m$ truncated SVD ${{𝒜}}_{[m]}^{m-1}\approx U_m{\mathrm{\Sigma }}_m{V}_m^T$, and store the top ${\overline{R}}_m$ left singular vectors $U_m\in F^{I_m\times {\overline{R}}_m}$. Unfortunately, truncation does not result in an optimal solution for the best low multilinear rank optimization problem,. However, both the classically and interleaved truncated M-mode SVD/HOSVD result in a quasi-optimal solution: if ${\overline{{A}}}_t$ denotes the classically or sequentially truncated M-mode SVD(HOSVD) and ${\overline{{A}}}^{*}$ denotes the optimal solution to the best low multilinear rank approximation problem, then$$\Vert {𝒜}-{\overline{{A}}}_t{\Vert }_F\le \sqrt{M}\Vert {𝒜}-{\overline{{A}}}^{*}{\Vert }_F;$$in practice this means that if there exists an optimal solution with a small error, then a truncated M-mode SVD/HOSVD will for many intended purposes also yield a sufficiently good solution.

Applications

The M-mode SVD (HOSVD/Tucker) is most commonly applied to the extraction of relevant information from multi-way arrays.

Starting in the early 2000s, Vasilescu addressed causal questions by reframing the data analysis, recognition and synthesis problems as multilinear tensor problems. The power of the tensor framework was showcased by decomposing and representing an image in terms of its causal factors of data formation, in the context of Human Motion Signatures for gait recognition, face recognition—TensorFaces and computer graphics—TensorTextures.

The M-mode SVD (HOSVD) has been successfully applied to signal processing and big data, e.g., in genomic signal processing. These applications also inspired a higher-order GSVD (HO GSVD) and a tensor GSVD.

A combination of M-mode SVD (HOSVD) and SVD also has been applied for real-time event detection from complex data streams (multivariate data with space and time dimensions) in disease surveillance.

It is also used in tensor product model transformation-based controller design.

The concept of M-mode SVD (HOSVD) was carried over to functions by Baranyi and Yam via the TP model transformation. This extension led to the definition of the M-mode SVD/HOSVD canonical form of tensor product functions and Linear Parameter Varying system models and to convex hull manipulation based control optimization theory, see TP model transformation in control theories.

M-mode SVD (HOSVD) was proposed to be applied to multi-way data analysis in an unsupervised manner and was successfully applied to in silico drug discovery from gene expression.

Robust L1-norm variant

L1-Tucker is the L1-norm-based, robust variant of Tucker decomposition. L1-HOSVD is the analogous of M-mode SVD(HOSVD) for the solution to L1-Tucker.