Randomized CP Tensor Decomposition

TL;DR

This paper introduces a randomized algorithm for efficient approximate CP tensor decomposition, enabling scalable analysis of large tensors with controlled error, supported by theoretical and empirical validation.

Contribution

It presents a novel randomized method that significantly reduces computation time for CP tensor decomposition on large-scale data.

Findings

Algorithm achieves faster decomposition on massive tensors.

Approximation error is controllable through oversampling and power iterations.

Empirical results validate the effectiveness of the proposed method.

Abstract

The CANDECOMP/PARAFAC (CP) tensor decomposition is a popular dimensionality-reduction method for multiway data. Dimensionality reduction is often sought after since many high-dimensional tensors have low intrinsic rank relative to the dimension of the ambient measurement space. However, the emergence of `big data' poses significant computational challenges for computing this fundamental tensor decomposition. By leveraging modern randomized algorithms, we demonstrate that coherent structures can be learned from a smaller representation of the tensor in a fraction of the time. Thus, this simple but powerful algorithm enables one to compute the approximate CP decomposition even for massive tensors. The approximation error can thereby be controlled via oversampling and the computation of power iterations. In addition to theoretical results, several empirical results demonstrate the performance of the proposed algorithm.

Figures

Tables

Table 1. Table 1 : Summary of the computational results for the noisy toy video.

	Parameters	Time (s)	Speedup	Iterations	Error
CP BCD	$k=4$	2.31	-	9	0.0171
rCP BCD	$k=4$, $p=10$, $q=0$	0.13	17	9	0.494
	$k=4$, $p=10$, $q=1$	0.14	16	10	0.0191
	$k=4$, $p=10$, $q=2$	0.15	15	10	0.0164
SVD	$k=4$	0.52	-	-	0.137

Table 2. Table 2 : Summary of the computational results for the noise-free cylinder flow.

	Parameters	Time (s)	Speedup	Iterations	Error
CP BCD	$k=30$	115.55	-	458	0.117
rCP BCD	$k=30$, $p=10$, $q=0$	1.27	91	533	0.122
	$k=30$, $p=10$, $q=1$	1.41	82	517	0.121
	$k=30$, $p=10$, $q=2$	1.56	74	437	0.118
SVD	$k=30$	0.57	-	-	4.25E-05

Table 3. Table 3 : Summary of the computational results for the noise-corrupted cylinder flow.

	Parameters	Time (s)	Speedup	Iterations	Error
CP BCD	$k=30$	64.01	-	239	0.191
rCP BCD	$k=30$, $p=10$, $q=0$	0.99	64	332	0.522
	$k=30$, $p=10$, $q=1$	1.23	52	414	0.189
	$k=30$, $p=10$, $q=2$	1.13	56	370	0.153
SVD	$k=30$	0.58	-	-	0.655
	$k=6$	0.58	-	-	0.311

Equations

$$\mathcal{X}\approx \sum _{r=1}^R{\mathbf{a}}_r\circ {\mathbf{b}}_r\circ {\mathbf{c}}_r,$$

$$\mathcal{X}\approx \sum _{r=1}^R{\lambda }_r\cdot {\mathbf{a}}_r\circ {\mathbf{b}}_r\circ {\mathbf{c}}_r.$$

$$\mathcal{X}\approx [[\mathbf{\lambda };\mathbf{A},\mathbf{B},\mathbf{C}]].$$

$$\begin{array}[]{cccc}\mathbf{A}&\approx&\mathbf{Q}&\mathbf{B}\\[-3.61371pt] m\times n&&m\times k&k\times n\end{array}$$

$$\mathbf{Y}=\mathbf{A\Omega },$$

$$\mathbf{A}\approx \mathbf{Q}{\mathbf{Q}}^{\top }\mathbf{A}$$

$$\mathbf{B}={\mathbf{Q}}^{\top }\mathbf{A},$$

$$\underset{\hat{\mathcal{X}}}{\mathrm{minimize}}$$

$$\mathcal{X}\approx \mathcal{X}{\times }_1{\mathbf{Q}}_1{\mathbf{Q}}_1^{\top }{\times }_2\cdots {\times }_N{\mathbf{Q}}_N{\mathbf{Q}}_N^{\top }.$$

$$\mathcal{M}=\mathcal{X}{\times }_n{\mathbf{Q}}_n{\mathbf{Q}}_n^{\top }$$

$${\mathbf{Y}}_n={\mathcal{X}}_{(n)}{\mathbf{\Omega }}_n,$$

$${\mathcal{B}}_n=\mathcal{X}{\times }_n{\mathbf{Q}}_n^{\top }$$

$${\mathcal{X}}_{(n)}^q:=({\mathcal{X}}_{(n)}{\mathcal{X}}_{(n)}^{\top }{)}^q{\mathcal{X}}_{(n)},$$

$${\mathbf{Y}}_n=({\mathcal{X}}_{(n)}{\mathcal{X}}_{(n)}^{\top }{)}^q{\mathcal{X}}_{(n)}{\mathbf{\Omega }}_n.$$

$$E\Vert \mathcal{E}{\Vert }_F=\Vert \mathcal{X}-\hat{\mathcal{X}}{\Vert }_F,$$

$$E\Vert \mathcal{X}-\mathcal{X}{\times }_1{\mathbf{Q}}_1{\mathbf{Q}}_1^{\top }{\times }_2\cdots {\times }_N{\mathbf{Q}}_N{\mathbf{Q}}_N^{\top }{\Vert }_F$$

$$\underset{\hat{\mathcal{B}}}{\mathrm{minimize}}$$

$$\mathbf{A}\approx {\mathbf{Q}}_1\tilde{\mathbf{A}},\mathbf{B}\approx {\mathbf{Q}}_2\tilde{\mathbf{B}},\mathbf{C}\approx {\mathbf{Q}}_3\tilde{\mathbf{C}},$$

$$\underset{\tilde{\mathbf{A}},\tilde{\mathbf{B}},\tilde{\mathbf{C}}}{\mathrm{minimize}}$$

$${\mathcal{B}}_{(1)}\approx \tilde{\mathbf{A}}(\tilde{\mathbf{C}}\odot \tilde{\mathbf{B}}{)}^{\top },{\mathcal{B}}_{(2)}\approx \tilde{\mathbf{B}}(\tilde{\mathbf{C}}\odot \tilde{\mathbf{A}}{)}^{\top },{\mathcal{B}}_{(3)}\approx \tilde{\mathbf{C}}(\tilde{\mathbf{B}}\odot \tilde{\mathbf{A}}{)}^{\top },$$

$${\tilde{\mathbf{A}}}^{j+1}$$

$${\tilde{\mathbf{B}}}^{j+1}$$

$${\tilde{\mathbf{C}}}^{j+1}$$

$$(\mathbf{A}\odot \mathbf{B}{)}^{†}=({\mathbf{A}}^{\top }\mathbf{A}\ast {\mathbf{B}}^{\top }\mathbf{B}{)}^{†}(\mathbf{A}\odot \mathbf{B}{)}^{\top },$$

$${\mathcal{Y}}_{\text{res}}=\mathcal{B}-\sum _{r=1}^{\tilde{R}}{\tilde{\mathbf{a}}}_r\circ {\tilde{\mathbf{b}}}_r\circ {\tilde{\mathbf{c}}}_r.$$

$${\tilde{\mathbf{a}}}_r^{j+1}$$

$${\tilde{\mathbf{b}}}_r^{j+1}$$

$${\tilde{\mathbf{c}}}_r^{j+1}$$

$$c_{SVD}$$

$$c_{CP}$$

$$\Vert \mathcal{E}{\Vert }_F=\Vert \mathcal{X}-\hat{\mathcal{X}}\Vert =\Vert \mathcal{X}-\mathcal{X}{\times }_1{\mathbf{Q}}_1{\mathbf{Q}}_1^{\top }{\times }_2\cdots {\times }_N{\mathbf{Q}}_N{\mathbf{Q}}_N^{\top }{\Vert }_F.$$

$$\Vert \mathcal{E}{\Vert }_F=\Vert \mathcal{X}-\mathcal{X}{\times }_1{\mathbf{P}}_1{\times }_2\cdots {\times }_N{\mathbf{P}}_N{\Vert }_F.$$

$$\Vert \mathcal{E}{\Vert }_F=\sum _{n=1}^N\Vert \mathcal{X}-\mathcal{X}{\times }_n{\mathbf{P}}_n{\Vert }_F=\sum _{n=1}^N\Vert \mathcal{X}{\times }_n(\mathbf{I}-{\mathbf{P}}_n){\Vert }_F,$$

Full text

Keywords: randomized algorithms; randomized least squares; dimension reduction; multilinear algebra; CP decomposition; canonical polyadic tensor decomposition.

1 Introduction

Advances in data acquisition and storage technology have enabled the acquisition of massive amounts of data in a wide range of emerging applications. In particular, numerous applications across the physical, biological, social and engineering sciences generate large multidimensional, multi-relational and/or multi-modal data. Efficient analysis of this data requires dimensionality reduction techniques. However, traditionally employed matrix decompositions techniques such as the singular value decomposition (SVD) and principal component analysis (PCA) can become inadequate when dealing with multidimensional data. This is because reshaping multi-modal data into matrices, or data flattening, can fail to reveal important structures in the data.

Tensor decompositions overcome the information loss from flattening. The canonical CP tensor decomposition expresses an $N$-way tensor as a sum of rank-one tensors to extract multi-modal structure. It is particularly suitable for data-driven discovery, as shown by (Hong et al., 2020) for various learning tasks on real world data. However, tensor decompositions of massive multidimensional data pose a tremendous computational challenge. Hence, innovations that reduce the computational demands have become increasingly relevant in this field. The idea of tensor compression, for instance, eases computational bottlenecks by constructing smaller (compressed) tensors, which are then used as a proxy to efficiently compute CP decompositions. Compressed tensors may be obtained, for example, using the Tucker decomposition of a tensor into one small core tensor and $N$ unitary matrices However, this approach requires the expensive computation of the left singular vectors for each mode.

Related Work.

This computational challenge can be tackled using modern randomized techniques developed to compute the SVD. Tsourakakis (2010) presents a randomized Tucker decomposition algorithm based on the random sparsification idea of Achlioptas and McSherry (2007) for computing the SVD. Zhou et al. (2014) proposed an accelerated randomized CP algorithm using the ideas of Halko et al. (2011) without the use of power iterations. An alternative randomized tensor algorithm proposed by Drineas and Mahoney (2007) uses random column selection. Related work by Vervliet et al. (2014) proposed a sparsity-promoting algorithm for incomplete tensors using compressed sensing.

Another efficient approach to building large-scale tensor decompositions is the subdivision of a tensor into a set of blocks. These smaller blocks can then be used to approximate the CP decomposition of the full tensor in a parallelized or distributed fashion (Phan and Cichocki, 2011). Sidiropoulos et al. (2014) fused random projection and blocking into a highly computationally efficient algorithm. More recently, Vervliet and Lathauwer (2016) proposed a block sampling CP decomposition method for the analysis of large-scale tensors using randomization, demonstrating significant computational savings while attaining near optimal accuracy. These block based algorithms are particularly relevant if the tensor is too large to fit into fast memory.

Contributions.

We present the randomized CP algorithm, which is closely related to the randomized methods of Halko et al. (2011). Our method proceeds in two stages. The first stage uses random projections with power iterations to obtain a compressed tensor (Algorithm 1). The second stage applies either alternating least squares (ALS) or block coordinate descent (BCD) to the compressed tensor (Algorithms 2 and 3), at significantly reduced cost. Finally, the CP decomposition of the original tensor is obtained by projecting the compressed factor matrices back using the basis derived in Algorithm 1.

Randomized algorithms for CP decomposition have been proposed (Zhou et al., 2014; Battaglino et al., 2018), albeit without incorporating power iterations. The power iteration concept is fundamental for modern high-performance randomized matrix decompositions, but to the best of our knowledge, has not been applied in the context of tensors. Without power iterations, the performance of randomized algorithms based on random projections can suffer considerably in the presence of white noise. We combine power iterations with oversampling to further control the error of the decomposition.

Embedding the CP decomposition into this probabilistic framework allows us to achieve significant computational savings, while producing near-optimal approximation quality. For motivation, Figure 1 shows the computational savings of a rank $k=20$ approximation for varying tensor dimensions. Here we compare the ALS and BCD solver for computing the deterministic and randomized CP decomposition. The computational advantage of the randomized algorithms is significant, while sacrificing a negligible amount of accuracy.

Organization.

The paper is organized as follows: Section 2 briefly reviews the CP decomposition and randomized matrix decompositions. Section 3 introduces the randomized CP tensor decomposition algorithm. Section 4 presents the evaluation of the computational performance, and examples. Finally, Section 5 summarizes the research findings.

2 Background

Ideas for multi-way factor analysis emerged in the 1920s with the formulation of the polyadic decomposition by Hitchcock (1927). However, the polyadic tensor decomposition only achieved popularity much later in the 1970s with the canonical decomposition (CANDECOMP) in psychometrics, proposed by Carroll and Chang (1970). Concurrently, the method of parallel factors (PARAFAC) was introduced in chemometrics by Harshman (1970). Hence, this method became known as the CP (CANDECOMP/PARAFAC) decomposition. However, in the past, computation was severely inhibited by available computational power. Today, tensor decompositions enjoy increasing popularity, yet runtime bottlenecks persist.

2.1 Notation

Scalars are denoted by lower case letters $x$, vectors as bold lower case letters $\mathbf{x}$, and matrices by bold capitals $\mathbf{X}$. Tensors are denoted by calligraphic letters $\bm{\mathcal{X}}$. The mode-$n$ unfolding of a tensor is expressed as $\bm{\mathcal{X}}_{(n)}$, while the mode-$n$ folding of a matrix is defined as $\mathbf{X}_{(n)}$. The vector outer product, the Kronecker product, the Khatri-Rao product and the Hadamard product are denoted by $\circ$, $\otimes$, $\odot$, and $\ast$, respectively. The inner product of two tensors is expressed as $\langle\cdot,\cdot\rangle$, and $\|\cdot\|_{F}$ denotes the Frobenius norm for both matrices and tensors.

2.2 CP Decomposition

The CP decomposition is the tensor equivalent of the SVD since it approximates a tensor by a sum of rank-one tensors. Here, tensor rank is defined as the smallest sum of rank-one tensors required to generate the tensor (Hitchcock, 1927). The CP decomposition approximates these rank-one tensors. Given a third order tensor $\bm{\mathcal{X}}\in\mathbb{R}^{I\times J\times K}$, the rank-$R$ CP decomposition is expressed as

[TABLE]

where $\circ$ denotes the outer product. Specifically, each rank-one tensor is formulated as the outer product of the rank-one components $\mathbf{a}_{r}\in\mathbb{R}^{I}$, $\mathbf{b}_{r}\in\mathbb{R}^{J}$, and $\mathbf{c}_{r}\in\mathbb{R}^{K}$. Components are often constrained to unit length with the weights absorbed into the vector $\bm{\lambda}=[\lambda_{1},...,\lambda_{R}]\in\mathbb{R}^{R}$. Equation (1) can then be re-expressed as (See Fig. 2)

[TABLE]

More compactly the components can be expressed as factor matrices, i.e., $\mathbf{A}=[\mathbf{a}_{1},\mathbf{a}_{2},...,\mathbf{a}_{R}]$, $\mathbf{B}=[\mathbf{b}_{1},\mathbf{b}_{2},...,\mathbf{b}_{R}]$, and $\mathbf{C}=[\mathbf{c}_{1},\mathbf{c}_{2},...,\mathbf{c}_{R}]$ . Using the Kruskal operator as defined by Kolda and Bader (2009), (2) can be more compactly expressed as

[TABLE]

2.3 Randomized Matrix Algorithms

The efficient computation of low rank matrix approximations is a ubiquitous problem in machine learning and data mining. Randomized matrix algorithms have been demonstrated to be highly competitive and robust when compared to traditional deterministic methods. Randomized algorithms aim to construct a smaller matrix (henceforth called sketch) designed to capture the essential information of the source matrix (Mahoney, 2011; Drineas and Mahoney, 2016). The sketch can then be used for various learning tasks. There exist several strategies for obtaining such a sketch, with random projections being the most robust off-the-shelf approach. Randomized algorithms have been in particular studied for computing the near-optimal low-rank SVD (Frieze et al., 2004; Liberty et al., 2007; Woolfe et al., 2008; Martinsson et al., 2011). Following the seminal work by Halko et al. (2011), a randomized algorithm computes the low-rank matrix approximation

[TABLE]

where the target rank is denoted as $k$, and assumed to be $k\ll\text{min}\{m,n\}$. Intuitively, we construct a sketch $\mathbf{Y}\in\mathbb{R}^{m\times k}$ by forming a random linear weighted combination of the columns of the source matrix $\mathbf{A}$. More concisely, we construct the sketch as

[TABLE]

where $\mathbf{\Omega}\in\mathbb{R}^{m\times k}$ is a random test matrix with entries drawn from, for example, a Gaussian distribution. If $\mathbf{A}$ has exact rank $k$, then the sketch $\mathbf{Y}$ spans with high probability a basis for the column space of the source matrix. However, most data matrices do not feature exact rank (i.e., the singular values $\{\sigma_{i}\}_{i=k+1}^{n}$ are non-zero). Thus, instead of just using $k$ samples, it is favorable to construct a slightly oversampled sketch $\mathbf{Y}\in\mathbb{R}^{m\times l}$ which has $l=k+p$ columns, were $p$ denotes the number of additional samples. In most situations small values $p=\{10,20\}$ are sufficient to obtain a good basis that is comparable to the best possible basis (Martinsson, 2016). An orthonormal basis $\mathbf{Q}\in\mathbb{R}^{m\times l}$ is then obtained via the QR-decomposition $\mathbf{Y}=\mathbf{Q}\mathbf{R}$, such that

[TABLE]

is satisfied. Finally, the source matrix $\mathbf{A}$ is projected to low-dimensional space

[TABLE]

where $\mathbf{B}\in\mathbb{R}^{l\times n}$. (Note, that Eq. 4 requires a second pass over the source matrix.) The matrix $\mathbf{B}$ can then be used to efficiently compute the matrix decomposition of interest, e.g., the SVD. The approximation error can be controlled by a combination of oversampling and power iteration (Rokhlin et al., 2009; Halko et al., 2011; Gu, 2015).

Randomized matrix algorithms are not only pass efficient, but they also have the ability to exploit modern computational parallelized and distributed computing architectures. Implementations in MATLAB, C and R are provided by Szlam et al. (2014), Voronin and Martinsson (2015), and Erichson et al. (2016).

3 Randomized CP Decomposition

Given a third order tensor $\bm{\mathcal{X}}\in\mathbb{R}^{I\times J\times K}$, the objective of the CP decomposition is to find a set of $R$ normalized rank-one tensors $\{\mathbf{a}_{r}\circ\mathbf{b}_{r}\circ\mathbf{c}_{r}\}_{r=1}^{R}$ which best approximates $\bm{\mathcal{X}}$, i.e., minimizes the Frobenius norm

[TABLE]

The challenge is that this problem is highly nonconvex and unlike PCA there is no closed-form solution. Solution methods for this optimization problem therefore rely on iterative schemes (e.g., alternating least squares). These solvers are not guaranteed to find a global minimum, yet for many practical problems, they can find high-quality solutions. Iterative schemes, however, can be computationally demanding if the dimensions of $\bm{\mathcal{X}}$ are large.

Fortunately, only the column spaces of the modes are of importance for obtaining the factor matrices $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$, rather then the individual columns of the mode matricizations $\bm{\mathcal{X}}_{(1)}$, $\bm{\mathcal{X}}_{(2)}$, $\bm{\mathcal{X}}_{(3)}$ of the tensor $\bm{\mathcal{X}}$. This is because the CP decomposition learns the components based on proportional variations in inter-point distances between the components. Therefore, a compressed tensor $\bm{\mathcal{B}}\in\mathbb{R}^{k\times k\times k}$ must preserve pairwise Euclidean distances, where $k\geq R$. This in turn requires that column spaces, and thus pairwise distances, are approximately preserved by compression - this can be achieved by generalizing the concepts of randomized matrix algorithms to tensors. We build upon the methods introduced by Martinsson et al. (2011) and Halko et al. (2011), as well as related work on randomized tensors by Drineas and Mahoney (2007), who proposed a randomized algorithm based on random column selection.

Figure 3 shows the schematic of the randomized CP decomposition architecture. Our approach proceeds in two stages. The first stage, detailed in Section 3.1, applies random projections with power iteration to obtain a compressed tensor $\mathcal{B}$, with expressivity analysis in Section 3.1.1. In Section 3.2, we describe two algorithms for performing CP on the compressed tensor $\mathcal{B}$ and approximating the original CP factor matrices. Section 3.3 provides additional details about our implementation in Python.

3.1 Randomized Tensor Algorithm

The aim is to use randomization as a computational resource to efficiently build a suitable basis that captures the action of the tensor $\bm{\mathcal{X}}$. Assuming an $N$-way tensor $\bm{\mathcal{X}}\in\mathbb{R}^{I_{1}\times\cdots\times I_{N}}$, the aim is to obtain a smaller compressed tensor $\bm{\mathcal{B}}\in\mathbb{R}^{k\times\cdots\times k}$, so that its $N$ tensor modes capture the action of the input tensor modes. Hence, we seek a natural basis in the form of a set of orthonormal matrices $\{\mathbf{Q}_{n}\in\mathbb{R}^{I_{n}\times R_{n}}\}_{n=1}^{N}$, so that

[TABLE]

Here the operator $\times_{n}$ denotes tensor-matrix multiplication defined as follows.

Definition 1

The $n$-mode matrix product $\bm{\mathcal{X}}\times_{n}\mathbf{Q}_{n}\mathbf{Q}_{n}^{\top}$ multiplies a tensor by the matrices $\mathbf{Q}_{n}\mathbf{Q}_{n}^{\top}$ in mode $n$, i.e., each mode-$n$ fiber is multiplied by $\mathbf{Q}_{n}\mathbf{Q}_{n}^{\top}$

[TABLE]

Given a fixed target rank $k$, these basis matrices can be efficiently obtained using a randomized algorithm. First, the approximate basis for the $n$-th tensor mode is obtained by drawing $k$ random vectors $\mathbf{\omega}_{1},\dots,\mathbf{\omega}_{k}\in\mathbb{R}^{\prod_{i\neq n}{I_{i}}}$ from a Gaussian distribution. (Note that if $N$ is large, it might be favorable to draw the entries of $\omega$ from a quasi random sequence, e.g., Halton or Sobol sequence. These so-called quasi-random numbers are known to be less correlated in high-dimensions.) Stacked together, the $k$ random vectors $\omega$ form the random test matrix $\mathbf{\Omega}_{n}\in\mathbb{R}^{\prod_{i\neq n}{I_{i}}\times k}$ used to sample the column space of $\bm{\mathcal{X}}_{(n)}\in\mathbb{R}^{I_{n}\times\prod_{i\neq n}{I_{i}}}$ as follows

[TABLE]

where $\mathbf{Y}_{n}\in\mathbb{R}^{I_{n}\times k}$ is the sketch. The sketch serves as an approximate basis for the range of the $n$-th tensor mode. Probability theory guarantees that the set of random vectors $\{\mathbf{\omega}_{i}\}_{i=1}^{k}$ are linearly independent with high probability. Hence, the corresponding random projections $\mathbf{y}_{1},\dots,\mathbf{y}_{k}$ efficiently sample the range of a rank deficient tensor mode $\bm{\mathcal{X}}_{(n)}$. The economic QR decomposition of the sketch $\mathbf{Y}_{n}=\mathbf{Q}_{n}\mathbf{R}_{n}$ is then used to obtain a natural basis, so that $\mathbf{Q}_{n}\in\mathbb{R}^{I_{n}\times k}$ is orthonormal and has the same column space as $\mathbf{Y}_{n}$. The final step restricts the tensor mode to this low-dimensional subspace

[TABLE]

Thus, after $N$ iterations a compressed tensor ${\bm{\mathcal{B}}}$ and a set of orthonormal matrices is obtained. Since this is an iterative algorithm, we set $\bm{\mathcal{X}}\leftarrow\mathbf{\bm{\mathcal{B}}}_{n}$ after each iteration. The number of columns of the basis matrices form a trade-off between accuracy and computational performance. We aim to use as few columns as possible, yet allow an accurate approximation of the input tensor. Assuming that the tensor $\bm{\mathcal{X}}$ exhibits low-rank structure, or equivalently, the rank $R$ is much smaller than the ambient dimensions of the tensor, the basis matrices will be an efficient representation. In practice, the compression performance is improved through using oversampling, i.e., drawing $l=k+p$ random vectors where $k$ is the target rank and $p$ the oversampling parameter.

The randomized algorithm as presented requires that the mode-$n$ unfolding of the tensor has a rapidly decaying spectrum in order to achieve good performance. However, this assumption is often not suitable, and the spectrum exhibits slow decay if the tensor is compressed several times. To overcome this issue, the algorithm’s performance can be substantially improved using power iterations (Rokhlin et al., 2009; Halko et al., 2011; Gu, 2015). Power iterations turn a slowly decaying spectrum into a rapidly decaying one by taking powers of the tensor modes. Thus, instead of sampling $\bm{\mathcal{X}}_{(n)}$ we sample from the following tensor mode

[TABLE]

where $q$ denotes a small integer. This power operation enforces that the singular values $\sigma_{j}$ of $\bm{\mathcal{X}}_{(n)}^{q}$ are $\{\!\sigma_{j}^{2q+1}\!\}_{j}$. Instead of using (7), an improved sketch is computed

[TABLE]

However, if (9) is implemented in this form the basis may be distorted due to round-off errors. Therefore in practice, normalized subspace iterations are used to form the sketch, meaning that the sketch is orthornormalized between each power iteration in order to stabilize the algorithm. For implementation, see Voronin and Martinsson (2015) and Szlam et al. (2014).

The combination of oversampling and additional power iterations can be used to control the trade-off between approximation quality and computational efficiency of the randomized tensor algorithm. Our results, for example, show that just $q=2$ subspace iterations and an oversampling parameter of about $p=10$ achieves near-optimal results. Algorithm 1 summarizes the computational steps.

3.1.1 Expressivity analysis

The average behavior of the randomized tensor algorithm is characterized using the expected residual error

[TABLE]

where $\hat{\bm{\mathcal{X}}}=\bm{\mathcal{X}}\times_{1}\mathbf{Q}_{1}\mathbf{Q}_{1}^{\top}\times_{2}\cdots\times_{N}\mathbf{Q}_{N}\mathbf{Q}_{N}^{\top}$. Theorem 1 generalizes the matrix version of Theorem 10.5 formulated by Halko et al. (2011) to the tensor case.

Theorem 1** (Expected Frobenius error)**

Consider a low-rank real $N$-way tensor $\bm{\mathcal{X}}\in\mathbb{R}^{I_{1}\times\cdots\times I_{N}}$. Then the expected approximation error, given a target rank $k\geq 2$ and an oversampling parameter $p\geq 2$ for each mode, is

[TABLE]

where $\sigma_{nj}$ denotes the $j$ singular value of the mode-$n$ unfolding of the source tensor $\bm{\mathcal{X}}$.

For the proof see Appendix A. Intuitively, the projection of each tensor mode onto a low-dimensional space introduces an additional residual. This is expressed by the double sum on the right hand side. If the low-rank approximation captures the column space of each mode accurately, then the singular values $j>k$ for each mode $n$ are small. Moreover, the error can be improved using the oversampling parameter. The computation of additional power (subspace) iterations can improve the error further. This result again follows by generalizing the results of Halko et al. (2011) to tensors. Sharper performance bounds for both oversampling and additional power iterations can be derived following, for instance, the results by Witten and Candes (2015).

3.2 Optimization Strategies

There are several optimization strategies for minimizing the objective function defined in (12), of which we consider, alternating least squares (ALS) and block coordinate descent (BCD). Both methods are suitable to operate on a compressed tensor $\bm{\mathcal{B}}\in\mathbb{R}^{k\times\cdots\times k}$, where $k\geq R$. The optimization problem (5) is reformulated

[TABLE]

Once the compressed factor matrices $\mathbf{\tilde{A}}\in\mathbb{R}^{k\times R}$, $\mathbf{\tilde{B}}\in\mathbb{R}^{k\times R}$, $\mathbf{\tilde{C}}\in\mathbb{R}^{k\times R}$ are estimated, the full factor matrices can be recovered

[TABLE]

where $\mathbf{Q}_{1}\!\in\!\mathbb{R}^{I\times k}$, $\mathbf{Q}_{2}\!\in\!\mathbb{R}^{J\times k}$, $\mathbf{Q}_{3}\!\in\!\mathbb{R}^{K\times k}$ denote the orthonormal basis matrices. For simplicity we focus on third order tensors, but the result generalizes to $N$-way tensors.

3.2.1 ALS Algorithm

Due to its simplicity and efficiency, ALS is the most popular method for computing the CP decomposition (Comon et al., 2009; Kolda and Bader, 2009). We note that the optimization (12) is equivalent to

[TABLE]

with respect to the factor matrices $\bf{\tilde{A},\tilde{B}}$ and $\bf{\tilde{C}}$. The tensor $\bm{\mathcal{B}}$ can further be expressed in matricized form

[TABLE]

where $\odot$ denotes the Khatri-Rao product. The optimization problem in this form is non-convex. However, an estimate for the factor matrices can be obtained using the least-squares method as follows. The ALS algorithm updates one component, while holding the other two components fixed, in an alternating fashion until convergence. It iterates over the following subproblems

[TABLE]

Each step therefore involves a least-squares problem which can be solved using the Khatri-Rao product pseudo-inverse. Algorithm 2 summarizes the computational steps.

Definition 2

The Khatri-Rao product pseudo-inverse is defined as

[TABLE]

where the operator $*$ denotes the Hadamard product, i.e., the elementwise multiplication of two equal sized matrices.

There exist few general convergence guarantees for the ALS algorithm (Uschmajew, 2012; Wang and Chu, 2014). Moreover, the final solution tends to depend on the initial guess $\bf{\tilde{A}}^{0},\bf{\tilde{B}}^{0}$ and $\bf{\tilde{C}}^{0}$. A standard initial guess uses the eigenvectors of $\bm{\mathcal{B}}_{(1)}\bm{\mathcal{B}}_{(1)}^{\top}$, $\bm{\mathcal{B}}_{(2)}\bm{\mathcal{B}}_{(2)}^{\top}$, $\bm{\mathcal{B}}_{(3)}\bm{\mathcal{B}}_{(3)}^{\top}$ (Bader and Kolda, 2015). Further, it is important to note that normalization of the factor matrices is necessary after each iteration to achieve good convergence. This prevents singularities of the Khatri-Rao product pseudo-inverse Kolda and Bader (2009). The algorithm can be further improved by reformulating the above subproblems as regularized least-squares problems, see for instance Li et al. (2013) for technical details and convergence results. The structure imposed by the ALS algorithm on the factor matrices permits the formulation of non-negative, or sparsity-constrained tensor decompositions (Cichocki et al., 2009).

3.2.2 BCD Algorithm

While ALS is the most popular algorithm for computing the CP decomposition, many alternative algorithms have been developed. One alternate approach is based on block coordinate descent, BCD (Xu and Yin, 2013). Cichocki and Phan (2009) first proposed this approach for computing nonnegative tensor factorizations. The BCD algorithm is based on the idea of successive rank-one deflation. Unlike ALS, which updates the entire factor matrix at each step, BCD computes the rank-1 tensors in a hierarchical fashion. Therefore, the algorithm treats each component $\mathbf{a}_{r},\mathbf{b}_{r},\mathbf{c}_{r}$ as a block. First, the most correlated rank-1 tensor is computed; then the second most correlated rank-1 tensor is learned on the residual tensor, and so on. Assuming that $\tilde{R}=r-1$ components have been computed, then the $r$-th compressed residual tensor $\bm{\mathcal{Y}}_{\text{res}}$ is defined

[TABLE]

The algorithm then iterates over the following subproblems

[TABLE]

Note that the computation can be more efficiently evaluated without explicitly constructing the residual tensor $\bm{\mathcal{Y}}_{\text{res}}$ (Kim et al., 2014). Algorithm 3 summarizes the computation.

3.3 Implementation Details

The algorithms we present are implemented in the programming language Python, using numerical linear algebra tools provided by the SciPy (Open Source Library of Scientific Tools) package (Jones et al., 2001). SciPy provides MKL (Math Kernel Library) accelerated high performance implementations of BLAS and LAPACK routines. Thus, all linear algebra operations are threaded and highly optimized on Intel processors. The implementation of the CP decomposition follows the MATLAB Tensor Toolbox implementation (Bader and Kolda, 2015). This implementation normalizes the components after each step to achieve better convergence. Furthermore, we use eigenvectors (see above) to initialize the factor matrices. Interestingly, randomly initialized factor matrices have the ability to achieve slightly better approximation errors, but re-running the algorithms several times with different random seeds can display significant variance in the results. Hence only the former approach is used for initialization. We note that the randomized algorithm introduces some randomness and slight variations into the CP decompositions as well. However, randomization can also act as an implicit regularization on the CP decomposition (Mahoney, 2011), meaning that the results of the randomized algorithm can be in some cases even ‘better’ than the results of the corresponding deterministic implementation. Finally, the convergence criterion is defined as the change in fit, following Bader and Kolda (2015). The algorithm therefore stops when the improvement of the fit $\rho$ is less then a predefined threshold, where the fit is computed using $\rho=1-({\|\bm{\mathcal{X}}\|_{F}^{2}+\|\hat{\bm{\mathcal{X}}}\|_{F}^{2}-2\cdot\langle\hat{\bm{\mathcal{X}}},\bm{\mathcal{X}}\rangle})/{\|\bm{\mathcal{X}}\|_{F}^{2}}$.

4 Numerical Results

The randomized CP algorithm is evaluated on a number of examples where the near optimal approximation of massive tensors can be achieved in a fraction of the time using the randomized algorithm. Approximation accuracy is computed with the relative error ${\|\bm{\mathcal{X}}-\hat{\bm{\mathcal{X}}}\|_{F}}/{\|\bm{\mathcal{X}}\|_{F}}$, where $\hat{\bm{\mathcal{X}}}$ denotes the approximated tensor.

4.1 Computational Performance

The robustness of the randomized CP algorithm is first assessed on random low-rank tensors. Here we illustrate the approximation quality in the presence of additive white noise. Figure 4 shows the average relative errors over $100$ runs for varying signal-to-noise ratios (SNR). In the case of high SNRs, all algorithms converge towards the same relative error. However, at excessive levels of noise (i.e., SNR$<4$) the deterministic CP algorithms exhibit small gains in accuracy over the randomized algorithms using $q=2$ power iterations, with which both the ALS and BCD algorithm show the same performance. The performance of the randomized algorithm without power iterations ($q=0$) is, however, poor, and stresses the importance of the power operation for real applications. The oversampling parameter for the randomized algorithms is set to $p=10$. Increasing $p$ can further improve accuracy.

Next, the reconstruction errors and runtimes for tensors of varying dimensions are compared. Figure 5 shows the average evaluation results over $100$ runs for random low-rank tensors of different dimensions, and for varying target ranks $k$. The randomized algorithms achieve near optimal approximation accuracy while demonstrating substantial computational savings. Interestingly, the relative error achieved by the BCD decreases sharply by about one order of magnitude when the target rank $k$ matches the actual tensor rank (here $R=50$).

The computational advantage becomes more pronounced with increasing tensor dimensions, and as the number of iterations required for convergence increase. Using random tensors as presented here, all algorithms rapidly converge after about $4$ to $6$ iterations. However, it is evident that the computational cost per iteration of the randomized algorithm is substantially lower. Thus, the computational savings can be even more substantial in real world applications that may require several hundred iterations to converge. Overall, the ALS algorithm is computationally more efficient than BCD, i.e., the deterministic ALS algorithm is faster than the BCD by nearly one order of magnitude, while the randomized algorithms exhibit similar computational timings.

Similar performance results are achieved for higher order tensors. Figure 5(c) shows the computational performance for a $4$-way tensor of dimension $100\times 100\times 100\times 100$. Again, the randomized algorithms achieve speedups of 1-2 orders of magnitude, while attaining good approximation errors.

4.2 Examples

The examples demonstrate the advantages of the randomized CP decomposition. The first is a multiscale toy video example, and the second is a simulated flow field behind a stationary cylinder. Due to the better and more natural interpretability of the BCD algorithm, only this algorithm is considered in subsequent sections.

4.2.1 Multiscale Toy Video Example

The approximation of the underlying spatial modes and temporal dynamics of a system is a common problem in signal processing. In the following, we consider a toy example presenting multiscale intermittent dynamics in the time direction. The data consists of four Gaussian modes, each undergoing different frequencies of intermittent oscillation, for $215$ times steps in the temporal direction on a two-dimensional spatial grid ($200\times 200$). The resulting tensor is of dimension $200\times 200\times 215$. Figure 6 shows the corresponding modes and the time dynamics. This problem becomes even more challenging when the underlying structure needs to be reconstructed from noisy measurements. Traditional matrix decomposition techniques such as the SVD are known to face difficulties capturing such intermittent, multiscale structure.

Figure 7 displays the decomposition results of the noisy (SNR=$2$) toy video for both the randomized CP decomposition and the SVD. The first subplot shows the results of a rank $k=4$ approximation computed using the rCP algorithm with $q=2$ power iterations, and a small oversampling parameter $p=10$. The method faithfully captures the underlying spatial modes and the time dynamics. For illustration, the second subplot shows the decomposition results without additional power iterations. It can clearly be seen that this approach introduces distinct artifacts, and the approximation quality is relatively poor. The SVD, shown in the last subplot, demonstrates poor performance at separating the modes and mixes the spatiotemporal dynamics of modes 2 & 3.

Table 1 further quantifies the observed results. Interestingly, the relative error using the randomized algorithm with $q=2$ power iterations is slightly better than the deterministic algorithm. This is due to the intrinsic regularizing behavior of randomization. However, the reconstruction error without power iterations is large, as is the error resulting from the SVD. The achieved speedup of randomized CP is substantial, with a speedup factor of about $15$.

A Note on Compression.

The CP decomposition provides a more parsimonious representation of the data. Comparing the compression ratios between the CP decomposition and SVD illustrates the difference. For a rank $R=4$ tensor of dimension $100\times 100\times 100$, the compression ratios are

[TABLE]

Note, the SVD requires the tensor to be reshaped in some direction. The comparison illustrates the striking difference between compression ratios. It is evident that the CP decomposition requires computing many fewer coefficients in order to approximate the tensor. Thus, less storage is required to approximate the data. While the CP decomposition yields a parsimonious approximation that is more robust to noise, the advantage of the SVD is that data can be approximated with an accuracy as low as machine precision.

4.2.2 Flow behind a cylinder

Extracting the dominant coherent structures from fluid flows helps to better characterize them for modeling and control (Brunton and Noack, 2015). The workhorse algorithm in fluid dynamics and for model reduction is the SVD. However, fluid simulations generate high-resolution spatiotemporal grids of data which naturally manifest as tensors. In the following we examine the suitability of the CP decomposition for decomposing flow data, and compare the results to those of the SVD.

The fluid flow behind a cylinder, a canonical example in fluid dynamics, is presented here. The data are a time-series generated from a simulation of fluid vorticity behind a stationary cylinder (Colonius and Taira, 2008). The corresponding flow tensor has dimension $199\times 449\times 151$, consisting of 151 snapshots on a $449\times 199$ spatial grid. Figure 8 shows a few example snapshots of the fluid flow. The flow is characterized by a periodically shedding wake structure and is inherently low-rank in the absence of noise. The characteristic frequencies of flow oscillations occur in pairs, reflecting the complex-conjugate pairs of eigenvalues corresponding to harmonics in the temporal direction.

Results in Absence of White Noise.

Figure 9 shows both the approximated spatial modes and the temporal dynamics for the randomized CP decomposition and the SVD. Observe that both methods extract similar patterns, although unlike SVD, rCP spatial modes are sparse on the domain. The reason is two-fold: (i) CP does not impose orthogonality constraints (which in general reveal dense structure), and (ii) CP imposes rank-one outer product structure in the $x$ and $y$ directions via the columns of the factor matrices. In doing so, CP isolates the contributions of single wavenumbers (spatial frequencies) to the steady-state vortex shedding regime. These are commonly analyzed to study energy transfer between scales in complex flows and turbulence (Meneveau and Katz, 2000; Sharma and McKeon, 2013), hence CP can help reveal new physically meaningful insights. In particular, randomization enables tensor decomposition of high-resolution vector fields that otherwise may not be tractable using deterministic methods. Thus rCP provides a novel decomposition of multimodal fields for downstream tasks in model reduction, pattern extraction and control.

Note, that for a fixed target rank of $k=30$ across all methods, the SVD achieves a substantially lower reconstruction error (see Table 2). However, the compression ratios for the CP and SVD methods are $c_{CP}\approx 562.17$ and $c_{SVD}\approx 5.02$, i.e., the CP compresses the data nearly two orders of magnitude more.

Results in Presence of White Noise.

Next, the analysis of the same flow is repeated in the presence of additive white noise. While this is not of concern when dealing with flow simulations, it is realistic when dealing with flows obtained from measurement. We chose a signal-to-noise ratio of 2 to demonstrate the robustness of the CP decomposition to noise.

Figure 10 shows again the corresponding dominant spatial modes and temporal dynamics. Both the SVD and the CP decomposition faithfully capture the temporal dynamics. However, comparing the modes of the SVD to Figure 9, it is apparent that the spatial modes contain a large amount of noise. The spatial modes revealed by the CP decomposition provide a significantly better approximation. Again, it is crucial to use power iterations to achieve a good approximation quality (see Table 3). By inspection, the relative reconstruction error using the SVD is poor compared to the error achieved using the CP decomposition. Here, we show the error for a rank $k=30$ and $k=6$ approximation. The target rank was determined using the optimal hard threshold for singular values (Gavish and Donoho, 2014).

The CP decomposition overcomes this disadvantage, and is able to approximate the first $k=30$ modes with only a slight loss of accuracy. Note that here the randomized CP decomposition performs better than the deterministic algorithm. We assume that this is due to the favorable intrinsic regularization effect of randomized methods.

5 Conclusion

The emergence of massive tensors require efficient algorithms for obtaining tensor decompositions. To address this challenge, we have presented a randomized algorithm which substantially reduces the computational demands of the CP decomposition. Indeed, randomized algorithms have established themselves as highly competitive methods for computing traditional matrix decompositions. A key advantage of the randomized algorithm is that modern computational architectures are fully exploited. Thus, the algorithm benefits substantially from multithreading in a multi-core processor. In contrast to previously proposed high-performance tensor algorithms which are based on computational concepts such as distributed computing, our proposed randomized algorithm provides substantial computational speedups even on standard desktop computers. Our proposed algorithm achieves these speedups by reducing the computational costs per iteration, which enables the user to decompose real world examples that typically require a large number of iterations to converge.

In addition to computational savings, the randomized CP decomposition demonstrates outstanding performance on several examples using artificial and real-world data, including decompositions of high-resolution flow fields that may not be tractable with deterministic methods. Moreover, our experiments show that the power iteration concept is crucial in order to achieve a robust tensor decomposition. Thus, our algorithm has a practical advantage over previous randomized tensor algorithms, at a slightly increased computational cost due to additional power iterations.

Acknowledgments

We would like to thank Alex Williams, and Michael W. Mahoney for insightful discussion on tensor decompositions and randomized numerical linear algebra. Further, we would like to express our gratitude to the two anonymous reviewers for their valuable feedback, which helped us greatly improve the manuscript. NBE would like to acknowledge DARPA and NSF for providing partial support of this work. KM acknowledges support from NSF MSPRF Award 1803663. JNK acknowledges support from the Air Force Office of Scientific Research (AFOSR) grant FA9550-17-1-0329. SLB acknowledges funding support from the Air Force Office of Scientific Research (AFOSR) grant FA9550-18-1-0200.

Data availability statement

The data that support the findings of this study are available upon request.

A Proof of Theorem 1

In the following, we sketch a proof for Theorem 1, which yields an upper bound for the approximate basis for the range of a tensor. To assess the quality of the basis matrices $\{\mathbf{Q}_{n}\}_{n=1}^{N}$, we first show that the problem can be expressed as a sum of subproblems. Defining the residual error

[TABLE]

Note that the Frobenius norm of a tensor and its matricized forms are equivalent. Defining the orthogonal projector $\mathbf{P}_{n}\equiv\mathbf{Q}_{n}\mathbf{Q}_{n}^{\top}$, we can reformulate (21) compactly as

[TABLE]

**Proof ** Assuming that $\mathbf{P}_{n}$ yields an exact projection onto the column space of the matrix $\mathbf{Q}_{n}$, we need to show first that the error can be expressed as a sum of the errors of the $n$ projections

[TABLE]

where $\mathbf{I}$ denotes the identity matrix. Following Drineas and Mahoney (2007), let us add and subtract the term $\bm{\mathcal{X}}\times_{N}\mathbf{P}_{N}$ in Equation (22) so that we obtain

[TABLE]

The bound (25) follows from the triangular inequality for a norm. Next, the common term $\mathbf{P}_{N}$ is factored out in Equation (26). Then, the bound (27) follows from the properties of orthogonal projectors. This is because the $range(\bm{\mathcal{X}}\times_{1}\mathbf{P}_{1}\times_{1}\cdots\times_{N-1}\mathbf{P}_{N-1})\subset range(\bm{\mathcal{X}}\times_{1}\mathbf{P}_{1}\times_{1}\cdots\times_{N}\mathbf{P}_{N})$, and then it holds that $\|\bm{\mathcal{X}}-\bm{\mathcal{X}}\times_{1}\mathbf{P}_{1}\times_{1}\cdots\times_{N}\mathbf{P}_{N}\|_{F}\leq\|\bm{\mathcal{X}}-\bm{\mathcal{X}}\times_{1}\mathbf{P}_{1}\times_{1}\cdots\times_{N-1}\mathbf{P}_{N-1}\|_{F}$. See Proposition 8.5 by Halko et al. (2011) for a proof using matrices. Subsequently the residual error $\bm{\mathcal{E}}_{N-1}$ can be bounded

[TABLE]

From this inequality, Equation (23) follows. We take the expectation of Equation (23)

[TABLE]

Recalling that Theorem 10.5 formulated by Halko et al. (2011) states the following expected approximation error (formulated here using tensor notation)

[TABLE]

assuming that the sketch in Equation (7) is constructed using a standard Gaussian matrix $\mathbf{\Omega}$. Here $\sigma_{j}$ denotes the singular values of the matricized tensor $\bm{\mathcal{X}}_{(n)}$ greater then the chosen target rank $k$. Combining Equations (29) and (30) then yields the results of the theorem (11).

Figures 11 evaluates the theoretical upper bound over 100 runs for both a third and fourth order random low-rank $R=25$ tensor. Here, we use a fixed oversampling parameter $p=2$. The results show that the empirical error is faithfully bounded by the theoretical upper bound for varying target ranks.