Sampling-Based Decomposition Algorithms for Arbitrary Tensor Networks

Tensor decomposition has emerged as an important tool in data mining and machine learning [24, 4, 5, 11]. Applications in data mining include network analysis, web mining, topic modeling and recommendation systems. Tensor decomposition is used widely in machine learning for things like parameter reduction in neural networks, understanding of deep neural network expressiveness, supervised learning and feature extraction.

Due to the multidimensional nature of tensors, they are inherently plagued by the curse of dimensionality. For example, representing a tensor $\bm{\mathscr{X}}\in\mathbb{R}^{I\times\cdots\times I}$ with $N$ modes requires $I^{N}$ numbers. This exponential dependence on $N$ makes its way into algorithms for computing tensor decompositions. Alternating least squares (ALS) is arguably the most popular and successful approach for computing a wide range of tensor decompositions. When decomposing a tensor $\bm{\mathscr{X}}$, each iteration of ALS involves solving a sequence of least squares problems for which the entries of $\bm{\mathscr{X}}$ feature as the dependent variables. The per-iteration cost of ALS therefore naturally inherits the exponential dependence on $N$.

A large number of papers have sought to reduce the cost of tensor decomposition by leveraging techniques from randomized numerical linear algebra. One particularly interesting line of work [3, 14, 9, 19, 18] seeks to construct ALS algorithms with a per-iteration cost which is sublinear in the number of tensor entries, i.e., $o(I^{N})$. Since any algorithm considering all entries of $\bm{\mathscr{X}}$ immediately incurs a cost of $\Omega(I^{N})$, these works have all resorted to sampling-based techniques. More precisely, they all sample the ALS subproblems according to the leverage score distribution (or an approximation thereof). This is done efficiently by taking advantage of the special structure of the design matrices in the ALS subproblems for tensor decomposition.

The previous works discussed above develop methods for specific tensor decompositions: The CP decomposition in [3, 14, 18], Tucker decomposition in [9], and the tensor ring decomposition in [19, 18]. In this paper, we consider all decompositions that can be expressed in tensor network (TN) format. The TN format allows for a very wide range of decompositions, including the CP, Tucker, tensor train and tensor ring decompositions. The following summarizes our contributions:

• •

  We first show how to efficiently sample rows of any tall-and-skinny matrix which is in TN format according to the exact leverage score distribution.

• •

  We then show how this sampling technique can be used to yield ALS algorithms with a per-iteration cost which is input sublinear for all TN decompositions that satisfy certain mild assumptions.

The decomposition framework we present builds on the work in [18]. That paper is notable since it provided the first sampling-based ALS methods for CP and tensor ring decomposition with a per-iteration cost depending polynomially on $N$; the dependence in earlier works was exponential [3, 14, 19]. We are able to substantially simplify the scheme in [18] by entirely avoiding the complicated recursive sketching procedure that it relies on. This makes implementation easier and is also what ensures that the sampling is done according to the exact leverage score distribution rather than an approximation of it. The simplification is also what paves the way for generalization to arbitrary TN formats.

Cheng et al. [3] develop the first ALS method for CP decomposition with an input sublinear per-iteration cost. Their method uses a mixture of leverage score and row-norm sampling applied to the matricized-tensor-times-Khatri–Rao product (MTTKRP) which arises as a key computational kernel in CP decomposition. Larsen and Kolda [14] use leverage score sampling to reduce the size of each ALS subproblem. Their method has improved theoretical guarantees compared to those in [3] as well as improved scalability due to various practical improvements. Malik and Becker [19] propose a sampling-based ALS method for tensor ring decomposition which also uses leverage score sampling to reduce the size of the ALS subproblems.

All three papers [3, 14, 19] require a number of samples which scales exponentially with the number of input tensor modes, $N$, for performance guarantees to hold. This translates into a per-iteration cost which is $\Omega(R^{N+1})$ for the CP decomposition [3, 14] and $\Omega(R^{2N+2})$ for the tensor ring decomposition [19], where $R$ is the relevant notion of rank. For low-rank decompositions ($R\ll I$), the cost will be input sublinear despite this exponential cost, but for higher rank it might not be (recall that, unlike in matrix decomposition, we can have $R>I$ in tensor decomposition). Malik [18] develop methods for both CP and tensor ring decomposition which avoid this exponential dependence on $N$, instead improving it to a polynomial dependence. This is achieved by sampling from a distribution much closer to the exact leverage score distribution than the previous works [3, 14, 19] do. Since our work builds on and improves the scheme in [18], it also has a polynomial dependence on $N$ when used for CP and tensor ring decomposition (see Tables 1 and 2).

Fahrbach et al. [9] develop a method for efficient sampling of ridge regression problems involving Kronecker product design matrices. They use these techniques to achieve an efficient sampling-based ALS method for regularized Tucker decomposition.

Efficient sketching of structured matrices is an active research area with applications beyond tensor decomposition. The recent work by Ma and Solomonik [16] is particularly interesting since it proposes structured sketching operators for general TN matrices. These operators take the form of TNs made up of Gaussian tensors and are therefore dense. When used in ALS for tensor decomposition, their sketching operators therefore need to access all entries of the data tensor in each least squares solve which leads to a per-iteration cost which is at least linear in the input tensor size (and therefore exponential in $N$).

Due to space constraints, we have focused on previous works that develop sampling-based ALS methods here. There is a large number of works that develop randomized tensor decomposition methods using other techniques. For a comprehensive list of such works, see [18, Sec. 2]. For an overview of work on structured sketching, see [20, Sec. 1.2].

In this section, we propose an efficient method for sampling rows of a tall-and-skinny TN matrix $\bm{A}\in\mathbb{R}^{I\times R}$ according to its leverage scores. The method is useful when $I$ is so large that costs and storage on the order $\Omega(I)$ are too expensive, and $R$ is small enough that costs and storage on the order $O(R^{3})$ are affordable. In this setting it is infeasible to sample the rows of $\bm{A}$ directly via the probability distribution in Definition 3.2 since doing so would require computing the pseudoinverse of $\bm{A}$ which costs $O(IR^{2})$. Moreover, once that pseudoinverse is computed, we would be required to compute, store and sample according to a length-$I$ vector of probabilities.

In view of Definition 3.1, note that

where $\bm{\Phi}\stackrel{{\scriptstyle\text{\tiny{def}}}}{{=}}(\bm{A}^{\top}\bm{A})^{+}$. This formulation is the basis for our sampling method. As we show in Section 4.1, the matrix $\bm{\Phi}$ can be computed efficiently when $\bm{A}$ is a TN matrix. Once $\bm{\Phi}$ is computed, we can again leverage the TN structure of $\bm{A}$ to draw samples from its leverage score distribution via the formula in (4) without forming a length-$I$ probability vector. We describe how this is done in Section 4.2.

TNs can be used to express a very wide range of tensor decompositions, including the CP [10], Tucker [27], tensor train [23] and tensor ring decompositions [31].

Suppose $\bm{\mathscr{X}}\in\mathbb{R}^{I_{1}\times\cdots I_{N}}$ is a data tensor that we want to decompose into some TN format consisting of tensors $\bm{\mathscr{A}}^{(1)},\ldots,\bm{\mathscr{A}}^{(M)}$. We use $\operatorname{TN}(\bm{\mathscr{A}}^{(1)},\ldots,\bm{\mathscr{A}}^{(M)})\in\mathbb{R}^{I_{1}\times\cdots\times I_{N}}$ to denote this TN. This decomposition problem can be formulated as the optimization problem

where $\|\cdot\|_{\textup{F}}$ is a straightforward generalization of the matrix Frobenius norm to tensors (see Section A). In general, this is a difficult non-convex optimization problem which is hard to solve. A popular approach to finding an approximate solution to (8) is to iteratively update one tensor $\bm{\mathscr{A}}^{(m)}$ at a time via alternating least squares (ALS). Minimizing the objective in (8) with respect to a single tensor at a time turns it into a linear least squares problem:

where $\bm{X}^{(m)}$ and $\bm{A}^{(m)}$ are appropriate matricizations of $\bm{\mathscr{X}}$ and $\bm{\mathscr{A}}^{(m)}$, respectively, and $\bm{A}^{\neq m}$ is a TN matrix which depends on the tensors $\bm{\mathscr{A}}^{(1)},\ldots,\bm{\mathscr{A}}^{(m-1)},\bm{\mathscr{A}}^{(m+1)},\ldots,\bm{\mathscr{A}}^{(M)}$. For example, for the CP and Tucker decompositions, $\bm{A}^{\neq m}$ takes the form of a Khatri–Rao and Kronecker product matrix, respectively [12]. Figure 7 shows examples of three tensor decompositions (top plots) and examples of what a TN matrix $\bm{A}^{\neq m}$ corresponding to each decomposition looks like (bottom plots). Algorithm 1 provides a high-level ALS algorithm for computing any TN decomposition.

ALS is the “workhorse approach” to tensor decomposition [12]. ALS-based methods for various decompositions are implemented in many leading tensor software packages, including Tensor Toolbox [2] and Tensorlab [28] for Matlab and TensorLy [13] for Python. ALS typically works well, especially when the data tensor $\bm{\mathscr{X}}$ is not too large. However, as previously discussed, each iteration of ALS depends on all entries of $\bm{\mathscr{X}}$. If $\bm{\mathscr{X}}$ has $N$ modes each of dimension $I$, then $\bm{\mathscr{X}}$ contains $I^{N}$ entries. Algorithms for tensor decomposition therefore inherit this exponential dependence on $N$. In the case of ALS, each iteration computed via (9) requires access to all entries of $\bm{\mathscr{X}}$. Moreover, for every iteration, it requires forming the matrix $\bm{A}^{\neq m}$ and then solving a linear system involving this matrix.

As discussed in Section 2, several previous works have sought to address the curse of dimensionality in ALS algorithms by sampling the least squares problem in (9). These works develop methods for specific tensor decompositions. By contrast, our sampling framework described in Section 4 can be used to develop a sampling-based decomposition scheme for any TN decomposition. This is done by independently sampling each least squares problem on line 1 in Algorithm 1 according to the leverage scores of the TN matrix $\bm{A}^{\neq m}$. The sampling is done by first determining the indices to sample by applying our techniques in Section 4. Letting $\bm{S}$ represent the sampling operation, the next step is to compute $\bm{S}\bm{X}^{(m)}$ and $\bm{S}\bm{A}^{\neq m}$. The final step is to solve the sampled least squares problem

In order to make the discussion in Sections 4 and 5 a bit more concrete, we now apply these ideas to the CP decomposition.

We run the same feature extraction and classification experiment as in [18]. The experiment considers 7200 images from the COIL-100 dataset [22]. These images depict 100 different objects (e.g., a toy car, an onion, a tomato), each photographed from 72 different angles. Each image is 128 by 128 pixels and has three color channels. We arrange these images into a tensor $\bm{\mathscr{X}}$ of size $7200\times 128\times 128\times 3$ and then decompose it using either a rank-25 CP decomposition or a rank-$(5,5,5,5)$ tensor ring decomposition. For the CP decomposition, we treat the rows of the factor matrix $\bm{A}^{(1)}\in\mathbb{R}^{7200\times 25}$ as feature vectors for each image. We treat 90% of the images as labeled and use them to classify the images in the remaining 10% using a $k$-nearest neighbor algorithm with $k=1$ and 10-fold cross validation. The tensor ring decomposition is used in the same way, but instead of the CP factor matrix we now use the core tensor $\bm{\mathscr{A}}^{(1)}\in\mathbb{R}^{5\times 7200\times 5}$ corresponding to the first mode and reshape the slices $(\bm{\mathscr{A}}^{(1)}_{:i:})_{i}$ into length-25 feature vectors. Further details on the experiment setup are given in Section C.

Table 3 shows the run time, relative decomposition error (Err.) and classification accuracy (Acc.) for a number of different CP and tensor ring decomposition methods. The algorithms for our TNS-CP and TNS-TR were implemented²²2Our code is available at https://github.com/OsmanMalik/TNS. by modifying the codes for CP-ALS-ES and TR-ALS-ES by [18] appropriately.

All methods achieve roughly the same decomposition error, with the randomized methods typically having a slightly higher error. All methods except CPD-NLS yield similar classification accuracy. It is clear that our methods are faster than CP-ALS-ES and TR-ALS-ES from [18], corroborating the improved complexity reported in Tables 1 and 2.

We have presented a sampling-based ALS approach for decomposing a tensor into any TN format. The generality of the framework is notable—all previous works on randomized tensor decomposition we are aware of develop methods for specific decompositions. Additionally, unlike most previous sampling-based ALS algorithms in the literature, our framework makes it possible to do the sampling according to the exact leverage scores. Our approach simplifies and generalizes the scheme developed in [18]. Both complexity analyses and numerical experiments confirm the improved run time we are able to achieve compared to [18] when we use our framework for CP and tensor ring decomposition.

There are many exciting avenues for future research. One direction is to adapt and implement our method for use in a high-performance distributed memory setting. Another interesting direction is designing an exact leverage score sampling scheme for general TN matrices $\bm{A}\in\mathbb{R}^{I\times R}$ that avoid the $\Omega(R^{3})$ cost that we incur when computing the pseudoinverse $\bm{\Phi}=(\bm{A}^{\top}\bm{A})^{+}$. Such methods are known for certain matrices with Kronecker product structure [20], but for general TN matrices such results are not known.

OAM was supported by the LDRD Program of Lawrence Berkeley National Laboratory under U.S. DOE Contract No. DE-AC02-05CH11231. VB was supported in part by the U.S. DOE, Office of Science, Office of Advanced Scientific Computing Research, DOE Computational Science Graduate Fellowship under Award No. DE-SC0022158. RM was supported by NSF Grant No. 2004235.