Fast and accurate randomized algorithms for low-rank tensor decompositions
Authors: Linjian Ma
Edgar Solomonik
Venue: NeurIPS 2021 Poster (2021)
Keywords: tensor decomposition
randomized numerical linear algebra
non-convex optimization
alternating least squares

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(\epsilon)$ 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. 

Fast and Accurate Randomized Algorithms for Low-rank Tensor Decompositions
Linjian Ma
Edgar Solomonik
Department of Computer Science University of Illinois at Urbana Champaign 1ma1601 mm166illino01 lliiiis....... edu
Department of Computer Science
University of Illinois Urbana Champaign
colomon28il11 weteeiiiiiiiseee.. edu
Abstract
Low-rank Tucker and CP tensor decompositions are powerful tools data analytics. The widely used alternating least squares (ALS) method, which solves sequence of over-determined least squares subproblems, is costly for large and sparse tensors. We propose fast and accurate sketched ALS algorithm for Tucker decomposi- tion, which solves a sequence of sketched rank-constrained linear least squares subproblems. Theoretical sketch size upper bounds are provided achieve OEe) 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 the standard higher-ords order orthogonal iteration (HOOI) algorithm. The new algorithm achieves up to 22.0% relative decomposition residual improvement compared the state-of-the- randomized algorithm for Tucker decomposition of various synthetic and real datasets. This Tucker-AL.S algorithm is further used to accelerate CP decompo- sition, by using randomized Tucker compression followed by CP decomposition the Tucker core tensor. Experimental results show that this algorithm not only converges faster, but also yields more accurate CP decompositions.
sketched
Introduction
Tensor decompositions [31] are general tools for compressing, approximating, well as extracting important features from high dimensional data, and are widely used both scientific computing [54, 25, 26] and machine learning [4, 59, 55J. In this paper, focus 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 fixed order. The ALS method for Tucker decomposition, called the higher-order orthogonal iteration (HOOI) [5, 16, 31]. updates 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 of the factor matrices. For both decompositions, solutions to each optimization subproblem guarantee decrease of the decomposition residual.
In this work, 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 < and R < nnz(T) decompositions for such tensors via ALS, which is often used for extracting principal component information from acale datasets. For Tucker decomposition, ALS is bottlenecked by the operation called the tensor times matrix-chain (TTMc). For CP decompo- sition, ALS is bottlenecked by the operation called the matricized tensor-tines Khatri-Ran product (MTTKRP). Both TTMc and MTTKRP have per-sweep cost of Q(nnz(T) (T) [62]. Consequently,
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

the 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 cost of ((nnz + s -²) 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 + subproblems, where only one of the factor matrices or the core tensor is updated in each subproblem. This scheme easier 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.
1.1 Our Contributions
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 x 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 (ce) relative error with two state-of-ther ahe sketching techniques, TensorSketch [52] and leverage score sampling [18].
With leverage score sampling, our analysis shows that with probability at least - 8, the sketch size of O ((R~-1/((26) is sufficient for results with ((e)- error. With TensorSketch, the sketch size upper bound is O( leverage score sampling. For both techniques, our bounds are at most O(1/e) times the sketch size
3N- (RN-1 + 1/62 least (33--1) times that for
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 con-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 < 8, 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-AL.S 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.
2

---

In summary, this paper makes the following contributions.
We introduce new sketched ALS algorithm for Tucker decomposition, which contains sequence of sketched con-ttontre 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 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 conk-contrrained linear least squares incurs (e) 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.
2 Background
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 ALSS algorithms for both Tucker and CP decompositions, and the previous work in Appendix A.
2.1 Notation
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 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 [ A 2) for order 4 tensor T. The ith column of M 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. (1) and (2) are unrelated tensors). Number of nonzeros of the tensor T is denoted by mnz(T). The pseudo-invedd inverse of matrix A is denoted with At. The Hadamard product of two matrices is denoted with +. The outer product of two or more vectors is denoted with o. The Kronecker product of two vectors/matrices denoted with @. For matrices A E Rmxk and B Raxk their Khatri-Rao product results in a matrix of size (mn) x k defined by A @ B = (A 1) & 1).,1),..... .,, A(: k) & ((: k)). The mode-n tensor times matrix of order N tensor T E R XsN with matrix A E R'Xx is denoted by T Xn A, whose output size is 81 X ... X Sn-1 x J x 8n+1 X ... x sn. Matricization is the process of unfolding a tensor into a matrix. The mode-n matricized version of is denoted by T(m) E aR x where K -II IIm=1.mm4 8m- We use [N] to denote [I....,N}. Ö denotes the asymptotic cost with logarithmic factors ignored.
Tucker decomposition with ALS. Throughout the analysis we assume the input tensor has order N X s, and the Tucker ranks are R x x R. Tucker decomposition approximates
and size 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, fC, (1), A (N)) :=
(1) X2 (2) XN A 1 The core tensor C is of order N with dimensions R x
X R. Each matrix A (n) E R for E [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 a time. The noth subproblem can be written as
((a) (P() CT A (n)T T 2 (n)
min (x)
(2.1)
(m)
F
3

---

@ ... AM)) This problem can be formulated as where P(n) = A @ ...... A ()-1) @ A (n+1)
rank-constrained linear least squares problem.
B() - T(, ||: such that rank(B (B()) S R.
min
(2.2)
A (m) corresponds the right singular vectors of the optimal B(a), while C (xx B A (w), Since
((n) contains orthonormal columns, the optimal ((x) can be obtained by calculating the Tensor Times Matrix-chain (TTMc),
y'n) = T 1 A (1) ... x--1 A (n-11 1)T x A (a+I)T XN A (N)T,
(2.3)
and taking (()) to be the transpose of the mode-n matricized y(u), Y Calculating y)() costs
(s R) for dense tensors and O (nnz RR) 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) ((())()) ((n) and sets A (m) U
for n E [M] If performing SVD via randomized SVD [22], updating (n) for each mode costs (s R)) for dense tensors, and costs (sN R² + (T) for sparse tensors.
CP decomposition with ALS. CP tensor decomposition [24, 23] decomposes the input tensor into sum of outer products of vectors. Throughout analysis we assume the input tensor has order N and size X ... X 8, and the CP rank is R. The goal of CP decomposition is to minimize the objective function, ((A()...,AA())
T - (()))) o.......... A (t.r) where A (6) R R for E [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 E [NV] is expressed as,
(m) = arg min A 2
p(m) AT
XT (n) F²
(2.4)
where P(n) = (1) ..... .... @ (n-1) O ( .... (N), Solving the linear least squares problem above has cost of Os R). For instance, when solving via normal equations the term p(n) (n) which is called MTTKRP, needs be calculated, and costs (ss R) for dense tensors and (unz R)) for sparse tensors, A major disadvantage of CP-ALS is its slow convergence. There are many cases where CP-ALS takes 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(R" (R*) elements.
TensorSketch and leverage score sampling. In this paper, 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 specifie 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 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 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.
3 Sketched Rank-constrained Linear Least Squares
Each subproblem of Tucker HOOI solves a linear least squares problem with the following properties, 1) the left-hand-side matrix is 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 (c) with least
probability. The detailed proofs are presented in Appendix
Theorem 3.1 (TensorSketch for Rank-constrained Linear Least Squares). Consider matrices P A (1)) @ A (2) @ ... @ A (1, where each A!i) E RaxR has orthonormal columns, >> R, and B E R="=1xx xn_ Let S E Rmxx be an order N - 1 TensorSketch matrix. Let X, be the best rank-R R approximation of the solution of the problem minx ||SPX SBI] and let X, arg milx (ann)))nn PX- B| With
(1 33W-1) (1 (II +1/2)//))
(3.1)
AX, B (1+ -1 AX, BlB holds with probability least 8.
the approximation,
Theorem 3.2 (Leverage Score Sampling for Rank-constrained Linear Least Squares). Given matrices P = A @ A @ ... @ A (N-1), where each A E I has orthonormal columns, s R. and B E R R="-1xn xn, Ler S E R==xs-1 be leverage score sampling matrix for P. Let X, be the best rank- R approximation of the solution of the problem minx ||SPX-SBII SBI] and let X, arg (onk()) PX B p. With m ((R(N-1)((25))) the approximation,
AX, B S (1 (1+
AX, B
holds with probability at least 1 1 8.
Therefore, for leverage score sampling, ((R(N-1)/((P))) number of samples are sufficient to get (1 + O(c))-accurate residual with probability at least - 8. The sketch size upper bound for TensorSketch higher than that for leverage score sampling, suggesting that leverage score sampling is better. As can be seen in (3.1), when R. < $1/2², the sketch size bound for TensorSketch is (33 times that for leverage score sampling. When N > 1/e², the ratio is even higher. The accuracy comparison of the two methods is discussed further in Section
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 derive sketch size (((e)))accuratt -accurate residual norm for linear least squares with general
bound that is sufficient to get (1 + (c)) (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 provide sample size bound for general constrained linear least squares.
Sketching method
Rank-constrained least squares
Unconstrained least squares
Leverage score sampling
o R(N (((N)))((28)) (Theorem 3.2)
o R (R(-1)/(E8)) (Theorem E.11) or o (1(-11) log log(18/2))))) [32]
TensorSketch
o ((3R)N-D)) (RR...((NN))) (R(N-1) 1/62)/6)
o (3R)(N-1) (R
(Theorem 3.1)
(Theorem F.7)
Table 1: Comparison of sketch size upper bounds for rank- constraidd linear least squares and unconstrained linear least squares. The upper bounds are sufficient for the relative residual norm error to be bounded by (e)) with least 1 probability.
We also compare the sketch size upper bounds for sunk-constrained linear least squares and uncon- strained linear least squares in Table 1. For both leverage score sampling and TensorSketch, the upper bounds for monscoontrai problems are at most O(1/c) times the upper bounds for unconstrained linear least squares problem. The error of sketched conk-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 relative emor bound of Oe)), we restrict the first error term to be relatively bounded by Oc2), incurring larger sketch size upper bound.
5

---

Main Algorithm
Our main algorithm is presented Algorithm To improve the robustness of leverage score sampling. an initialization scheme that uses the randomized range finder (RRF) [22] to initialize the factor matrices (lines 3-5). In this scheme, the composition CountSketch and Gaussian random matrix is used 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 Sketch- Sketch-TTucker-ALSS Sketched ALS procedure for Tucker decomposition
Input: Input tensor E x Tucker ranks {Ri,..., Ry }, maximum number of sweeps /mxx. sketching tolerance
2:
3: for [2, do
4: A + Init RRF (n), Rm.
5: end for
6: for (I, mmxx do
7: for n E fI,..., NY. do
8:
Build the sketching matrix (())
9:
Y < S (x)
10:
- ((n) [ @ ... & A (x-1) @ (1+1) 8....8 A())
11:
CT (n) A (n) < RSVD-LALS(Z, Y, R)
12: end for
end for return (C,A , (8))
14:
Algorithm for Tucker
subproblem cost
Sketch size (m)
Prep cost
ALS
Q(nnz(T)R
ALS+TensorSketch [39]
(msR mR
((3R))~1100 /8. (R" 1/c))
ONN nnz(T))
ALS+TTMTS [39]
O (msR
Not shown
((Nnnz())
ALS TensorSketch
0 mss + ml)
((3R)N-118 (NN-1 1/c2))
ON mnz(T))
(Theorem 3.1)
ALS+leverage scores
O (msR + m m
O(RN /((")) (Theorem 3.2)
Table Comparison of algorithm complexity between Tucker-ALS (HOOI), ALS with the TensorS- ketch/leverage score sampling, and the sketched Tucker-ALSS algorithms introduced in [39]. The third column shows the sketch size sufficient for the sketched linear least squares to be (1 + O()) accurate with probability at least 1 - 8. Underlined algorithms are our new contributions.
Algorithm/ 2 RSVD-LRIS: Low-rank approximation of least squares solution via randomized SVD
Input: Left-hand-sidh matrix Z E R right-hand hand-side matrix E mgxx rank
2: Initialize S (=xO(R) as a random Gaussian sketching matrix
3: B t- Z)²
4: C BZ"YY
5: Q,R + qr(C)
Q
6: D
7: U.E.V +- -vs
8: return QU(:: R)E(: R, R), V(t,: 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 ((") dependent on the factor matrices. On the other hand, the TensorSketch embedding is oblivious to the state of the factor matrices,
choose to use the same (()) 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 is listed below, where we assume 81 N s and R
R R.
For line 3-5, the cost is ON (T) + N by the analysis in Appendix C. For line 8, if using leverage score sampling, the cost is (sttt) per subproblem (for computing the leverage scores of the previously updated (.), If using TensorSketch, the cost is O(Ns), (N)., which is only incurred for the first sweep. For line 9, if using leverage score sampling, the cost is (ms) per subproblem; if using TensorSketch, the cost is (N (T)) and is only incurred for the first sweep. For line 10, if using leverage score sampling, the cost is OmmRN-1) per subproblem; if using TensorSketch, the cost (NNs + m log (m)) per subproblem, by the analysis in Appendix B. For line 11, the cost is O (mssR + 12)) per subproblem, under the condition that m 2 R and using randomized SVD as detailed Algorithm 2.
Therefore, the cost for each subproblem (lines 8-11) is O + m m22(NN1) for both leverage score sampling and TensorSketch. For TensorSketch, another cost of O(N mnz(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 (NN (T) + NNRR /e).
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 For the ALS + Skntch algorithm in [39], NN + 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 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 RM-1 and hence m = 12(R2--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 m 1)) However, as will be discussed in Appendix 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 get a (1 + ((e))-accuraee 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 eтror bound for the approximate matrix multiplication in ALS + TTMTS, but the relative error of the overall linear least squares problem not given. For the ALS + TensorSketch algorithm [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 order N tensor with dimensions s1 x x n and the Tucker ranks Rí X ... x Ry, the least squares subproblem cost for the ith mode for both ALS with TensorSketch and ALS with leverage score sampling generalize from O(msR + m m (shown in Table 2) to (ms, R4 + m II Live R2). For ALS with leverage score sampling, the sketch size bound changes to O O ( j=1.74 (=)() For ALS with TensorSketch, the sketch size bound changes to 0 3 11, 1ل ii R2 (II N
ijy R2 1/6))/8).
Algorithm 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.
Experiments
In this section, we compare our randomized algorithms with reference algorithms for both Tucker and CP decompositions several synthetic and real tensors. We evaluate accuracy based on the final
7

---

fitness f for each algorithm, defined as f= T- where is the input tensor and T the 1
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 Intel Core 17 2.9 GHz Quad-Core machine using NumPy [50] routines in Python.
5.1 Experiments for Tucker Decomposition
We compare five different algorithms for Tucker decomposition. Two baselines from previous work are considered, standard HOOI and the original TensorSketch-basee 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 well 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.
Dense tensors with specific Tucker rank. We create order 3 tensors based on randomly- generated factor matrices (z)) E gexBow and core tensor C,
T=C x1 (()) X2 ( X3 B(3)
(5.1)
Each element in the core tensor and the factor matrices are i.i.d. normally distributed random variables . (00 1). The ratio R R, where R is the decomposition rank, is denoted as a. Dense tensors with strong low-rank signal. We also test on dense tensors with strong low-rank signal, (b) T + Σi= ).a oa; a T is generated based on (5.1), and each vector (1) has unit 2-norm. The magnitudes A for i E [nx] are constructed based on the power-law distribution, X; = C I|TII._ In our experiments, we set n = 5, C 3 and 1) = 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.
a,
3. Tensors with large coherence. We also test on tensors with large coherence, T = T
generated based on 1), and N contains n < s elements with random positions and same large magnitude. In our experiments,
set 10, and each nonzero element W has the i.i.d. normal distribution N ([ |r /V, 1). which means the expected norm ratio [[[[[]|||| /|||||||]]] = 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 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 tensor of size 128 X 128 x 7200 for COIL-100, and 1024 x 1344 x 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 KR The constant factor 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 x (R + 5), where the oversampling size is 5. We find that this yields accurate randomized SVD solutions. Let C (n))
A (n) be the output of Line 11, Algorithm via calling accurate SVD, and
let ÖT [n))
A (n) be the output via calling randomized SVD. We observe that the error
|C A (n)
(n)
Cg Á [1)]] 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 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, 1c include comparison between random initialization and the initialization scheme based on RRF detailed in Algorithm 5. For Tensor 1, both initialization schemes have similar performance.
8

---

0.7
0.71
0.70
0.5-
8
0.25
0.68
0.3-
0.87-
0.2
HOEVO/RF
random
0.88
HDSVD/RRF
random
0.1
random
HOSVDIRRF
HOOI
Lev
Leu-xx Method
15
T5-w
H0O
Lav
Leviis Method
18
TS-mn
HOOI
Law
Lew-fix Method
TS
T3-mef
Tensor with
200
(b) Tensor with 200
Tensor 3 with 1000
0.65
0.80
Fitness
Fitness
0.60
HOOI
HOOI
Lev
0.55
TS
TS-ref
TS-ref
0.50
0.70
0
6
8
10
2
4
6
10
Sweeps
Sweeps
(d) COILI00 dataset
(e) Time-Lapse dataset
Figure I: Experimental results for Tucker decomposition. For all experiments, the Tucker rank is x 5 5 and the sketch size parameter = 16. synthetic tensors, we set a = 1.6. HOSVD/RRF means HOOI is initialized with HOSVD, and all other methods are initialized with RRF Lev, Lev-fix TS denote our new sketched Tucker-ALS scheme with leverage score random sampling. leverage score deterministic sampling, and TensorSketch, respectively. TS-ref denotes the reference ALS- TensorSketch algorithm [39]. (a)( Box plots of the final fitness for each algorithm with different synthetic tensors. Each box is based the 25th-75th quartiles, the median is indicated by a horizontal line inside the box, and outliers are displayed dots. (d)(e) Detailed fitness- sweess relation for real image datasets. HOOI is initialized with HOSVD, and all other methods are initialized with RRF.
10 experiments with different random seeds. Each box shows
0.6
0.6
0.4-
Fitness
0.4-
Fitness
Cucker+CP
0.2
Tucker+CP
Tucker+CP
Lev CP
Lev CP
CP
Lev CPcker+
Lev Tucker+CP
Tucker+CP
0.0
0.0
10
10
10
Sweeps
Sweeps
Sweeps
(a) p 0.5
(b) 0.1
(c) p = 0.02
Figure 2: Detailed fitness-sweeps relation for CP decomposition of three tensors with different
parameters. For all the experiments, we set 2000, R = 10, 1.2 and K 16. Markers represent the results per sweep. For Tucker CP algorithms, the fitness shown for ith sweep is the output fitness after running Tucker sweeps along with 25 CP-ALS sweeps core tensors afterwards.
For the deterministic leverage score sampling on Tensor (Fig. 1b), using RRF-based initialization substantially decreases variability of attained accuracy. For leverage score sampling on Tensor (Fig. Ic), 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 all the ALS sweeps. With the RRF-based initialization, the output fitness of the algorithms based on leverage score sampling is close to HOOL

---

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. 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 the baseline HOOI.
Although the analysis Section shows leverage score sampling has 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.
5.2 Experiments for CP Decomposition
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 sub- problem (2.4), as well the sketched TLS 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, Ruse
T D.uu a, a, a,
(3) where each element in a, is an i.i.d normally distributed random variable (0,1) 1) with probability p and is zero otherwise. This guarantees that the expected sparsity of T is lower-bounded by Rmupp" The ratio Rave/ / R. where R is the decomposition rank, is denoted as a.
(1) (2)
We show the detailed fitness-sweeps relation in Fig. 2. The detailed experimental set-up and additional results with different parameter a are presented 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 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.
Conclusions
In this work, we propose fast and accurate sketching based ALS algorithm for Tucker decompo- sition, which consists of sequence of sketched rank-constrained linear least squares subproblems. Theoretical sketch size upper bounds provided to achieve Oellrrelatvve residual norm error for each subproblem with two sketching techniques, TensorSketch and leverage score sampling. For both techniques, our bounds are most (11/e) 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 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.
10

---

Acknowledgments
This work is supported by the US NSF OAC via award No. 1942995.
References
[1] Evrim Acar, Daniel M Dunlavy, and Tamara G Kolda. A scalable optimization approach for fitting canonical tensor decompositions. Journal of Chemometrics, 25(2):67-86 2011.
[2] Kareem S Aggour, Alex Gittens, and Bülent Yener. Adaptive sketching for fast and convergent canonical polyadic decomposition. In International Conference on Machine Learning. PMLR, 2020.
[3] Salman Ahmadi-Asl, Stanislav Abukhovich, Maame Asante-Mensah, Andrzej Cichocki, Anh Huy Phan, Tohishisa Tanaka, and Ivan Oseledets. Randomized algorithms for computation of Tucker decomposition and higher order SVD (HOSVD). IEEE Access, 9:28684-28706, 2021.
[4] Animashree Anandkumar, Rong Ge, Daniel Hsu, Sham M Kakade, and Matus Telgarsky. Tensor decompositions for learning latent variable models. Journal of Machine Learning Research, 15:2773-2832 -8322, 2014.
[5] Claus A Andersson and Rasmus Bro. Improving the speed of multi-way algorithms: Part I. Tucker3. Chemometrics and intelligent laboratory systems, 42(1-2))93-103, 1998.
[6] Haim Avron, Kenneth I. Clarkson, and David P Woodruff. Sharper bounds for regularized data fitting. arXiv preprint 161 2016.
[7] Casey Battaglino, Grey Ballard, and Tamara Kolda. A practical randomized CP tensor decomposition. SIAM Journal on Matrix Analysis and Applications, 39(2):876-901, 2018.
[8] Christos Boutsidis and D. Woodruff. Communication-optinall distrrbbttee analysis in the aolumn-prttitio model. arXiv preprint arXiv: 2015.
distributed principal component
[9] Rasmus Bro and Claus A Andersson. Improving the speed of multiway algorithms: Part II: Compression. Chemometrics and intelligent laboratory systems, 42(1-2):105-113, 1998.
[10] Emmanuel J Candès and Benjamin Recht. Exact matrix completion via convex optimization. Foundations of Computational mathematics, 9(6):717, 2009.
[11] J Douglas Carroll, Sandra Pruzansky, and Joseph B Kruskal. CANDELINC: A general ap- proach to multidimensional analysis of many-way arrays with linear constraints on parameters. Psychometrika, 45(1):3-24 1980.
[12] Maolin Che and Yimin Wei. Randomized algorithms for the approximations of Tucker and the tensor train decompositions. Advances in Computational Mathematics, 45(1):395-428, 2019.
[13] Maolin Che, Yimin Wei, and Hong Yan. Randomized algorithms for the low multilinear rank approximations of tensors. Journal of Computational and Applied Mathematics, page 113380, 2021.
[14] Dehua Cheng, Richard Peng, Yan Liu, and loakeim Perros. SPALS: Fast alternating least squares via implicit leverage scores sampling. Advances in neural information processing systems, 29:721-729, 2016.
[15] Lieven De Lathauwer, Bart De Moor, and Joos Vandewalle. A multilinear singular value decomposition. SIAM journal on Matrix Analysis and Applications. 21(4):1253-122 2000.
[16] Lieven De Lathauwer, Bart De Moor, and Joos Vandewalle. On the best rank-1 and rank- (r',, "2,,r,.... rw)) approximation of higher-order tensors. SIAM journal on Matrix Analysis and Applications, 21(4):1324-132 2000.

---

[17] Huaian Diao, Zhao Song, Wen Sun, and David Woodruff. Sketching for Kronecker product regression and p-splines. In International Conference on Artificial Intelligence and Statisties, pages 1299-1308. PMLR, 2018.
[18] Petros Drineas, Malik Magdon- Ismail, Michael W Mahoney, and David P Woodruff. Fast approximation matrix coherence and statistical leverage. The Journal of Machine Learning Research, 13(1 (31)5:375533 2012.
[19] Petros Drineas, Michael W Mahoney, Shan Muthukrishnan, and Tamás Sarlós. Faster least squares approximation. Numerische mathematik, 117(2):219 249, 2011.
[20] Benjamin Erichson, Krithika Manohar, Steven I. Brunton, and Nathan Kutz. Randomized
CP tensor decomposition. Machine Learning: Science and Technology, (2):022512,, 2020.
[21] Ming Gu and Stanley C Eisenstat. Efficient algorithms for computing strong rank rank-revealing QR factorization. SIAM Journal on Scientific Computing, 17(4):848-86) 869, 1996.
[22] Nathan Halko, Per- Gunnar Martinsson, and Joel Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM review, 53(2):217-288, 2011.
[23] Richard Harshman. Foundations of the PARAFAC procedure: models and conditions for an explanatory multimodal factor analysis. 1970.
[24] Frank L Hitchcock. The expression of a tensor or polyadic as a sum of products. Studies in Applied Mathematics, 6( 6(1-4):164-11, 1927.
[25] Edward G Hohenstein, Robert M Parrish, and Todd Martínez. Tensor hypercontraction density fitting. I. Quartic scaling second-and third-order Moller-Plesse perturbation theory. The Journal of chemical physics, 137(4))044103, 2012.
[26] Felix Hummel, Theodoros Tsatsoulis, and Andreas Grüneis. Low rank factorization of the Coulomb integrals for periodic coupled cluster theory. The Journal of chemical physics, 146(12):124105, 2017.
[27] Ruhui Jin, Tamara Kolda, and Rachel Ward, Faster Johnson-LLindeestrauss transforms via Kronecker products. arXiv preprint arXiv: 09090, 2019.
(28] lan T Jolliffe. Discarding variables in principal component analysis. I: Artificial data. Journal of the Royal Statistical Society: Series C (Applied Statistics), 21(2):160-173, 1972.
[29] Oguz Kaya and Bora Uçar. High performance parallel algorithms for the Tucker decomposition of sparse tensors. In 2016 45th International Conference on Parallel Processing (ICPP), pages 103-112. IEBE, 2016.
[30] Oguz Kaya and Bora Uçar. Parallel CANDECOMPPARAFFAC decomposition of sparse tensors using dimension trees. SIAM Journal on Scientific Computing, 40(1):C99-C130 2018.
[31] Tamara G Kolda and Brett W Bader. Tensor decompositions and applications. SIAM review, 51(3):555-000 2009.
[32] Brett W Larsen and Tamara G Kolda, Practical leverage-based sampling for low-rank tensor decomposition. urXiv preprini arXiv:2006. 16438, 2020.
[33] Hao Li, Zixuan Li, Kenli Li, Jan S Rellermeyer, Lydia Chen, and Keqin Li. SGD_Tuckes Tucker: A novel stochastic optimization strategy for parallel sparse Tucker decomposition. arXiv preprint urXiv:2012.03550, 2020.
[34] Jiajia Li, Jee Choi, Toakeim Perros, Jimeng Sun, and Richard Vuduc. Model-driven sparse CP decomposition for higher-order tensors. In 20/7 IEEE international parallel and distributed processing symposium (IPDPS), pages 1048-1057. IEEE, 2017.
[35] Na Li, Stefan Kindermann, and Carmeliza Navasca. Some convergence results on the regularized alternating squares method for tensor decomposition. Linear Algebra and its Applications, (38(2):7968812. 2013. 2013.
12

---

[36] Linjian Ma and Edgar Solomonik. Accelerating alternating least squares for tensor decomposi- tion by pairwise perturbation. arXiv preprint arXiv:1811.10577 2018.
[37] Linjian Ma and Edgar Solomonik. Efficient parallel CP decomposition with pairwise pertur- bation and multi-sweep dimension tree. In 2021 IEEE International Parallel and Distributed Processing Symposium (IPDPS). pages 412-421. IEEE, 2021.
[38] Michael W Mahoney. Randomized algorithms for matrices and data. arXiv preprint arXiv:1104.5557 2011.
[39] Osman Asif Malik and Stephen Becker. Low-rank Tucker decomposition of large tensors using Tensorsketch. Advances in neural information processing systems, 11:10096100166 2018.
[40] Rachel Minster, Arvind Saibaba, and Misha Kilmer. Randomized algorithms for low-rank tensor decompositions in the Tucker format. SIAM Journal on Mathematics of Data Science, 2(1):189-215, 2020.
[41] Leon Mirsky. Symmetric gauge functions and unitarily invariant norms. The quarterly journal of mathematics, 11(1):50-59, 1960.
[42] Ben C Mitchell and Donald S Burdick. Slowly converging PARAFAC sequences: swamps and two-factor degeneracies. Journal of Chemometrics, 8(2):155-168, 1994.
[43] Drew Mitchell, Nan Yc, and Hans De Sterck. Nesterov acceleration of alternating least squares for canonical tensor decomposition. arXiv preprint arXiv:1810.05847 2018.
[44] Sérgio MC Nascimento, Kinjiro Amano, and David H Foster. Spatial distributions of local illumination color in natural scenes. Vision Research, 120:39-44, 2016.
[45] Carmeliza Navasca, Lieven De Lathauwer, and Stefan Kindermann. Swamp reducing technique for tensor decomposition. In 2008 16th European Signal Processing Conference, pages 1-5. IEEE, 2008.
[46] Sameer A Nene, Shree K Nayar, and Hiroshi Murase. Columbia object image library (COIL- 100).
[47] Dimitri Nion and Lieven De Lathauwer. An enhanced line search scheme for complex-valuee tensor decompositions. Application in DS-CDMA. Signal Processing, 88(3):749-755, 2008.
[48] Israt Nisa, Jiajia Li, Aravind Sukumaran-Rapp Rajam, Richard Vuduc, and P Sadayappan. Load- balanced sparse MTTKRP on GPUs In 2019 IEEE International Parallel and Distributed Processing Symposium (IPDPS), pages 123- 133. IEEE, 2019.
[49] Sejoon Oh, Namyong Park, Sael Lee, and Uksong Kang. Scalable Tucker factorization for sparse tensors- algorrihss and discoveries. In 2018 IEEE 34th International Conference on Data Engineering (ICDE), pages 1 1120-1131. IEEE, 2018.
[50] Travis E Oliphant. A guide NumPy, volume 1. Trelgol Publishing USA, 2006.
[51] Pentti Paatero. A weighted non-negative least squares algorithm for three-way PARAFAC factor analysis. Chemometrics and Intelligent Laboratory Systems, 38(2):223-242, 1997.
[52] Rasmus Pagh. Compressed matrix multiplication. ACM Transactions on Computation Theory (TOCT), 5(3):1177, 2013.
[53] Dimitris Papailiopoulos. Anastasios Kyrillidis, and Christos Boutsidis. Provable deterministic leverage score sampling. In Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 997-1006, 2014.
[54] Will Pazner and Per-Olof Persson. Approximate tensor-produce preconditioners for very high order discontinuous Galerkin methods. Journal of Computational Physics, 354:344-333 --944446,,,,- 2018.
[55] Anh Huy Phan and Andrzej Cichocki. Tensor decompositions for feature extraction and classification of high dimensional datasets. Nonlinear theory and its applications, IFICE, 1(1):37-68, 2010.
13

---

[56] Anh-Huy Phan, Petr Tichavsky, and Andrzej Cichocki. Low complexity damped Gauss-Newtor Newton algorithuss for CANDECOMPPARAFAC. SIAM Journal on Matrix Analysis and Applications. 34(1):126-147, 2013.
[57] Mert Pilanci and Martin J Wainwright. Iterative Hessian sketch: fast and accurate solution approximation for constrained least-squares. The Journal of Machine Learning Research, 17(1):1842-1879 2016.
[58] Myriam Rajih, Pierre Comon, and Richard A Harshman. Enhanced line search: A novel method to accelerate PARAFAC. SIAM journal on matrix analysis and applications, 30(3): 128-1147, 2008.
[59] Nicholas D Sidiropoulos, Lieven De Lathauwer, Xiao Fu, Kejun Huang, Evangelos E Papalex- akis, and Christos Faloutsos. Tensor decomposition for signal processing and machine learning. IEEE Transactions on Signal Processing, 65(13):3551-33)
[60] Navjot Singh, Linjian Ma, Hongru Yang. and Edgar Solomonik. Comparison of accuracy and scalability of geussnn and alternating least squares for CANDECOMC/PARATAC decomposition. SIAM Journal on Scientific Computing, 43(4))C200-C311, 2021.
[61] Shaden Smith and George Karypis. A medium medium-grained algorithm for sparse tensor factorization. In 20/6 IEEE International Parallel and Distributed Processing Symposium (IPDPS), pages 902-911, IEEE, 2016.
[62] Shaden Smith, Niranjay Ravindran, Nicholas D Sidiropoulos, and George Karypis. SPLATT: Efficient and parallel sparse tensor-matrix multiplication. In 2015 IEEE International Parallel and Distributed Processing Symposium, pages 61-70. IEEE, 2015.
[63] Zhao Song, David P Woodruff, and Peilin Zhong. Relative error tensor low rank approximation. In Proceedings of the Thirtieth Annual ACM-SIAM Symposium on Discrete Algorithms, pages 2772-2789. SIAM, 2019.
[64] Yiming Sun, Yang Guo, Charlene Luo, Joel Tropp, and Madeleine Udell. Low-rank Tucker approximation of a tensor from streaming data. SIAM Journal on Mathematics of Data Science, ((4):1123111) 2020.
[65] Yiming Sun, Yang Guo, Joel A Tropp, and Madeleine Udell. Tensor random projection for low memory dimension reduction. In NeurlPS Workshop on Relational Representation Learning, 2018.
[66] Petr Tichavský, Anh Huy Phan, and Andrzej Cichocki. A further improvement of fast damped Gauss-Newton algorithm for SANDECOMP-PARAEES tensor decomposition, In 2013 IEEE International Conference on Acoustics, Speech and Signal Processing, pages 5964-5969999 5968. IEEE, 2013.
[67] Ledyard R Tucker. Some mathematical notes on three-mode factor analysis. Psychometrika,
31 31(3)177--11, 1966.
[68] Shusen Wang. A practical guide to randomized matrix computations with MATLAB implemen- tations. arXiv preprint arXiv:/505.07570 2015.
[69] David P Woodruff. Sketching tool for numerical linear algebra. arXiv preprint arXiv:1411.4335 2014.
[70] Guoxu Zhou, Andrzej Cichocki, and Shengli Xie. Decomposition of big tensors with low multilinear rank. arXiv preprint arXiv:14122 1885, 2014.
14


Reviews:
[{"review_id": "zl0pRx2x0Zi", "reviewer": "Reviewer_UMK5", "summary": "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.", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Tensor factorization and approximation is a fundamental subroutine in learning algorithms. The use of randomization to speed up numerical routines is also a topic that has received much attention. However, I'm not aware of too many hollistic examinations of how to brining randomized tools for matrices to tensor settings, and this work shows significantly more can be done than simply accelerating the linear subroutines on matricizations. I believe some of the ideas introduced may be broadly applicable to other numerical problems on tensors, or even other non-linear numerical problems. So this result should be of wide interest to those working on efficient numerical primitives.\n\nI thank the authors for addressing my concerns, and discussing avenues for incorporating them in future works. However, I feel obligated to take these concerns into account in evaluating the current submitted version.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "A general issue with empirically comparing tensor factorization routines is that the ground truth is unclear, unless the data is synthetic. The experiments seem to show that the proposed new algorithms converge to the same ground truth as non-randomized CP/Tucker. So I'm a bit worried that the 22% increase in relative decomposition residue might be misinterpretted as final residue, instead of at the intermediate iterations.\n\nAlso, usually a key motivation for using randomized methods is to improve time/space efficiency. This seems to be underemphasized (at least compared to residues).", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "vQ5wnWF6YmR", "reviewer": "Reviewer_WEvf", "summary": "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\nby $O(\\epsilon)$ with at least $1-\\delta$ probability.\n\nThe 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.", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper is well-structured and the algorithm is well supported by the theorems. The paper compared the complexity thoroughly, though there is no improvement. The numerical results verified that the proposed algorithm is more accurate than the baselines. \n\nThe paper has the following issues.\n\n1. In Table 2, the sketch size in ALS+TensorSketch [37] could be much smaller than proposed ALS+TensorSketch if $\\epsilon$ is very small (very accurate fitting). When $\\epsilon$ is small enough,  the  sketch size in ALS+TensorSketch [37]  is also much smaller than that of the proposed one.\n\n2. Observed from Table 2, the LS subproblem cost of ALS+TensorSketch [37] is less than proposed two methods,  especially when the rank of ground-truth tensor is high.\n\n3. The abbreviations of the proposed algorithms and the baselines are very confusing. What does TS-ref stand for?\n\n4. In the experiments, where are the results of the two algorithms (ALS+TensorSketch and ALS+TTM) proposed in [37]?\n\n5. The paper did not report the per-sweep cost numerically in the experiment.\n\n6. In Figure (a), the fitness of Lev Tucker+CP dropped significantly at sweep 4. Is the algorithm unstable?\n\n7. It is not clear why the fitness on the synthetic data are always very low, e.g. <0.4 in Figure 1(a).", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "urGF0IxJ2ru", "reviewer": "Reviewer_mFWR", "summary": "\nThis work proposes sketch based tensor decomposition algorithms for accelerating alternating least squares (ALS) based Tucker decomposition and canonical polyadic decomposition. The claim is that new sketch sample complexity bounds are derived, with improved complexity, which can help improve Tucker and CPD reduce fitting errors under the same number of samples used as some baseline sketching algorithms.\n\n", "questions": "", "limitations": "", "rating": 5, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "\n1.\tThe proposed method is feels a bit incremental. The algorithmic structure is a standard one in the tensor community. Applying sketching based techniques to reduce the dimension of the unfolded least squares problems was used in [Battaglino and Kolda 2018] for CPD, and was exploited by a number of follow-up works. Improving sketching complexity upon existing framework is indeed a reasonable direction, but the techniques used in this work seem to be standard.\n\n2.\tThe claimed contribution is a bit thin---and the main claim is likely wrong. Essentially, there is one clear claimed contribution. That is, a sketching complexity characterization for rank-constrained least squares. The authors claimed that there is no papers discussing this problem without using nuclear norm approximation. This may be true. But the derivation in Appendix F is wrong.  Consider the rank-constrained least squares problem:\n\n            X_r =arg min_{X} ||PX-B||^2,  subject to  rank(X)<=r,\n\nThe problem cannot be solved by first computing unconstrained least squares and then projecting the solution to a rank-r matrix via SVD. The entire proof seems to have been built upon this misconception. For example, in line 870, the authors stated that “X_r is the best rank-r approximation of X_opt” where X_opt is the unconstrained least squares solution. \n\n3.\tFollowing the above comment, the RSVD-LRLS algorithm embedded in the ALS algorithm does not make much sense. The rank-constrained least squares cannot be solved by first doing least squares and then doing truncated SVD. In general, any constrained least squares problem, even with convex constraints, cannot be tackled this way. This procedure does not have any guarantees of the overall ALS algorithm. In general, ALS would need the subproblems to be either solved to a certain accuracy or at least decrease the subproblems cost with a sufficiently large amount to ensure the overall convergence; see [R1] and [R2].\n\n[R1] Razaviyayn, Meisam, Mingyi Hong, and Zhi-Quan Luo. \"A unified convergence analysis of block successive minimization methods for nonsmooth optimization.\" SIAM Journal on Optimization 23.2 (2013): 1126-1153.\n\n[R2] Xu, Yangyang, and Wotao Yin. \"A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion.\" SIAM Journal on imaging sciences 6.3 (2013): 1758-1789.\n\n\n4.\tThe writing of this paper is a bit hard to parse or follow. It seems that too many details are buried in the supplementary material. There is a claim that “One can derive a sketch size bound that is sufficient to get … residual norm for linear least squares with general  (not necessarily rank-based) constraints (detailed in Appendix G). This bound is looser than (3.1), while applicable for general constraints.” But a close inspection on Appendix G finds that the proof there never used the constraint {C}---which is not easy to understand.\n\n5.\tAlso on the representations. The bounds on Tensor sketch is also not easy to understand since R^N-1 appears in the denominator of the bounds. This essentially means that when R and N gets large, the sketching difficulty does not rely on the tensor rank or the tensor order, which is quite counter-intuitive.\n\n6.\tThe authors seemed to have missed many recent developments in stochastic gradient/proximal gradient based tensor decomposition algorithms, but merely focused on sketched least squares. It was shown that if one takes a stochastic gradient view point to compute tensor decomposition, then there is no need to guarantee sketching size/complexity for each subproblem, and thus each subproblem could be updated rather inexactly using lightweight computations. There is a serious lack of discussion on this line of methods.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "n/a", "review_text": "", "extra_scores": {"ethical_concerns": "n/a", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "YTeYmoQ8hZ", "reviewer": "Reviewer_AP2P", "summary": "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.", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper presents an interesting new randomized sketching approach to solve the low rank tensors decomposition problem. The paper is well-written and the results are interesting. A randomized  approach that approximately solves the rank-constrained problem (instead of the relaxed variants) seems novel and will be useful in many applications. Numerical results demonstrate improve performance results using the proposed method. \n\nHowever, I have the following comments about the paper:\n1. Certain aspects of the paper are not clear and can be included to improve the paper. E.g., \nDoes the analysis extend to different dimensions and ranks along different modes? How does the sampling complexity and computational cost compare to others existing randomized methods?  How large in the constant in O(\\epsilon) relative error?\nHow accurately is line 11 in the Algorithm solved. How does the proposed method differ from the other randomized approaches such as the one in [38]. Note that the proposed method will be more expensive (by a factor of mR^n) than some of the other randomized method that solve relaxed versions of the LS problem.\n\n\n2. Numerical results: The numerical experiment results section can be further improved. Computational cost comparison can be provided. Also, results are presented wrt. fitness. It would be interesting to see the approximation error for the factor computed using the different methods (as in solve the ALS completely and report the approximation error for the decomposition).\n\n3. Minor Comments:\ni. Many notations are not defined before usage. Having a small subsection on Notation will be helpful. Readers not familiar with tensor algebra might find it difficult to follow otherwise.\nii. Page 3, line 128, mostly widely --> most widely", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "There does not seem to be a discussion on the limitations and potential negative societal impact of their work in the paper.\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}, {"review_id": "VGpvaNE8qKq", "reviewer": "Reviewer_Vycx", "summary": "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.", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Strengths:\n- The problem is definitely an important one\n- Several contributions across both theoretical and methodological results\n\nLimitations:\n- Very limited experimental evaluation and discussions on real data\n- More extensive intuition behind the RRF algorithm would be necessary\n- No discussion exists regarding what is the real-world challenge being tackled through the use of tensor decomposition on the real image datasets\n\nMy primary concern related to the submission is that the experimental evaluation related to real datasets is not convincing enough. Very limited discussions (if any) exist regarding the Figures 1d and 1e and it's not clear what should the reader consider as the main message being conveyed.\n\n----\nI acknowledge having read the authors' response and my score remains the same; while the authors acknowledge the limitations raised around the real-data experiments, those are not sufficiently addressed based on their response.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "There is no such Section or discussions around that in the paper.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This work contributes sketching-based alternating least squares algorithms for fast and accurate low-rank Tucker tensor decomposition. The main novel contributions are relative-error analyses of the error in the solution of the least squares systems when leverage-score sampling or TensorSketches are used to sketch the systems. The experimental evaluation is somewhat lacking, as uses small and/or synthetic data sets, but the algorithms are shown to result in lower error than previous randomized algorithms for Tucker decompositions. 
Decision: Accept (Poster)