Fast and Accurate Randomized Algorithms for Low-rank Tensor Decompositions

Abstract:
Low-rank Tucker and CP tensor decompositions are powerful tools in data analytics. The widely used alternating least squares (ALS) method, which solves a sequence of over-determined least squares subproblems, is costly for large and sparse tensors. We propose a fast and accurate sketched ALS algorithm for Tucker decomposition, which solves a sequence of sketched rank-constrained linear least squares subproblems. Theoretical sketch size upper bounds are provided to achieve O( ) relative error for each subproblem with two sketching techniques, TensorSketch and leverage score sampling. Experimental results show that this new ALS algorithm, combined with a new initialization scheme based on the randomized range finder, yields decomposition accuracy comparable to the standard higher-order orthogonal iteration (HOOI) algorithm. The new algorithm achieves up to 22.0% relative decomposition residual improvement compared to the state-of-the-art sketched randomized algorithm for Tucker decomposition of various synthetic and real datasets. This Tucker-ALS algorithm is further used to accelerate CP decomposition, by using randomized Tucker compression followed by CP decomposition of the Tucker core tensor. Experimental results show that this algorithm not only converges faster, but also yields more accurate CP decompositions.

text: Low-rank Tucker and CP tensor decompositions are powerful tools in data analytics. The widely used alternating least squares (ALS) method, which solves a sequence of over-determined least squares subproblems, is costly for large and sparse tensors. We propose a fast and accurate sketched ALS algorithm for Tucker decomposition, which solves a sequence of sketched rank-constrained linear least squares subproblems. Theoretical sketch size upper bounds are provided to achieve O( ) relative error for each subproblem with two sketching techniques, TensorSketch and leverage score sampling. Experimental results show that this new ALS algorithm, combined with a new initialization scheme based on the randomized range finder, yields decomposition accuracy comparable to the standard higher-order orthogonal iteration (HOOI) algorithm. The new algorithm achieves up to 22.0% relative decomposition residual improvement compared to the state-of-the-art sketched randomized algorithm for Tucker decomposition of various synthetic and real datasets. This Tucker-ALS algorithm is further used to accelerate CP decomposition, by using randomized Tucker compression followed by CP decomposition of the Tucker core tensor. Experimental results show that this algorithm not only converges faster, but also yields more accurate CP decompositions.

heading: 1 Introduction
text: Tensor decompositions [31] are general tools for compressing, approximating, as well as extracting important features from high dimensional data, and are widely used in both scientific computing [54, 25, 26] and machine learning [4, 59, 55]. In this paper, we focus on Tucker decomposition [67] and CANDECOMP/PARAFAC (CP) decomposition [24, 23]. The alternating least squares (ALS) method is widely used to compute both decompositions. The ALS algorithm consists of sweeps, and each sweep updates every factor matrix once in a fixed order. The ALS method for Tucker decomposition, called the higher-order orthogonal iteration (HOOI) [5, 16, 31], updates one of the factor matrices along with the core tensor at a time. Similarly, each update procedure in the ALS algorithm for CP decomposition (CP-ALS) updates one of the factor matrices. For both decompositions, solutions to each optimization subproblem guarantee decrease of the decomposition residual.
In this work, we consider decomposition of order N tensors (T ) that are large in dimension size (s) and can be potentially sparse. We focus on the problem of computing low-rank (with target rank R s and R nnz(T )) decompositions for such tensors via ALS, which is often used for extracting principal component information from large-scale datasets. For Tucker decomposition, ALS is bottlenecked by the operation called the tensor times matrix-chain (TTMc). For CP decomposition, ALS is bottlenecked by the operation called the matricized tensor-times Khatri-Rao product (MTTKRP). Both TTMc and MTTKRP have a per-sweep cost of Ω(nnz(T )R) [62]. Consequently,
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
the per-sweep costs of both HOOI and CP-ALS are proportional to the number of nonzeros in the tensor, which are expensive for large tensors with billions of nonzeros.
Recent works have applied different randomized techniques to accelerate both CP [32, 2, 14, 70] and Tucker decompositions [13, 3, 12, 40, 70, 65]. For Tucker decomposition, these randomized algorithms apply sketching techniques to the higher-order singular value decomposition (HOSVD). To do so, they calculate each factor matrix by applying randomized SVD on each matricization of the input tensor. These methods calculate the core tensor via TTMc among the input tensor and all the factor matrices, which incurs a cost of O(nnz(T )R + sN−1R2) for sparse tensors and is still expensive. In addition, HOSVD generates decompositions that are generally less accurate compared to HOOI.
Becker and Malik [39] introduce a sketched ALS algorithm for Tucker decomposition, which avoids the expensive cost of TTMc. Unlike the traditional HOOI, each sweep of this ALS scheme contains N + 1 subproblems, where only one of the factor matrices or the core tensor is updated in each subproblem. This scheme is easier to analyze theoretically, since each subproblem is an unconstrained linear least squares problem, which can be efficiently solved via sketching. However, the scheme produces decompositions that are generally less accurate than HOOI.

heading: 1.1 Our Contributions
text: In this work, we propose a new sketched ALS algorithm for Tucker decomposition. Different from Becker and Malik [39], our ALS scheme is the same as HOOI, where one of the factor matrices along with the core tensor are updated in each subproblem. This guarantees the algorithm can reach the same accuracy as HOOI with sufficiently large sketch size. Experimental results show that it provides more accurate results compared to those in [39].
In this scheme, each subproblem is a sketched rank-constrained linear least squares problem, with the left-hand-side matrix with size sN−1 ×RN−1 composed of orthonormal columns. To the best of our knowledge, the relative error analysis of sketching techniques for this problem have not been discussed in the literature. Existing works either only provide sketch size upper bounds for the relaxed problem [57], where rank constraint is relaxed with a nuclear norm constraint, or provide upper bounds for general constrained problems [69]. We provide tighter sketch size upper bounds to achieve O( ) relative error with two state-of-the-art sketching techniques, TensorSketch [52] and leverage score sampling [18].
With leverage score sampling, our analysis shows that with probability at least 1 − δ, the sketch size of O ( RN−1/( 2δ) ) is sufficient for results with O( )-relative error. With TensorSketch, the
sketch size upper bound is O ( (RN−1 · 3N−1) · (RN−1 + 1/ 2)/δ ) , at least O ( 3N−1 ) times that for leverage score sampling. For both techniques, our bounds are at most O(1/ ) times the sketch size upper bounds for the unconstrained linear least squares problem.
The upper bounds suggest that under the same accuracy criteria, leverage score sampling potentially needs smaller sketch size for each linear least squares problem and thus can be more efficient than TensorSketch. Therefore, with the same sketch size, the accuracy with leverage score sampling can be better. However, with the standard random initializations for factor matrices, leverage score sampling can perform poorly on tensors with high coherence [10] (the orthogonal basis for the row space of each matricization of the input tensor has large row norm variability), making it less robust than TensorSketch. To improve the robustness of leverage score sampling, we introduce an algorithm that uses the randomized range finder (RRF) [22] to initialize the factor matrices. The initialization scheme uses the composition of CountSketch and Gaussian random matrix as the RRF embedding matrix, which only requires one pass over the non-zero elements of the input tensor. Our experimental results show that the leverage score sampling based randomized algorithm combined with this RRF scheme performs well on tensors with high coherence.
For R s, our new sketching based algorithm for Tucker decomposition can also be used to accelerate CP decomposition. Tucker compression is performed first, and then CP decomposition is applied to the core tensor [70, 9, 20]. Since the per-sweep costs for both sketched Tucker-ALS and sketched CP-ALS are comparable, and Tucker-ALS often needs much fewer sweeps than CP-ALS (Tucker-ALS typically converges in less than 5 sweeps based on our experiments), this Tucker + CP method can be more efficient than directly applying randomized CP decomposition [32, 14] on the input tensor.
In summary, this paper makes the following contributions.
• We introduce a new sketched ALS algorithm for Tucker decomposition, which contains a sequence of sketched rank-constrained linear least squares subproblems. Experimental results show that the algorithm yields decomposition accuracy comparable to HOOI, and provides up to 22.0% relative decomposition residual improvement compared to the previous randomized Tucker algorithm.
• We provide theoretical upper bounds for the sketch size of both leverage score sampling and TensorSketch, which ensure that each sketched rank-constrained linear least squares incurs O( ) relative error with high probability.
• We provide detailed comparison of TensorSketch and leverage score sampling in terms of efficiency and accuracy. Our theoretical analysis shows that leverage score sampling is better in terms of both metrics.
• We propose an initialization scheme based on RRF that improves the accuracy of leverage score sampling based sketching algorithm on tensors with high coherence.
• We show that CP decomposition can be more efficiently and accurately calculated based on the sketched Tucker + CP method, compared to directly performing sketched CP-ALS on the input tensor.

heading: 2 Background
text: We introduce the notation used throughout the paper, and briefly review ALS algorithms for Tucker and CP decompositions, and TensorSketch as well as leverage score sampling in this section. We present additional backgrounds, including the pseudo-codes of ALS algorithms for both Tucker and CP decompositions, and the previous work in Appendix A.

heading: 2.1 Notation
text: Our analysis makes use of tensor algebra in both element-wise equations and specialized notation for tensor operations [31]. Vectors are denoted with bold lowercase Roman letters (e.g., v), matrices are denoted with bold uppercase Roman letters (e.g., M), and tensors are denoted with bold calligraphic font (e.g., T ). An order N tensor corresponds to an N -dimensional array. Elements of vectors, matrices, and tensors are denoted in parentheses, e.g., v(i) for a vector v, M(i, j) for a matrix M, and T (i, j, k, l) for an order 4 tensor T . The ith column of M is denoted by M(:, i), and the ith row is denoted by M(i, :). Parenthesized superscripts are used to label different vectors, matrices and tensors (e.g. T (1) and T (2) are unrelated tensors). Number of nonzeros of the tensor T is denoted by nnz(T ). The pseudo-inverse of matrix A is denoted with A†. The Hadamard product of two matrices is denoted with ∗. The outer product of two or more vectors is denoted with ◦. The Kronecker product of two vectors/matrices is denoted with⊗. For matrices A ∈ Rm×k and B ∈ Rn×k, their Khatri-Rao product results in a matrix of size (mn) × k defined by A B = [A(:, 1) ⊗ B(:, 1), . . . ,A(: , k)⊗B(:, k)]. The mode-n tensor times matrix of an order N tensor T ∈ Rs1×···×sN with a matrix A ∈ RJ×sn is denoted by T ×n A, whose output size is s1 × · · · × sn−1 × J × sn+1 × · · · × sN . Matricization is the process of unfolding a tensor into a matrix. The mode-n matricized version of T is denoted by T(n) ∈ Rsn×K where K = ∏N m=1,m6=n sm. We use [N ] to denote {1, . . . , N}. Õ denotes the asymptotic cost with logarithmic factors ignored.
Tucker decomposition with ALS. Throughout the analysis we assume the input tensor has order N and size s× · · · × s, and the Tucker ranks are R× · · · ×R. Tucker decomposition approximates a tensor by a core tensor contracted along each mode with matrices that have orthonormal columns. The goal of Tucker decomposition is to minimize the objective function, f(C ,A(1), . . . ,A(N)) := 1 2
∥∥T −C ×1 A(1) ×2 A(2) · · · ×N A(N)∥∥2F . The core tensor C is of order N with dimensions R× · · · × R. Each matrix A(n) ∈ Rs×R for n ∈ [N ] has orthonormal columns. The ALS method for Tucker decomposition [5, 16, 31], called the higher-order orthogonal iteration (HOOI), proceeds by updating one of the factor matrices along with the core tensor at a time. The nth subproblem can be written as
min C ,A(n)
1
2 ∥∥∥P(n)CT(n)A(n)T −TT(n)∥∥∥2 F , (2.1)
where P(n) = A(1) ⊗ · · · ⊗A(n−1) ⊗A(n+1) ⊗ · · · ⊗A(N). This problem can be formulated as a rank-constrained linear least squares problem,
min B(n)
1
2 ∥∥∥P(n)B(n) −TT(n)∥∥∥2 F , such that rank(B(n)) ≤ R. (2.2)
A(n) corresponds to the right singular vectors of the optimal B(n), while CT(n) = B (n)A(n). Since P(n) contains orthonormal columns, the optimal B(n) can be obtained by calculating the Tensor Times Matrix-chain (TTMc),
Y (n) = T ×1 A(1)T · · · ×n−1 A(n−1)T ×n+1 A(n+1)T · · · ×N A(N)T , (2.3)
and taking B(n) to be the transpose of the mode-n matricized Y (n), Y(n)T(n) . Calculating Y (n) costs O ( sNR ) for dense tensors andO ( nnz(T )RN−1 ) for sparse tensors. Before the HOOI procedure, the factor matrices are often initialized with the higher-order singular value decomposition (HOSVD) [15, 67]. HOSVD computes the truncated SVD of each T(n) ≈ U(n)Σ(n)V(n)T , and sets A(n) = U(n) for n ∈ [N ]. If performing SVD via randomized SVD [22], updating A(n) for each mode costs O ( sNR ) for dense tensors, and costs O ( sN−1R2 + nnz(T )R ) for sparse tensors.
CP decomposition with ALS. CP tensor decomposition [24, 23] decomposes the input tensor into a sum of outer products of vectors. Throughout analysis we assume the input tensor has order N and size s × · · · × s, and the CP rank is R. The goal of CP decomposition is to minimize the objective function, f(A(1), . . . ,A(N)) := 12 ∥∥∥T −∑Rr=1 A(1)(:, r) ◦ · · · ◦A(N)(:, r)∥∥∥2 F , where A(i) ∈ Rs×R for i ∈ [N ] are called factor matrices. CP-ALS is the mostly widely used algorithm to get the factor matrices. In each ALS sweep, we solve N subproblems, and the objective for the update of A(n), with n ∈ [N ], is expressed as,
A(n) = arg min A
1
2 ∥∥∥P(n)AT −XT(n)∥∥∥2 F , (2.4)
where P(n) = A(1) · · · A(n−1) A(n+1) · · · A(N). Solving the linear least squares problem above has a cost of O ( sNR ) . For instance, when solving via normal equations the term
P(n)TXT(n), which is called MTTKRP, needs to be calculated, and it costsO ( sNR ) for dense tensors andO(nnz(T )R) for sparse tensors. A major disadvantage of CP-ALS is its slow convergence. There are many cases where CP-ALS takes a large number of sweeps to converge when high resolution is necessary [42]. When R < s, the procedure can be accelerated by performing Tucker-ALS first, which typically converges in fewer sweeps, and then computing a CP decomposition of the core tensor [11, 70, 9], which only has O ( RN ) elements.
TensorSketch and leverage score sampling. In this paper, we sketch the linear least squares problems using both TensorSketch and leverage score sampling. TensorSketch is a special type of CountSketch, where the hash map is restricted to a specific format to allow fast multiplication of the sketching matrix with the chain of Kronecker products. Leverage score sampling picks important rows based on leverage scores to form the sampled/sketched problem. Both algorithms can be efficiently applied to a chain of Kronecker products, and the detailed analysis is presented in Appendix B.
In the paper, we test two forms of leverage score sampling, random sampling, where we perform importance random sampling based on leverage scores, and deterministic sampling [28], where we deterministically sample rows having the largest leverage scores. Both ideas are also used in [32] for randomized CP decomposition. Papailiopoulos et al. [53] show that if the leverage scores follow a moderately steep power-law decay, then deterministic sampling can be provably as efficient and even better than the random sampling. We compare both leverage score sampling techniques in Section 5.

heading: 3 Sketched Rank-constrained Linear Least Squares
text: Each subproblem of Tucker HOOI solves a linear least squares problem with the following properties, 1) the left-hand-side matrix is a chain of Kronecker products of factor matrices, 2) the left-hand-side matrix has orthonormal columns, since each factor matrix has orthonormal columns, 3) the rank of the output solution is constrained to be less than R, as is shown in (2.2). To the best of our knowledge,
the relative error analysis of sketching techniques for this problem have not been discussed in the literature. In the following two theorems, we will show the sketch sizes of TensorSketch and leverage score sampling that are sufficient for the relative residual norm error of the problems to be bounded by O( ) with at least 1− δ probability. The detailed proofs are presented in Appendix F. Theorem 3.1 (TensorSketch for Rank-constrained Linear Least Squares). Consider matrices P = A(1) ⊗ A(2) ⊗ · · · ⊗ A(N−1), where each A(i) ∈ Rs×R has orthonormal columns, s > R, and B ∈ RsN−1×n. Let S ∈ Rm×sN−1 be an order N − 1 TensorSketch matrix. Let X̃r be the best rank-R approximation of the solution of the problem minX ‖SPX− SB‖F , and let Xr = arg minX,rank(X)=R‖PX−B‖F . With
m = O ( (R(N−1) · 3N−1) · (R(N−1) + 1/ 2)/δ ) , (3.1)
the approximation, ∥∥∥AX̃r −B∥∥∥2
F ≤ (1 +O( )) ∥∥∥AXr−B∥∥∥2 F , holds with probability at least 1− δ.
Theorem 3.2 (Leverage Score Sampling for Rank-constrained Linear Least Squares). Given matrices P = A(1) ⊗ A(2) ⊗ · · · ⊗ A(N−1), where each A(i) ∈ Rs×R has orthonormal columns, s > R, and B ∈ RsN−1×n. Let S ∈ Rm×sN−1 be a leverage score sampling matrix for P. Let X̃r be the best rank-R approximation of the solution of the problem minX ‖SPX− SB‖F , and let Xr = arg minX,rank(X)=R‖PX−B‖F . With m = O ( R(N−1)/( 2δ)
) , the approximation,∥∥∥AX̃r −B∥∥∥2
F ≤ (1 +O( )) ∥∥∥AXr −B∥∥∥2 F , holds with probability at least 1− δ.
Therefore, for leverage score sampling, O ( R(N−1)/( 2δ) ) number of samples are sufficient to get (1 + O( ))-accurate residual with probability at least 1 − δ. The sketch size upper bound for TensorSketch is higher than that for leverage score sampling, suggesting that leverage score sampling is better. As can be seen in (3.1), when RN−1 ≤ 1/ 2, the sketch size bound for TensorSketch is O ( 3N−1 ) times that for leverage score sampling. When RN−1 > 1/ 2, the ratio is even higher. The accuracy comparison of the two methods is discussed further in Section 5.
While TensorSketch has a worse upper bound compared to leverage score sampling, it is more flexible since the sketching matrix is independent of the left-hand-side matrix. One can derive a sketch size bound that is sufficient to get (1 +O( ))-accurate residual norm for linear least squares with general (not necessarily rank-based) constraints based on existing proof techniques (detailed in Appendix G). Although that bound is applicable for general constraints, it is looser than (3.1). For leverage score sampling, we do not provide a sample size bound for general constrained linear least squares.
We also compare the sketch size upper bounds for rank-constrained linear least squares and unconstrained linear least squares in Table 1. For both leverage score sampling and TensorSketch, the upper bounds for rank-constrained problems are at most O(1/ ) times the upper bounds for unconstrained linear least squares problem. The error of sketched rank-constrained solution consists of two parts, the error of the sketched unconstrained linear least squares solution, and the error from low-rank approximation of the unconstrained solution. To make sure the second error term has a relative error bound of O( ), we restrict the first error term to be relatively bounded by O ( 2 ) , incurring a larger sketch size upper bound.

heading: 4 Main Algorithm
text: Our main algorithm is presented in Algorithm 1. To improve the robustness of leverage score sampling, we use an initialization scheme that uses the randomized range finder (RRF) [22] to initialize the factor matrices (lines 3-5). In this scheme, the composition of CountSketch and Gaussian random matrix is used as the RRF embedding matrix, which only requires one pass over the non-zero elements of the input tensor. The detailed initialization algorithm and its cost analysis is detailed in Appendix C.
Algorithm 1 Sketch-Tucker-ALS: Sketched ALS procedure for Tucker decomposition 1: Input: Input tensor T ∈ Rs1×···×sN , Tucker ranks {R1, . . . , RN}, maximum number of sweeps Imax, sketching tolerance
2: C ← O 3: for n ∈ {2, . . . , N} do 4: A(n) ← Init-RRF(T(n), Rn, ) 5: end for 6: for i ∈ {1, . . . , Imax} do 7: for n ∈ {1, . . . , N} do 8: Build the sketching matrix S(n)
9: Y ← S(n)T(n) 10: Z← S(n)(A(1) ⊗ · · · ⊗A(n−1) ⊗A(n+1) ⊗ · · · ⊗A(N)) 11: CT(n),A
(n) ← RSVD-LRLS(Z,Y, R) 12: end for 13: end for 14: return { C ,A(1), . . . ,A(N) }
Algorithm 2 RSVD-LRLS: Low-rank approximation of least squares solution via randomized SVD 1: Input: Left-hand-side matrix Z ∈ Rm×r, right-hand-side matrix Y ∈ Rm×s, rank R 2: Initialize S ∈ Rs×O(R) as a random Gaussian sketching matrix 3: B← (ZTZ)−1 4: C← BZTYS 5: Q,R← qr(C) 6: D← QTBZTY 7: U,Σ,V← svd(D) 8: return QU(:, : R)Σ(: R, : R),V(:, : R)
We provide detailed cost analysis for Algorithm 1. Note that for leverage score sampling, lines 8 and 9 need to be recalculated for every sweep, since S(n) is dependent on the factor matrices. On the other hand, the TensorSketch embedding is oblivious to the state of the factor matrices, so we can choose to use the same S(n) for all the sweeps for each mode n to save cost. This strategy is also
used in [39]. Detailed cost analysis for each part of Algorithm 1 is listed below, where we assume s1 = · · · = sN = s and R1 = · · · = RN = R. For line 3-5, the cost is O ( N nnz(T ) +NsR3/ ) by the analysis in Appendix C. For line 8, if using leverage score sampling, the cost is O(sR) per subproblem (for computing the leverage scores of the previously updated A(i)). If using TensorSketch, the cost is O(Ns), which is only incurred for the first sweep. For line 9, if using leverage score sampling, the cost is O(ms) per subproblem; if using TensorSketch, the cost is O(N nnz(T )), and is only incurred for the first sweep. For line 10, if using leverage score sampling, the cost is O ( mRN−1 ) per subproblem; if using TensorSketch, the cost is
O ( NsR+m log(m)RN−1 ) per subproblem, by the analysis in Appendix B. For line 11, the cost is
O ( msR+mR2(N−1) ) per subproblem, under the condition that m ≥ RN−1 and using randomized SVD as detailed in Algorithm 2.
Therefore, the cost for each subproblem (lines 8-11) is O ( msR+mR2(N−1) ) , for both leverage score sampling and TensorSketch. For TensorSketch, another cost of O(N nnz(T )) is incurred at the first sweep to sketch the right-hand-side matrix, which we refer to as preparation cost. Using the initialization scheme based on RRF to initialize the factor matrices would increase the cost of both sketching techniques by O ( N nnz(T ) +NsR3/ ) .
We compare the cost of each linear least squares subproblem between our sketched ALS algorithms with both HOOI and the sketched ALS algorithms introduced in [39] in Table 2. For the ALS + TensorSketch algorithm in [39],N+1 subproblems are solved in each sweep, and in each subproblem either one factor matrix or the core tensor is updated based on the sketched unconstrained linear least squares solutions. For the ALS + TTMTS algorithm, TensorSketch is simply used to accelerate the TTMc operations, and it has been shown to be less accurate compared to the reference ALS + TensorSketch algorithm in [39].
For the solutions of sketched linear least squares problems to be unique, we need m ≥ RN−1 and hence m = Ω(RN−1). With this condition, the cost of each linear least squares subproblem of our sketched ALS algorithms is less than that for ALS + TTMTS, but is more expensive with related to R compared to the ALS + TensorSketch algorithm in [39], since our cost involves a term mR2(N−1). However, as will be discussed in Appendix E.1, this term does not dominate in the low-rank decomposition regime. In addition, as shown in Section 5, our algorithms provide better accuracy as a result of updating more variables at a time. We also show the sketch size upper bound sufficient to get a (1 +O( ))-accurate approximation in residual norm. As can be seen in the table, our sketching algorithm with leverage score sampling has the smallest sketch size, making it the best algorithm considering both the cost of each subproblem and the sketch size. In [39], the authors give an error bound for the approximate matrix multiplication in ALS + TTMTS, but the relative error of the overall linear least squares problem is not given. For the ALS + TensorSketch algorithm in [39], the sketch size upper bound in Table 2 comes from the upper bound for the unconstrained linear least squares problem.
Note that the analysis generalizes to the case with non-uniform input tensor dimensions and Tucker ranks. For the decomposition of an orderN tensor with dimensions s1×· · ·×sN and the Tucker ranks R1 × · · · ×RN , the least squares subproblem cost for the ith mode for both ALS with TensorSketch and ALS with leverage score sampling generalize from O(msR+mR2(N−1)) (shown in Table 2) to O(msiRi + m ∏N j=1,j 6=iR 2 j ). For ALS with leverage score sampling, the sketch size bound
changes to O( ∏N
j=1,j 6=iRj/( 2δ)). For ALS with TensorSketch, the sketch size bound changes to O ( 3N−1 ∏N j=1,j 6=iRj · ( ∏N j=1,j 6=iRj + 1/ 2)/δ ) .
Algorithm 1 can also be used to accelerate CP decomposition when R s. Tucker compression is performed first, and then CP decomposition is applied to the core tensor. The detailed algorithm and the cost analysis is presented in Appendix D.

heading: 5 Experiments
text: In this section, we compare our randomized algorithms with reference algorithms for both Tucker and CP decompositions on several synthetic and real tensors. We evaluate accuracy based on the final
fitness f for each algorithm, defined as f = 1− ‖T −T̃ ‖F‖T ‖F , where T is the input tensor and T̃ is the reconstructed low-rank tensor. For Tucker decomposition, we focus on the comparison of accuracy and robustness of attained fitness across various synthetic datasets for different algorithms. For CP decomposition, we focus on the comparison of accuracy and sweep count. Our experiments are carried out on an Intel Core i7 2.9 GHz Quad-Core machine using NumPy [50] routines in Python.

heading: 5.1 Experiments for Tucker Decomposition
text: We compare five different algorithms for Tucker decomposition. Two baselines from previous work are considered, standard HOOI and the original TensorSketch-based randomized Tucker-ALS algorithm, which optimizes only one factor in Tucker decomposition at a time [39]. We compare these to our new randomized algorithm (Algorithm 1) based on TensorSketch, random leverage score sampling, and deterministic leverage score sampling. For each randomized algorithm, we test both random initialization for factor matrices as well as the initialization scheme based on RRF. For the baseline HOOI algorithm, we report the performance with both random and HOSVD initializations. We use the following four synthetic tensors and real datasets to evaluate these algorithms.
1. Dense tensors with specific Tucker rank. We create order 3 tensors based on randomlygenerated factor matrices B(n) ∈ Rs×Rtrue and a core tensor C ,
T = C ×1 B(1) ×2 B(2) ×3 B(3). (5.1) Each element in the core tensor and the factor matrices are i.i.d. normally distributed random variables N (0, 1). The ratio Rtrue/R, where R is the decomposition rank, is denoted as α. 2. Dense tensors with strong low-rank signal. We also test on dense tensors with strong low-rank signal, T (b) = T + ∑n i=1 λia (1) i ◦ a (2) i ◦ a (3) i . T is generated based on (5.1), and each vector
a (j) i has unit 2-norm. The magnitudes λi for i ∈ [n] are constructed based on the power-law distribution, λi = C ‖T ‖F i1+η . In our experiments, we set n = 5, C = 3 and η = 0.5. This tensor is used to model data whose leading low-rank components obey the power-law distribution, which is common in real datasets.
3. Tensors with large coherence. We also test on tensors with large coherence, T (b) = T +N . T is generated based on (5.1), andN contains n s elements with random positions and same large magnitude. In our experiments, we set n = 10, and each nonzero element inN has the i.i.d. normal distribution N (‖T ‖F / √ n, 1), which means the expected norm ratio E[‖N ‖F /‖T ‖F ] = 1. This
tensor has large coherence and is used to test the robustness of sketching techniques. 4. Real image datasets. We test on two image datasets, COIL-100 [46] and a Time-Lapse hyper-
spectral radiance images dataset called “Souto wood pile” [44], both have been used previously as a tensor decomposition benchmark [7, 70, 36]. Transferring the data into tensor format results in a tensor of size 128× 128× 7200 for COIL-100, and 1024× 1344× 33 for the Time-Lapse dataset.
For all the experiments, we run 5 ALS sweeps unless otherwise specified, and calculate the fitness based on the output factor matrices as well as the core tensor. We observe that 5 sweeps are sufficient for both HOOI and randomized algorithms to converge. For each randomized algorithm, we set the sketch size to be KR2. The constant factor K reveals the accuracy of each subproblem. For the randomized SVD routine in Algorithm 2, we set the dimension sizes of the random matrix S as s × (R + 5), where the oversampling size is 5. We find that this yields accurate randomized SVD solutions. Let CT(n),A (n) be the output of Line 11, Algorithm 1 via calling accurate SVD, and let ĈT(n), Â (n) be the output via calling randomized SVD. We observe that the error ||CT(n)A
(n) − ĈT(n)Â (n)||F is always smaller than 10−10 for all experiments.
We show the experimental results in Fig. 1. As can be seen in the figure, our new randomized ALS scheme, with either leverage score sampling or TensorSketch, outperforms the reference randomized algorithm for all the synthetic and real tensors. The relative fitness improvement ranges from 4.5% (Fig. 1b,3b) to 22.0% (Fig. 1a,3a) when K = 16 for synthetic tensors. With our new randomized scheme, the relative final fitness difference between HOOI and the randomized algorithms is less than 8.5% when K = 16, indicating the efficacy of our new scheme.
Fig. 1a,1b,1c include a comparison between random initialization and the initialization scheme based on RRF detailed in Algorithm 5. For Tensor 1, both initialization schemes have similar performance.
For the deterministic leverage score sampling on Tensor 2 (Fig. 1b), using RRF-based initialization substantially decreases variability of attained accuracy. For leverage score sampling on Tensor 3 (Fig. 1c), we observe that the random initialization is not effective, resulting in approximately zero final fitness. This is because the random initializations are far from the accurate solutions, and some elements with large amplitudes are not sampled in all the ALS sweeps. With the RRF-based initialization, the output fitness of the algorithms based on leverage score sampling is close to HOOI.
Therefore, our proposed initialization scheme is important for improving the robustness of leverage score sampling.
We present additional experiments on dense synthetic tensors in Appendix E.1, where we show the computational cost comparison of different algorithms, the relation between the sketch size and the least squares subproblem accuracy, as well as the perturbation of factor matrices for each randomized algorithm relative to the baseline HOOI.
Although the analysis in Section 3 shows leverage score sampling has a better sketch size upper bound, the random leverage score sampling scheme performs similar to TensorSketch for the tested dense tensors. In Appendix E.1 and Appendix E.2, we provide additional experimental results on sparse tensors, and results with other sketch sizes. Results show that for multiple sparse tensors and several experiments with smaller sketch sizes, leverage score sampling performs better than TensorSketch.

heading: 5.2 Experiments for CP Decomposition
text: We show the efficacy of accelerating CP decomposition via performing Tucker compression first. We compare four different algorithms, the standard CP-ALS algorithm, the Tucker HOOI + CP-ALS algorithm, sketched CP-ALS, where the sketching matrix is applied to each linear least squares subproblem (2.4), as well as the sketched Tucker-ALS + CP-ALS algorithm, where Tucker compression is performed first, and then CP decomposition is applied to the core tensor. Random leverage score sampling is used for sketching, since it has been shown to be efficient for both Tucker (Section 5.1) and CP (reference [32]) decompositions. We use the synthetic tensor to evaluate these four algorithms, T = ∑Rtrue i=1 a (1) i ◦ a (2) i ◦ a (3) i , where each element in a (j) i is an i.i.d normally distributed random variable N (0, 1) with probability p and is zero otherwise. This guarantees that the expected sparsity of T is lower-bounded by 1 − Rtruep3. The ratio Rtrue/R, where R is the decomposition rank, is denoted as α.
We show the detailed fitness-sweeps relation in Fig. 2. The detailed experimental set-up and additional results with different parameter α are presented in Appendix E.3. We observe that for (sketched) CP-ALS, more than 10 sweeps are necessary for the algorithms to converge. On the contrary, less than 5 Tucker-ALS sweeps are needed for the sketched Tucker + CP scheme, making it more efficient. In summary, we observe CP decomposition can be accurately calculated with fewer passes over the tensor data based on the sketched Tucker + CP method.

heading: 6 Conclusions
text: In this work, we propose a fast and accurate sketching based ALS algorithm for Tucker decomposition, which consists of a sequence of sketched rank-constrained linear least squares subproblems. Theoretical sketch size upper bounds are provided to achieve O( )-relative residual norm error for each subproblem with two sketching techniques, TensorSketch and leverage score sampling. For both techniques, our bounds are at most O(1/ ) times the sketch size upper bounds for the unconstrained linear least squares problem. We also propose an initialization scheme based on randomized range finder to improve the accuracy of leverage score sampling based randomized Tucker decomposition of tensors with high coherence. Experimental results show that this new ALS algorithm is more accurate than the existing sketching based randomized algorithm for Tucker decomposition. This Tucker decomposition algorithm also yields an efficient CP decomposition method, where randomized Tucker compression is performed first, and CP decomposition is applied to the Tucker core tensor afterwards. Experimental results show this algorithm not only converges faster, but also yields more accurate CP decompositions.
We leave high-performance implementation of our sketched ALS algorithm as well as testing its performance on large-scale real sparse datasets for future work. Additionally, although our theoretical analysis shows a much tighter sketch size upper bound for leverage score sampling compared to TensorSketch, their experimental performance under the same sketch size for multiple tensors are similar. Therefore, it will be of interest to investigate potential improvements to sketch size bounds for TensorSketch.

heading: Acknowledgments
text: This work is supported by the US NSF OAC via award No. 1942995.

Review summaries:
The authors propose an Alternating Least Squares (ALS) scheme for fitting the Tucker tensor decomposition, based on sketching techniques. In particular, the proposed strategy is to update one of the factor matrices and the core tensor during each iteration, through the use of sketched rank-constrained linear least squares solvers. The theoretical analysis suggests that leverage score sampling may be more accurate than TensorSketch technique, given the same sketch size budget. However, the authors expose that with standard random initializations for factor matrices, leverage score sampling performs poorly on tensors with high coherence; and introduce the randomized range finder (RRF) algorithm so as to provide a better initialization for the factor matrices, leading to better empirical performance.

This paper studies faster tensor decomposition algorithms using randomization. It studies both Tucker and CP tensor decomposition algorithms, both of which solve a sequence of least squares regression problems. Each of these problems are over-constrained, and thus amenable to randomized dimensionality reduction approaches. These methods were implemented using standard numerical packages (NumPy), and tested on both real and synthetic data sets. They obtain significantly faster convergence, as well as lower costs per iteration.

This is paper presents a new randomized sketching algorithm for computing a low rank Tucker and CP decomposition of tensors. The prospered method is based on using sketching to solve the rank-constrained linear least squares subproblems that appear in the popular alternating least squares (ALS) method. Theoretical analysis is presented that give the sampling complexity needed to solve each rank-constrained subproblems to a certain relative error guarantee. Results are presented for two types of sketching, TensorSketch and leverage score sampling. The rank-constrained linear least squares is solved using the randomized SVD. Many numerical results are presented that illustrate the performance of the proposed method.

The paper proposed a fast and accurate sketching based ALS algorithm for Tucker decomposition. The method is composed of a sequence of sketched rank-constrained linear least squares subproblems. In addition, the paper proved that the sketch sizes of TensorSketch and leverage score sampling that are sufficient for the relative residual norm error of the problems to be bounded by $O(\epsilon)$ with at least $1-\delta$ probability.  The experiments show that proposed algorithm is more accurate than the existing sketching based randomized algorithm for Tucker decomposition. The algorithm can also be extended to CP decomposition effectively.