Subquadratic Kronecker Regression with Applications to Tensor Decomposition
Authors: Matthew Fahrbach
Gang Fu
Mehrdad Ghadiri
Venue: NeurIPS 2022 Accept (2022)
Keywords: Linear Regression
Kronecker Product
Leverage Score Sampling
Sketching
Iterative Methods
Tensor Decomposition

Abstract: Kronecker regression is a highly-structured least squares problem $\min_{\mathbf{x}} \lVert \mathbf{K}\mathbf{x} - \mathbf{b} \rVert_{2}^2$, where the design matrix $\mathbf{K} = \mathbf{A}^{(1)} \otimes \cdots \otimes \mathbf{A}^{(N)}$ is a Kronecker product of factor matrices. This regression problem arises in each step of the widely-used alternating least squares (ALS) algorithm for computing the Tucker decomposition of a tensor. We present the first subquadratic-time algorithm for solving Kronecker regression to a $(1+\varepsilon)$-approximation that avoids the exponential term $O(\varepsilon^{-N})$ in the running time. Our techniques combine leverage score sampling and iterative methods. By extending our approach to block-design matrices where one block is a Kronecker product, we also achieve subquadratic-time algorithms for (1) Kronecker ridge regression and (2) updating the factor matrix of a Tucker decomposition in ALS, which is not a pure Kronecker regression problem, thereby improving the running time of all steps of Tucker ALS. We demonstrate the speed and accuracy of this Kronecker regression algorithm on synthetic data and real-world image tensors.

Subquadratic Kronecker Regression with Applications to Tensor Decomposition
Matthew Fahrbach" Google Research aahrbach@goggee. com
Gang Fu Google Research thomasfu@google com
Mehrdad Ghadiri Georgia Tech
ghadiri@gatech. edu
Abstract
Kronecker regression a highly-structured least squares problem minx || b| where the design matrix K A (1) matrices. This regression problem arises in each step of the widely-use alternating least squares (ALS) algorithm for computing the Tucker decomposition of tensor. We present the first subquadratic-tiin to a (1++)-app --aaproxximation that avoids the exponential term O( ((e-N) the running time. Our techniques combine leverage score sampling and iterative methods. By extending our approach block-design matrices where one block is a Kronecker product, we also achieve subquadratic-tine algorithms for (1) Kronecker ridge regression and (2) updating the factor matrices Tucker decomposition in ALS, which is not pure Kronecker regression problem, thereby improving the running time of all steps of Tucker ALS. We demonstrate the speed and accuracy of this Kronecker regression algorithm synthetic data and real-world image tensors.
A (N) is a Kronecker product factor
algorithm for solving Kronecker regression
Introduction
Tensor decomposition has rich multidisciplinary history with countless applications in data mining, machine learning, and signal processing [35, 55, 58, 31]. The most widely-used tensor decom- positions are the CP decomposition and the Tucker decomposition. Similar the singular value decomposition of matrix, both decompositions have natural analogs of low-rank structure. Unlike matrix factorization, however, computing the rank of a tensor and the best rank-one tensor are NP- hard [27]. Therefore, most low-rank tensor decomposition algorithms decide on the rank structure in advance, and then optimize the variables of the decomposition to fit the data. While conceptually simple, this approach is extremely effective in practice for many applications.
The alternating least squares (ALS) algorithm is the main workhorse for low-rank tensor decom- position, e.g., it is the first algorithm mentioned in the MATLAB Tensor Toolbox [7]. For both CP and Tucker decompositions, ALS cyclically optimizes disjoint blocks of variables while keeping all others fixed. As the name suggests, each step solves linear regression problem. The core tensor update step in ALS for Tucker decompositions notoriously expensive but highly structured. In fact, the design matrix of this regression problem is the Kronecker product of the factor matrices of the Tucker decomposition K = A @ ... & A (N), Our work builds on line of Kronecker regression algorithms [17, 18, 47] give the first subquadratic-time algorithm for solving Kronecker regression to (1 + s) apprrxxxiatiion while avoiding an exponential term of O NN) ) in the running time.
We combine leverage sampling, iterative methods, and novel way of multiplying sparsified Kronecker product matrices to fully exploit the Kronecker structure of the design matrix. We also extend approach block-design matrices where one block is Kronecker product, achieving
Authors are listed alphabetically. A preliminary version of this work that facuses on efficient sketching for Tucker decompositions appears in arXiv:2107 10654 [21].
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

subquadratic- time algorithms for (1) Kronecker ridge regression and (2) updating the factor matrix of a Tucker decomposition ALS, which not a pure Kronecker regression problem. Putting everything together, this work improves the running time of all steps of ALS for Tucker decompositions and runs in time that is sublinear in the size of the input tensor, linear in the error parameter & g-1, and subquadratic in the number of columns of the design matrix in each step. Our algorithms support 1.2 regularization in the Tucker loss function, so the decompositions can readily be used in downstream learning tasks, e.g. using the factor matrix rows as embeddings for clustering [67]. Regularization also plays a critical role in the more general tensor completion problem to prevent overfitting when data is missing and has applications in differential privacy [10, 8].
The fastent Kronecker regression algorithm of Diao et al. [18] uses leverage score sampling and achieves the following running times for A (n) R..xR,, R.. with I 2 Rg for all [N], where
R = II Rn and w < 2.373 denotes the matrix multiplication exponent [4]:
1. Ö (nnz( (A + by sampling Ö(R= rows of K by their leverage Scores.
2. Õ (nnz(A (")) + !) Re Re~N) by sampling Ö rows from each factor matrix A (n) and taking the Kronecker product of the sampled factor matrices.
Note that the second approach is linear in R, but the error parameter has exponential cost in the number of factor matrices. In this work, we show that the running time of the first approach can be improved subquadratic in R without increasing the running time dependence on in the dominant term, simultaneously improving on both approaches.
Theorem 1.1. For E [N], let (1) E RR,, In 2 Rn, and b E Rh.., There is (1 + e)- approximation algorithm for solving arg min (rggix((((()))
(1) @
A (N)) blß that runs in time
N
))
Ô nnz((A [
min SC[N]
MM
II RES
Ro, Rs
II RA
(1)
nE[NTyS
where MM(a, b, e) the running time of multiplying an a x matrix with a e matrix.
If we do not use fast matrix multiplication (Gall and Urrutia [24] and Alman and Williams [4]), the last term in (1) is -1R2=-1). which is already an improvement over the standard - running time. With fast matrix multiplication, MM Inc Rr, Rs-1 [n] R.,) ) is subquadratic in R for any nontrivial subset S & {0, [N]]. which is an improvement over ((H===)) x 0 If there exists a "balanced" subset S such that Ines Rr & vR, our running time goes as low as Ö (1..626-11) using [24]. For ease of notation, denote the subquadratic improvement by the constant 8" where 2 mingg[N] MM R,, R. II R.,.
Updating the core tensor in the ALS algorithm for Tucker decomposition is pure Kronecker product regression as described Theorem 1. 1.1, but updating the factor matrices is regression problem of the form arg min KMx b 1 where K is Kronecker product and M is a matrix without any particular structure. We show that such problems can be converted to block regression problems where one the blocks is K. We then develop sublinear time leverage score sampling techniques for these block matrices, which leads to the following theorem that accelerates all of the ALS steps.
Theorem 1.2. There ALS algorithm for 1.2- 1.2-regllrizeed Tucker decompositions that takes tensor X E RRIXX-IINN and returns N factor matrices A (n) E RRax XRx and a core tensor S E RR, x--x Ra. such that each factor matrix and core update is (1 + (l++e)-ap appprooximtaion to the optimum with high probability. The running times of each step are:
Factor matrix A (k); O (1)) R. N2g 22 - R. Rg).
Core tensor 9: (
(n)) + + N2 + +22200122000 122-0-=-1)
= Pm, Ryp = R/Ru. RR, and 0* > a constant derived from fast rectangular matrix multiplication.
where R m=

---

Table 1: Running times of TuckerALS factor matrix and core tensor updates with different Kronecker regression methods. The factor matrices are denoted by A (n) E R The input tensor has size I= ..... IN and the core tensor has size R R Rv. Let 1yk l and RyR R/ Ra. We use w < 2.373 to denote the matrix-multipliationn exponent and the constant 0 > 0 for the optimally balanced fast rectangular matrix multiplication as stated in Theorem 4.5, i.e., R2- minrgiN] MM Ra, R. Iner R.).
Algorithm
Factor matrix A (5)
Core tensor g
Naive
O + IRRR& + RE + IR + 1a I.Rt Ol =-= Rn)+ + R + I[ Rn) + (2)) Ól g- I+R([=== R,) ) + RiEe
O + IR) ORR + +I Ra) --=
This paper (Lemma 4.4) This paper (Theorem 1.2)
Ó(R E E Ō(R²
Diao et al. [18]
For tensors of even modest order, the superlinear term in R is the bottleneck in many applications since R is exponential in the order of the tensor. follows that our improvements are significant in both theory and practice as illustrated in our experiments in Section 6.
1.1 Our Contributions and Techniques
We present several new results about approximate Kronecker regression and the ALS algorithm for Tucker decompositions. Below is a summary of our contributions:
1. Our main technical contribution is the algorithm FastKronnckerreegressiio in Section 4. This Kronecker regression algorithm builds on the block-sketching tools introduced in Sec- tion 3, and combines iterative methods with fast novel Kronecker-matie multiplication for sparse vectors and matrices and fast rectangular matrix multiplication to achieve a running time that is subquudratic in the number of columns in the Kronecker matrix. A key insight is to use the original (non-sketched) Kronecker product as the preconditioner in the Richardson iterations when solving the sketched problem. This, by itself, improves the running time to quadratic. Then to achieve subqudratic running time, we exploit the singular value decomposition of Kronecker products and present novel method for multiplying sparsified Kronecker product matrix (Lemma 4.4 and Theorem 4.5).
2. We generalize our Kronecker regression techniques to work for Kronecker ridge regression and the factor matrix updates in ALS for Tucker decomposition. We show that a factor matrix update is equivalent to solving an cqustiiy-ooss/rained Kronecker regression problem with a low-rank update to the preconditioner in the Richardson iterations. We can implement these new matrix-vector products nearly as fast by using the Woodbury matrix identity. Thus, we provably speed up each step of Tucker ALS, i.c., the core tensor and factor matrices.
3. We give a block sketching toolkit in Section 3 that states we can sketch blocks of matrix by their leverage scores, i.c., their leverage scores in isolation, not with respect to the entire block matrix. This is one of the ways we exploit the Kronecker product structure of the design matrix. This approach can be useful for constructing spectral approximations and for approximately solving block regression problems. One corollary is that we can use the "sketch-ame -nd-solve" method for any ridge regression problem (Corollary 3.5).
4. We compare FastkrroneckerRegressi Kronecker regression task studied in [17, 18] and a subroutine in ALS for computing the Tucker decomposition of various image tensors [44, 50, 51]. Our results demonstrate the importance of reducing the running time dependence on the number of columns in the Kronecker product.
with Diao et al. [18, Algorithm 1] on synthetic
1.2 Related Work
Kronecker Regression. Diao et al. [17] recently gave the first Kronecker regression algorithm based on TensorSketch [53] that is faster than forming the Kronecker product. Diao et al. [18] improved this by removing the dependence on ((nz(bb) from the running time, where b E II-.. is the response vector. Reddy, Song, and Zhang [56] recently initiated the study of dynamic Kronecker
3

---

regression, where the factor matrices (()) undergo updates and the solution vector can be efficiently queried. Marco, Martínez, and Viaña [47] studied the generalized Kronecker regression problem.
Ridge Leverage Scores. Alaoui and Mahoney [3] extended the notion of statistical leverage scores to account for L2 regularization. Sampling from approximate ridge leverage score distributions has since played key role in sparse low-rank matrix approximation [16]. the Nyström method [49], bounding statistical risk in ridge regression [48], and ridge regression [14, 48, 41, 33). Fast recursive algorithms for computing approximate leverage scores [15] and for solving overconstrained least squares [40] are also closely related.
Tensor Decomposition. Cheng et al. [13] and Larsen and Kolda [38] used leverage score sampling to speed up ALS for CP decomposition." Song et al. [59] gave a polynomial time, relative- error approximation algorithm for several low-rank tensor decompositions, which include CP and Tucker. Frandsen and Ge [23] showed that if the tensor has an exact Tucker decomposition, then all local minima are globally optimal. Randomized low-rank Tucker decompositions based on sketching have become increasingly popular, especially in streaming applications: [45, 61, 1 11, 60, 31, 46, 44, 2]. The more general problem of low-rank tensor completion is also fundamental approach for estimating the values of missing data [1, 43, 29, 28, 22]. Fundamental algorithms for tensor complction are based on ALS [68, 25, 42], Riemannian optimization [37, 34, 52], or projected gradient methods [65].
Preliminaries
Notation. The order of tensor is the number of its dimensions. We de-
Algorithm TuckerALS
Input: X E 1x'xx-xXxx (Rg, Ro..... Rx), A
note scalars by normal lowercase let- ters x E IR, vectors by boldface lower- case letters x E R", matrices by bold- face uppercase letters X E Rmxn and higher-order tensors boldface script letters X E xR4xXxx---xx, We
1: Initialize core tensor 9 R2 x Rax....xxx XR.
2: Initialize factors A (n) E RlaxRR for E [N]
3: repeat
4: for 1 to N do
5:
K + A (1) 8- . (a 1) 8 . A
use normal uppercase letters to de- note the size of an index set (e.g.. 7:
6:
B x X (n)
for i i= 1 I, do
[N] = {1, (1. N).. The i-th en- try of a vector x is denoted by x,, the (i. entry of a matrix X by Tay, and the (i, -))th entry of third-
8:
y" arg miny KG (m) y bTll2 + A|ly11
9:
Update factor row a t y
10:
K A @ (2) @ @ A (N)
order tensor X by Zisk.
11:
g" + arg ming Kg vec
2 + A||g|
12: Update core tensor G +- vec 1 (g")
Linear Algebra. Let I, denote the n x n identity matrix and 0mxn de- note the X zero matrix. The trans-
13:
until convergence
14:
return 9, A (1),
A (N)
pose of A E Rmxn is AT, the Moore-Penrose inverse (also called pseudoinverse) A+, and the spectral norm is || ||All The singular value decomposition (SVD) of A is a factorization of the form UEVT, where U E Rmxxm and V R=** are orthogonal matrices, and Σ E mx non-negative real diagonal matrix. The entries (A) of Σ are the singular values of A, and rank(A). The compact SVD is related
the number of non-zero singular values is equal to decomposition where E "x" is diagonal matrix containing the non-zero singular values. The Kronecker product of two matrices E R mm and B E I2200 is denoted by A B E R (m9))
Tensor Products. Fibers of tensor are the vectors we get by fixing all but one index. If is third-order tensor, we denote the column, row, and tube fibers by X-yk- Xisk, and Xij:, respectively. The mode-n unfolding of tensor X E x4xx1xx--IIx is the matrix X(n) E 14
that arranges the mode-n fibers of X as columns of X (n) ordered lexicographically by index. The vectorization of X E R //xx1xx- xxxxxxxxxxxxx xly is the vector vec(X)) R Ia formed by vertically stacking
2 The design matrix in each step of ALS for CP decomposition is a Khatri Rao product, not a Kronecker product. CP decomposition does not suffer from a bottleneck step like ALS for Tucker decomposition since it is sparser decomposition, i.e., CP decomposition does not have core tensor jjust factor matrices.

---

the entries of X ordered lexicographically by index. For example, this transforms X E R into tall vector vec(X) by stacking its columns. We use vec (xx to undo this operation when is clear from context what the shape of the output tensor should be.
The n-mode product of tensor X E RR x1xxxxxxx and matrix A RR. is denoted by y X A where y E RR x-/1x--X1-1x Xicss-xxxxnn---a-- XII This operation multiplies each mode-n fiber of X by the matrix A. This operation is expressed elementwise
(XC x - 1.8
21,=1 Pipig...y Aji,."
The Frobenius norm ||0||p of tensor X is the square root of the sum of the squares of its entries.
Tucker Decomposition. The Tucker decomposition decomposes tensor X R into core lensor G E RR x XRxx-XRN and N factor matrices A A E xR,, Given regularization
parameter À E R20, we compute a Tucker decomposition by minimizing the nonconvex loss function
N
L., A x() = |X 9 X1 A
XN A 18 + A
|SS|
Σ n=1
A
(r)
Entries of the reconstructed tensor X
del
9x
(()) x XN (N) are
R1
Rx
Σ 11
Trira...
(N) OiNN
(2)
TN=1
Equation (2) demonstrates that 5 is the sum of R1 ... RN rank- I tensors. The tuple (R1, Rg...., RN the multilinear rank of the decomposition. The multilinear rank is typically chosen in advance and much smaller than the dimensions ..
Alternating Least Squares. We present TuckerALS in Algorithm 1 and highlight its connections to Kronecker regression. The core tensor update (Lines 10-12) is a ridge regression problem where the design matrix Keore E IR (R RR is Kronecker product of the factor matrices. Each factor matrix update (Lines 5-9) also has Kronecker product structure, but there are additional subspace constraints must account for. We describe these constraints in more detail in Section 5.
3 Row Sampling and Approximate Regression
Here we establish our sketching toolkit. The X-ridge leverage score of the i-th row A E Raxd
e (A) def a, (ATA *I)+
The matrix of cross X-ridge leverage scores is A(ATA + *I)+AA AT. We denote its diagonal by (A) because it contains the X-ridge leverage scores of A. Ridge leverage scores generalize statistical leverage scores in that setting À 0 gives the leverage scores of A. We denote the vector of statistical leverage scores by E(A). If A UEVT is the compact SVD of A, then for all E [n]]. we have
(AA) = E
G((A) d[(A)+X uie.
where (A). It follows that every (A) < since U is orthogonal matrix. We direct the reader to Alaoui and Mahoney [3] or Cohen et al. [15] for further details.
The main results in this paper build on approximate leverage score sampling for block matrices. The X-ridge leverage scores of A E gaxd can be computed by appending VXL to the bottom of A to get A E (n++d)xxd xd and considering the leverage scores of A, so we state the following results in terms of statistical leverage scores without loss of generality.
Definition 3.1. For any A E Raxd, the vector Ê(AA) E I a 3-overestimate for the leverage score distribution of A if, for all i [n], it satisfies
(A( [[((A)]11
2 B ((A)
(((A)
[((Aa)]i rank(A) (A)
5

---

Next we describe the approximate leverage score sampling algorithm in Woodruff [64, Section 2.4]. The core idea here is that if we sample ((d/3) 3) rows and reweight them appropriately, this smaller sketched matrix can be used instead of A give provable guarantees for many problems.
Definition 3.2 (Leverage score sampling). Let A E Rn and p E [o, 1 be 3-overestimate for the leverage score distribution of A such that l!pl|, = 1. CampleRows(A,.. 9, p) denotes the following procedure. Initialize sketch matrix xsxn. For each row of S, independently and with replacement, select an index E [n] with probability Pj and set 8üj = 1/ /P5s. Return sketch S.
The main result in this section is that we can choose to sketch single block of a matrix by the leverage scores of that block in isolation. This sketched submatrix can then be used with the other (non-sketched) block to give a spectral approximation the original matrix or for approximate linear regression. The notation A & B is the Loewner order and means B A is positive semidefinite.
Lemma 3.3. Let A = [A1; A2] be vertically stacked with A I xd and A₂ E R Let E [0, 1) be a 3-overestimate for the leverage score distribution of A If 14hdl((2/0)/((6e)) the sketch S returned by SampleRows A 11 8, p) guarantees, with probability at least 6, that
(1 =)AAA< (SA1) + 3 (1
Lemma 3.4 (Approximate block regression). Consider the problem arg min || - blß2 where
A = [A A2] and b [b1: b2] are vertically stacked and A E Ruidd A E Rraad b E Rm, b E 22 Let p E [0. 11 be a B-overestimate for the leverage score distribution of A Let > 1680d 1680d (400d) (Bs) and let S be the ouiput of SampleRows( 81, S.).P). If
(ANEERR(SSSAA (IS(AiX- bi)l A2x b2 12
x arg min
then, with probability at least 9/10, we have
AA bli2 5 (1+
min Ax bll
We defer the proofs of these results Appendix A. The key idea behind Lemma 3.4 is that leverage scores do not increase if rows are appended to the matrix. This then allows us to prove sketched submatrix version of Drineas al. [19, Lemma 8] for approximate matrix multiplication and satisfy the structural conditions for approximate least squares in Drineas et al. [20]. One consequence is that we can "sketch and solve" ridge regression, which was shown in [63, Theorem 1] and [6, Theorem 2].
Corollary 3.5. For any A Igad b Rd X N 0, consider
arg minnmI Ax bl12 + x11x112)
XE3d
Let E 10m be a (--oversstimate for the leverage scores of A and 2 1680d (n0d))(4403))))) (Bee). IfS the output of "SampleRows A, 8, P)) then, with probability at least 9/10, the sketched solution
x arg min (asimiiinnS(aaCCC b)ll 2
XERd
gives (1 + -ppppopoxinaaiion to the original problem.
Remark 3.6. The success probability of the sketch be boosted from 9/10 to 1 - by sampling factor of (log(1/8)) more rows. See the discussion in Chen and Price [12, Section 2] about matrix concentration bounds for more details.
4 Kronecker Regression
Now describe the key ingredients that allow to design an approximate Kronecker regression algorithm whose running time is subquadratic in the number of columns in the design matrix.
1. The leverage score distribution of a Kronecker product matrix K (1) Ø ... @ A (N) is product distribution of the leverage score distributions of its factor matrices. Therefore, we can sample rows of K from ((K) with replacement in ((V) time after preprocessing step.

---

2. The normal matrix KTK + XI in the ridge regression problem miny| | Kx bll2 + >||x 2 is [speccrlaa approximation of the sketched matrix (SK) AI by Lemma 3.3. Thus we can use Richardson iteration with (KTK + XI)+ as the preconditioner to solve the sketched instance, which guarantees (1 + -pprrxxxxxiitoo Using (KTK + XT)+ as the preconditioner allows to heavily exploit the Kronecker structure with fast matrix-vector multiplications.
3. At this point, Kronecker matrix-vector multiplications factor matrices into two groups by their number of columns and use our novel way of multiplying sparsified Kronecker product matrices as well as fast rectangular matrix multiplication to get subquadratic running time.
still the bottleneck,
partition the
This first result shows how X-ridge leverage scores of Kronecker product matrix decompose accord- ing the SVDs of its factor matrices. All missing proofs in this section are deferred Appendix B.
Lemma 4.1. Let K A A (2) @.... A (N), where each factor matrix A (m) R.xxBs, Let (11, 22...., ..). be the natural row indexing of K by its factors. Let the factor SVDs be A (n)
(())) For any X 2 O, the X-ridge leverage scores of K are
D- of (A
N II
(N)....)) (N) (K)
Σ
(5)
n= of (A()) +
tET
11
where the sum is over T [81]
[2g]
x
[x For statistical leverage scores, this simplifies
to (i (K) = II L,, (A (n)..
This proof repeatedly uses the mixed-product property for Kronecker products and the definition of X-ridge leverage scores Equation (3).
4.1 Iterative Methods
Now we state a result for the convergence rate of preconditioned Richardson iteration [57], which uses the notation ||x|| >
Lemma 4.2 (Preconditioned Richardson iteration).
Let M be any matrix such that ATA M M+ (xx)x()) -Abb). Then,
k. ATA for some k 2 1. Let (k++1) = x (1)
xx x" I (1- 1/k) (0) - I'I
where x" arg aargminnggs
| - AA- bll²
Remark 4.3. The ridge regression algorithm in Chowdhury et al. [14] is also based on sketching and preconditioned Richardson iteration. They consider short and wide matrices where d >> n and use the sketched normal matrix as the preconditioner to solve the original problem. One of our main technical contributions is to use the original normal matrix as the preconditioner solve the sketched problem. Reversing this is advantageous because computing the pseduoinverse and matrix-vector products with the original Kronecker matrix substantially less expensive due to its Kronecker structure. However, this still motivates the need for faster Kronecker matrix-vector multiplications.
4.2 Fast Kronecker-Matrix Multiplication
The next result is simple but useful observation about extracting the rightmost factor matrix from the Kronecker product and recursively computing new less expensive Kronecker-mand multiplication.
Lemma 4.4. Let (n) RI,.J,, for n E [N), and E RJI-J There is an algorithm KronMatMul ([A (1) ...., A A B)) that computes (1 in O(K) EM J1 In IN) time.
(1) 6 (2) ... 68 A R((I--xxxX
The following theorem is more sophisticated. We write the statement in terms of rectangular matrix multiplication time (M b, c), which is the time to multiply an a b matrix by e matrix.
Theorem 4.5, Let (n) E IR. Rm, for n E [N], I I ... IN. R = R1 ... Rv, E R cE R and S E I x be a diagonal matrix with ((R=-1) nonzeros. The vectors
(A1 2 @ Ab and S(A1 @ .... 8 AN)c can be computed in time Ó [[N MM NET Ru, Re IIInT Rn)).
7

---

Algorithm 2 FastikroneckeRRegressionn
Input:
Factor matrices A A(") E RRaxRR response vector b E Rh L2 regularization strength X,
enror E, failure probability
Set R < R2 R2 ... R...
R
2: for 1 to N do
3:
Compute a spectral approximation Á (n) with ÖR,, rows by Lemma 3.3 such that
(n) (m)
z (m)]
(n) < (1 log(1 + e/4) (n) (
(6)
4: Compute À (m)T (m) and the SVD of À ( A (z)(()()) ((") V(()))
5: Compute (1 + log(1 + ( leverage scores E( (n)) using Lemma B.4 by
applying random Johnson- Lindenssn-Lnnddssrruusss projection
6: Initialize product distribution data structure P sample indices from (E(A ... _((A
" ())
7: Set D + (Σ()) (11) @ ..... (NN) E(N) + Rl)
8: Let M = (V(1) ... & D( @.... V(N))T
9: Set -- I1680 (40R) ln(1/)//e]
10: Set S + SampleRowa s,
11: Let = SK and = Sb On
12: Initialize
13: repeat
14: X x (1 (1- + xx - kTb) using fast Kronecker-matrix multiplication
until convergence
16: return
The core idea behind Theorem 4.5 that the factor matrices can be partitioned into two groups to achieve a good "column-
"columm-producct" product" balance, i.e., minrg[N] max (II Rn, Ingr R., is close to
VR. Then we use the fact that nnz(S) = Ö(R=1)) with a sparsity- KronMatMul to solve each part of this partition separately, and combine them with fast rectangular matrix multiplication. If achieve perfect balance, the rumning time is Ö D[R.6626--1) using results of Gall and Urrutia [24], which are explained in detail in van den Brand and Nanongkai [62, Appendix CJ. If one of these two factor matrix groups has at most 0.9 of the "columm-produe mass," the running time is -
4.3 Main Algorithm
We are now ready to present our main algorithm for solving approximate Kronecker regression.
Theorem 4.6. For any Kronecker product matrix K A (1) @ & (N) R Ix Ri Rx
Rd IN A 0. 104, and 8 > FastKronnecer/legresssoo returns R....r
in
( + 2) + minsgIM) MM [Inss Ru, Re- IINENNIS Rs) Ra)).
time such that, with probability at least 8,
||xx bll2 + A||x S (1+ min Kx - bl2 A||x||
We defer the proof to Appendix B.2 and sketch how the ideas in Algorithm 2 come together. First, we do not compute the pseudoinverse K- but instead use iterative Richardson iteration (Lemma 4.2), which allows avoid e(H''=1) running time. This technique by itself, however, only allows us to reduce the running time to Ō since all of the matrix-vector products (e.g., KTb, Kx, and multiplication against M+) naively take 2( time. To achieve subquadratic time, we need three more ideas: (1) compute an approximate SVD of each Gram matrix A (n) A (T)) in order to construct the decomposed preconditioner M (2) use fast Kronecker-veest vector multiplication (e.g., Lemma 4.4) to exploit the Kronecker structure of the decomposed preconditioner; (3) noting that Lemma 4.4 for the Kronecker-vector products KTb and KT (Kx) is insufficient because the intermediate vectors can be large, we develop novel multiplication algorithm in Theorem 4.5 that fully exploits the sparsity, Kronecker structure, and fast rectangular matrix multiplication of Gall and Urrutia [24].

---

Applications to Low-Rank Tucker Decomposition
Now we apply our fast Kronecker regression algorithm to TuckerALS and prove Theorem 1.2. We list the running times of different factor matrix and core update algorithms in Table I, and we analyze these subroutines in Appendix C.3.
Core Tensor Update. The core update running time in Theorem 1.2 is a direct consequence of algorithm for fast Kronecker regression in Theorem 4.6. The only difference is that avoid recomputing the SVD and Gram matrix of each factor since these are computed at the end of each factor matrix update and stored for future use.
Factor Matrix Update. The factor matrix updates require more work because of the G (n) term
in Line 8 of TuckerALS. To overcome this, we substitute variables and recast each factor update as an ggnstity- Kronecker regression problem with an appended low-rank block account for the L2 regularization of the original variables. To support this new low-rank block, we use the Woodbury matrix identity to extend the technique of using Richardson iterations with fast Kronecker matrix-vecter multiplication for solving sketched regression instances.
The next result formalizes this substitution and reduces the problem block Kronecker regression with subspace constraint. This result relies on the fact that the least squares solution ||xx 2 with minimum norm is M*z.
Lemma 5.1. Let A E RR M E Rm xd E R", and X 2 0. For any ridge regression problem of the form arg minlER' (|| bll2 + 21)x we can solve
arg min Az bll2 + A|| Nx
Zupt
where N = Im - MM and return vector M Zopt instead.
Proof. Let z= = Mx RTT For any x E R is in the column space of M and hence orthogonal to
any vector in the left null space of M. Therefore, we can optimize over z E RT subject to Nz instead because for any x E Rd, NMx = (Im MM MM = (M - M)x = O. Using this substitution, we can also replace the term \||x||2 by A||M+ 21|2 because for any z, the least squares solution to z Mx with minimum norm is M 2 [54].
To solve this constrained regression problem, we can add scaled version of the constraint matrix N as block to the approximate regression problem and take the projection of the resulting solution. regression). Let M Raxd N E mmd b e R", b12 to (1 + =)paprooxinttion it suffices to solve
Lemma 5.2 (Approximate
and 0 < E < 1/3. To solve minNx=o MM
M
b
min
XER
V
0
-a)po///apoimmattion with w (1 12/e) |MN+|
Letting z = GT y y in Line 8 of TuckerALS and modifying FastKroneckerleerresiinn to support additional low-rank updates to the preconditioner, we get the FastFactorVatrixUpdas algorithm, presented as Algorithm in Appendix C.2. The analysis is similar the proofs Theorem 4.6. The factor matrix updates benefit in the same way as before from fast Kronecker matrix-vecton products, and new low-rank block updates are supported via the Woodbury identity. We defer the proofs of the next two results to Appendix C.
Theorem 5.3. For any À 2 0, e (0, 1/3). and > 0. the Fattaaaaatattttatiiiipppaaa algorithm updates A Blaxity in TuckerALS with (1 + =appproxximittoonn with probability at least 1 - 8. in time
l1/(1/)) + I R. + R
Corollary 5.4. FastFactorVatriitppattt updates A (k) R Ru in Ó E- log(1/8)+ (og(1/)))))))
4h
IRYN Rr + Iife time, where 0* > 0 is the optimally balanced MM exponent in Theorem 4.5.

---

30
30
1200
Time (SeconoS)
20
20
20
800
20
10
10
400
0-
5200
10000
15000
3000
10003
15000
3000
10000
15000)
3000
10000
15000
0
1250018
Favikrrr www.kkrrebbbellliieooo berlllggessiii
Figure 1: Running times of Kronecker regression algorithms with design matrix of size n² x
6 Experiments
All experiments were run using NumPy [26] with an Intel Xeon W-2135 processor (8.25MB cache, 3.70 GHz) and 128GB of RAM. The battkroneckerkegresssioo obbbased ALS experiments for low-rank Tucker decomposition on image tensors are deferred to Appendix D.2. All of our code available at https httpt/ttttttttttt////// http//farraacssuuuuuuddddee haapppuadrrattc--oonneeee----eee
Kronecker regression. We build on the numerical experiments in [17, 18] for Kronecker regression
that use two random factor matrices. We generate matrices A (1) A (2) Rnxd where each entry is drawn i.i.d. from the normal distribution (11, 0.001) and compare several algorithms for solving minx]](A (1) & A x 1 1 + +| 1| as we increase d. The running times are plotted in Figure 1.
The algorithms we compare are: (I) baseline that solves the normal equation (KTK+ A KTb and fully exploits the Kronecker structure of KTK before calling mp. linalg pinv (); (2) an enhanced baseline that combines the SVDs of A (m) with Lemma 4.4, e.g., KKwoMaaMualalla (((U)))) (U(2)) b), using only Kronecker-vector products; (3) the sketching algorithm of Diao et [18, Algorithm 1]; and (4) our FastKrroncckerBegrees on algorithm in Algorithm 2. For both sketching algorithms, we use = 0.1 and 0.01. We reduce the number of row samples in both algorithms by a = 10-5 so that the algorithms more practical and comparable to the earlier experiments in [17, 18). Lastly, we set À = 10 We discuss additional parameter choice details and the full results in Appendix D.1.
The running times in Figure demonstrate several different behaviors. The naive baseline quickly becomes impractical for moderately large values especially since it is exact method. The runtimes of the sketching algorithms are nearly-independent of Diao al. [18] works well for small d, deteriorates tremendously as d grows because
d. KronMatMul is competitive for n S 10
computes (SK) TSK + AI)+ E Te xd² and cannot exploit the Kronecker structure K, which takes (a))) time. PastKrronackerRRgrrsssii
the other hand, runs in O(d*) time because uses quadratic-time Cronecker-vector products in each Richardson iteration step (Line 14).
Table 2: Kronecker regression losses for = 64. OPT denotes the loss of the KronMatMul algorithm, DJSSW19 is Diao et al. [18, Algorithm 1], and Algorithm is FaskKroneckerrRgresssii record the relative error of each algorithm and the number of rows sampled from A @ A (2),
We also
7
OPT
Algorithm 2
Approx
DJSSW19
Approx
Rows sampled (%)
1024
0.031
0.032
1.051
0.035
1.138
0.0370
2048
0.123
0.126
1,026
1.577
12,792
0.0093
4096
0.507
0.520
1.026
275.566
543.776
0.0023
8192
2.073
2.136
1,030
333.430
160.809
0.0006
16384
8.238
8.608
1.045
546391.728
66329.791
0.0001
These experiments also show that combining sketching with iterative methods can give better sketch efficiency. Table compares the loss of [18, Algorithm 1] and FastKronecke ceereeneerrerrssssisc. to an exact baseline OPT for d 64. Both algorithms use the exact same sketch SK for each value of Our algorithm uses the original (KTК + AI) + as a preconditioner to solve the sketched problem, whereas Diao et al. [18, Algorithm 1] computes (SK))SK
(SK) exactly and becomes numerically unstable for 2 2048 when d E (16, 32, 64). This raises the question about how to combine sketched information with the original data to achieve more efficient algorithms, even when solving sketched instances. We leave this question of sketch efficiency interesting future work.
10

---

References
[1] Evrim Acar, Daniel M Dunlavy, Tamara G Kolda, and Morten Mørup. Scalable tensor factoriza- tions for incomplete data. Chemometrics and Intelligent Laboratory Systems, 106(():41-56. 2011.
[2] Salman Ahmadi-Asl, Stanislav Abukhovich, Maame G 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.
[3] Ahmed Alaoui and Michael W. Mahoney. Fast randomized kernel ridge regression with statistical guarantees. In Advances in Neural Information Processing Systems, volume 28, 2015.
[4] Josh Alman and Virginia Vassilevska Williams. A refined laser method and faster matrix multiplication. In Proceedings of the 2021 ACM-SIAM Symposium on Discrete Algorithms (SODA), pages 522-539. SIAM, 2021.
[5] Rosa I Arriaga and Santosh Vempala. An algorithmic theory of learning: Robust concepts and random projection. Machine learning, 63(2):161-182, 2006.
[6] Haim Avron, Kenneth L. Clarkson, and David P. Woodruff. Sharper Bounds for Regularized Data Fitting. In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques (APPROX/RANDOM 2017), volume 81, pages 27:1-27:22, 2017.
[7] Brett W. Bader and Tamara G. Kolda. Tensor toolbox for MATLAB, version 3.2.1. https:
/ / www. www.ortooorrooooxooooooooooooo.. org/, 2021.
[8] Raghavendran Balu and Teddy Furon. Differentially private matrix factorization using sketching techniques. In Proceedings of the 4th ACM Workshop on Information Hiding and Multimedia Security, pages 57-62, 2016.
[9] Casey Battaglino, Grey Ballard, and Tamara G Kolda. A practical randomized cp tensor decomposition. SIAM Journal on Matrix Analysis and Applications, 39(2):876-901, 2018.
[10] Kamalika Chaudhuri, Claire Monteleoni, and Anand D Sarwate. Differentially private empirical risk minimization. Journal of Machine Learning Research, 12(3), 2011.
[11] 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.
[12] Xue Chen and Eric Price. Active regression via linear-sample sparsification. Conference on Learning Theory, pages 663-695. PMLR, 2019.
[13] 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.
[14] Agniva Chowdbury, Jiasen Yang, and Petros Drineas. An iterative, sketching- based framework for ridge regression. In International Conference on Machine Learning, pages 989-998. PMLR, 2018.
[15] Michael B Cohen, Yin Tat Lee, Cameron Musco, Christopher Musco, Richard Peng, and Aaron Sidford. Uniform sampling for matrix approximation. In Proceedings of the 2015 Conference on Innovations in Theoretical Computer Science, pages 181-190, 2015.
[16] Michael B Cohen, Cameron Musco, and Christopher Musco. Input sparsity time low-rank approximation via ridge leverage score sampling. In Proceedings of the Twenty-Eight Annual ACM-SIAM Symposium on Discrete Algorithms, pages 1758-1777. SIAM, 2017.
[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 Statistics, pages 1299-1308. PMLR, 2018
11

---

[18] Huaian Diao, Rajesh Jayaram, Zhao Song, Wen Sun, and David Woodruff. Optimal sketching for kronecker product regression and low rank approximation. Advances in neural information processing systems, 32, 2019.
[19] Petros Drineas, Ravi Kannan, and Michael W Mahoney. Fast Monte Carlo algorithms for matrices I: Approximating matrix multiplication. SIAM Journal on Computing, 36( 36(1):132-157, 2006.
[20] Petros Drineas, Michael W. Mahoney, Shan Muthukrishnan, and Tamás Sarlós. Faster least squares approximation. Numerische Mathematik, 117(2):219-25 249, 2011.
[21] Matthew Fahrbach, Mehrdad Ghadiri, and Thomas Fu. Fast low-rank tensor decomposition by ridge leverage score sampling. arXiv preprint 1rXiv:2107.10644 2021.
[22] Marko Filipović and Ante Jukić. Tucker factorization with missing data with application to low-n-rank tensor completion. Mulridimensional systems and signal processing, 26(3):677-692,
2015.
[23] Abraham Frandsen and Rong Ge. Optimization landscape of tucker decomposition, Mathemati- cal Programming, pages 1-26, 2020.
[24] François Le Gall and Florent Urrutia. Improved rectangular matrix multiplication using powers of the coppersmith-wing Symposium on Discrete Algorithms, pages 1029-1046. SIAM, 2018
tensor. In Proceedings of the Twenty-Ninth Annual ACM-SIAM
[25] Lars Grasedyck, Melanie Kluge, and Sebastian Kramer. Variants of alternating least squares tensor completion in the tensor train format. SIAM Journal on Scientific Computing, 37(5): A2424-A2450, 2015.
[26] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Vir- tanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Hal- dane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre GGrarddMarrchat Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant. Array programming with NumPy. Nature, 385-7825) September 2020. doi: 10.1038/<41586-222 000-1-88--000--9---- URL https://dost org/ 10 00.00////41888
[27] Christopher Hillar and Lek- lek-Heng Lim. Most tensor problems are NP-hard. Journal of the ACM (JACM), 60(6):1-39, 2013.
[28] Prateek Jain and Sewoong Oh. Provable tensor factorization with missing data. Advances in Neural Information Processing Systems, 27, 2014.
[29] Prateek Jain, Praneeth Netrapalli, and Sujay Sanghavi. Low-rank matrix completion using alternating minimization. Proceedings of the forty-fifth annual ACM symposium on Theory of computing, pages 665-674, 2013.
[30] M James. The generalised inverse. The Mathematical Gazette, 62(420): 1978.
[31] Jun-Gi Jang and Kang. Fast and memory-efficient tucker decomposition for answering diverse time range queries. In Proceedings of the ACM SIGKDD Conference on Knowledge Discovery & Data Mining, pages 725-735, 2021.
[32] William B Johnson and Joram Lindenstrauss, Extensions of lipschitz mappings into hilbert space 26. Contemporary mathematics, 26, 1984.
[33] Praneeth Kacham and David Woodruff. Sketching algorithms and lower bounds for ridge regression. In International Conference Machine Learning, pages 10539-105556. PMLR, 2022.
[34] Hiroyuki Kasai and Bamdev Mishra. Low-rank tensor completion: a Riemannian manifold preconditioning approach. In International Conference on Machine Learning, pages 1012-1021. PMLR, 2016.
12

---

[35] Tamara G. Kolda and Brett W. Bader. Tensor decompositions and applications. SIAM Review, 51(3):55-500, 2009.
[36] Jean Kossaifi, Yannis Panagakis, Anima Anandkumar, and Maja Pantic. Tensorly: Tensor learning in Python. Journal of Machine Learning Research (JMLR), 20(26), 2019.
[37] Daniel Kressner, Michael Steinlechner, and Bart Vandereycken. Low-rank tensor completion by Riemannian optimization. BIT Numerical Mathematicx, 54(2):447-468, 2014.
[38] Brett W Larsen and Tamara G Kolda. Practical leverage-based sampling for low-rank tensor decomposition. SIAM Journal on Matrix Analysis and Applications, 43(3):1488 1517, 2022.
[39] Brett W. Larsen and Tamara G. Kolda, Sketching matrix least squares via leverage scores estimates. arXiv preprint 1036: 2022.
[40] Mu Li, Gary Miller, and Richard Peng. Iterative row sampling. In 2013 IEEE 54th Annual Symposium Foundations of Computer Science, pages 127-136. IEEE, 2013.
[41] Zhu Li, Jean-Francois Ton, Dino Oglic, and Dino Sejdinovic. Towards unified analysis of random Fourier features. In International Conference on Machine Learning, pages 3905-3914. PMLR, 2019.
[42] Allen Liu and Ankur Moitra. Tensor completion made practical. In Advances in Neural Information Processing Systems, volume 33. pages 11905-18916, 2020.
[43] Ji Liu, Przemyslaw Musialski, Peter Wonka, and Jieping Ye. Tensor completion for estimating missing values in visual data. IEEE transactions on pattern analysis and machine intelligence, 35(1):208-22) 200,22202222222 2012.
[44] Linjian Ma and Edgar Solomonik. Fast and accurate randomized algorithms for low-rank tensor decompositions. Advances in Neural Information Processing Systems, 34, 2021.
[45] Osman Asif Malik and Stephen Becker. Low-rank Tucker decomposition of large tensors using TensorSketch. Advances in Neural Information Processing Systems, 31: 10106-100666, 2018.
[46] Osman Asif Malik and Stephen Becker. A sampling-basee method for tensor ring decomposition. In International Conference on Machine Learning, pages 7400-7411. PMLR, 2021.
[47] Ana Marco, José-Javier Martínez, and Raquel Viaña. Least squares problems involving gen- cralized kronecker products and application to bivariate polynomial regression. Numerical Algorithms, 82(1):21-39, 2019.
[48] Shannon McCurdy. Ridge regression and provable deterministic ridge leverage score sampling. In Advances in Neural Information Processing Systems, volume 31, 2018.
[49] Cameron Musco and Christopher Musco. Recursive sampling for the Nyström method. In Advances in Neural Information Processing Systems, volume 30, 2017.
[50] 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.
[51] Sameer A Nene, Shree K Nayar, Hiroshi Murase, al. Columbia object image library (coil-20). 1996.
[52] Madhay Nimishakavi, Pratik Kumar Jawanpuria, and Bamdev Mishra. A dual framework for low-rank tensor completion. In Advances in Neural Information Processing Systems, volume 31, 2018.
[53] Rasmus Pagh. Compressed matrix multiplication. ACM Transactions on Computation Theory (TOCT), 5(3):1177, 2013.
[54] M. Planitz. Inconsistent systems of linear equations. The Mathematical Gazetlee 63(425): 181-185, 1979.
13

---

[55] Stephan Rabanser, Oleksandr Shchur, and Stephan Günnemann. Introduction to tensor de- compositions and their applications machine learning. arXiv preprint arXiv:1711.10711, 2017.
[56] Aravind Reddy, Zhao Song, and Lichen Zhang. Dynamic tensor product regression. arXiv preprint arXiv:2210.03961, 2022.
[57] Yousef Saad. Iterative methods for sparse linear systems. SIAM, 2003.
[58] 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((3):3551-3582 2017.
[59] Zhao Song, David P Woodruff, and Peilin Zhong. Relative enror tensor low rank approximation. In Proceedings of the Thirtieth Annual ACM-SIAM Symposium on Discrete Algorithms, pages SIAM, 2019.
[60] 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, 2(4):1123-115) 2020.
[61] Abraham Traore, Maxime Berar, and Alain Rakotomamoniy. Singleshot: scalable tucker tensor decomposition. In Advances in Neural Information Processing Systems, volume 32, 2019.
[62] Jan van den Brand and Danupon Nanongkai. Dynamic approximate shortest paths and beyond: Subquadratic and worst-case update time. In 2019 JEEE 60th Annual Symposium on Foundations of Computer Science (FOCS), pages 436 455, IEEE, 2019.
[63] Shusen Wang. Alex Gittens, and Michael W Mahoney. Sketched ridge regression: Optimization perspective, statistical perspective, and model averaging. In International Conference Machine Learning, pages 3608-3616. PMLR, 2017.
[64] David P. Woodruff. Sketching as a Tool for Numerical Linear Algebra. 2014.
[65] Rose Yu and Yan Liu. Learning from multiway data: Simple and efficient tensor regression. In International Conference Machine Learning, pages 373-381. PMLR, 2016.
[66] Huamin Zhang and Feng Ding. On the kronecker products and their applications. Journal of Applied Mathematics, 2013, 2013.
[67] Guoxu Zhou, Andrzej Cichocki, and Shengli Xie. Decomposition of big tensors with low multilinear rank. arXiv preprint arXiv: 1885, 2014.
[68] Hua Zhou, Lexin Li, and Hongtu Zhu. Tensor regression with applications in neuroimaging data analysis. Journal of the American Statistical Association, 108(502):540-555 2013.
14


Reviews:
[{"review_id": "yj3pnVjIiYF", "reviewer": "Reviewer_FzUX", "summary": "The work considered the fundamental problem of Kronecker regression and proposed the first sketching-based subquadratic algorithm to obtain a (1+\\epsilon)-approximation solution. Both rigorous proofs and empirical validations were shown. The block-sketching methodology can be of independent interests itself.", "questions": "In the experiment part, was fast matrix multiplication methods used? How does the runtime complexity compared to the theoretical results in Theorem 1?", "limitations": "Yes.", "rating": 7, "confidence": 2, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The work focused on an accelerated solution for a fundamental problem which makes the work very meaningful. The technique itself is quite interesting. The math is well presented and validated by empirical results.", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "x6uq2daMUxv", "reviewer": "Reviewer_Ubgq", "summary": "The authors propose the first subquadratic-time algorithm for solving Kronecker regression. Kronecker regression is a structured problem where the design matrix is a Kronecker product of multiple matrices. This problem appears when trying to compute the Tucker decomposition of a Tensor. The use score sampling and iterative methods to achieve a complexity reduction for Kronecker regression problems. Then, they extend their approach to block design matrices. They use synthetic and real examples to demonstrate the computational advantages of their approach. ", "questions": " Can you also provide examples on real data? \nWhat is nnz on page two? This is now explained.\nPerhaps the paper would read better if the theorems can be pushed forward outside of the introduction. \nP5L201 what is s (the scalar s).\n", "limitations": "The limitations of the proposed work are properly detailed in section F. ", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper reads well, and the English level is satisfactory. I am not familiar with how important the addressed problem is, but the method seems novel with solid justification, both theoretical and empirical. The runtime improvement seems substantial, and the relative error of the algorithm on synthetic data seems low. The related works are reviewed well and proper background to the problem is provided. ", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "sn9TFe9_6aH", "reviewer": "Reviewer_yEQL", "summary": "In this paper, the authors study fast sketching based algorithms for Kronecker Product Regression and show how this algorithm can be used to speedup Alternating Least Squares (ALS) algorithm for computing Tucker Decomposition. The main improvement over previous work of Diao et. al. [16], which introduced leverage score sampling for Kronecker Product Regression seems to be recognizing that the matrix $M = (K^{\\intercal}K +\\lambda  I$) can be used as a preconditioner for solving the sketched ridge regression using Richardson's Iteration and showing a way to quickly compute $(K^{\\intercal}K +\\lambda  I)^+v$ for a vector $v$ using fast matrix multiplication by appropriately dividing the matrices in the Kronecker Product into two disjoint subsets.. Although updating the core tensor in Tucker Decomposition requires one to solve a simple Kronecker Product Ridge Regression problem, the factor matrix update requires one to solve a more complicated Kronecker Product Ridge Regression with a subspace constraint. The authors show that their techniques can be used to solve this version of the Kronecker Product Ridge Regression and therefore all the steps in ALS algorithm for computing Tucker Decomposition are made faster.", "questions": "- Can you clarify my doubts about the input sparsity claims in Theorem 1.1?\n- I am confused about where Theorem 4.5 is used in $\\texttt{FastKroneckerRegression}$. Isn't the only product of a vector with a kronecker product in line 13? And the vector there isn't sparse. Am I missing anything in regards to where it is used?\n- If I understand correctly, the main speedup over [16] is obtained by not computing the product $\\tilde{K}^{\\intercal} \\times \\tilde{K}$ and instead using the Richardson's iteration which using Fast Kronecker Product you can perform quickly? Can you explain how the term $\\tilde{O}(\\text{MM}(\\Pi_{n \\in T} R_n, R\\varepsilon^{-1}, \\Pi_{n \\notin T}R_n))$ arises in line 13?\n- How does solving intermediate problems of ALS to a small approximation affect the overall convergence? I'd think it shouldn't be affected too much but do you have any small experiments where you observe anything significant?\n- Please let me know if I am misunderstanding or missing anything. I can take a closer look at sections in the appendix if they are explained there better.", "limitations": "NA", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- Shows how speeding up Kronecker Product Ridge regression helps in running ALS for Tucker Decomposition. This wasn't observed in the previous works [15,16].\n- Improving the running times over previous work by observing that solving the sketched problem using Richardson Iteration with $(K^{\\intercal}K + \\lambda I)^{+}$ is preconditioner leads to faster running times.\n\nWeaknesses:\n- The paper doesn't have any new contributions for sketching literature. Note that sketching a submatrix $A_1$ by leverage score sampling rows of $A_1$ using leverage scores with respect to $A_1$ obviously works because of the property that leverage scores only decrease when the new rows are added. This observation is not new and can be seen in the paper \"Sharper Bounds for Regularized Data Fitting\".\n- It is not clear how the authors claim input sparsity running times for Kronecker Product Ridge Regression when in their algorithm, they compute singular value decompositions of $A_1,\\ldots, A_N$. If one doesn't compute singular value decompositions, how is $(K^{\\intercal}K + \\lambda I)^{+}$ used as a preconditioner? Am I missing something here? I understand input sparsity time algorithms are probably not necessary for intermediate problems of ALS for tucker decomposition as one cannot guarantee that the intermediate matrices will be sparse but how is the claim in Theorem 1.1 correct?\n- Writing could be a bit better by introducing Tucker Decomposition in the introduction. I found it a bit hard to follow without knowing what the core tensor and factor matrices are. This is probably not that big an issue as the target audience are aware of the terminology.\n\nSome typos:\n- In definition 3.2, you also need to state that $\\||p\\||_1 = 1$ -- not a big problem\n- Missing $\\||\\cdot\\||_2$ at a few places like in Corollary 3.5\n\n———————————————————————————————————————\nUpdate after author responses:\n- Authors answered my questions and made a few changes to the manuscript which address most of my concerns.\n- Although the paper does not contribute to the sketching literature as I said earlier, there are a lot of details and new techniques required to apply the existing techniques to the important problem of Tucker Decomposition which lead to my acceptance decision.", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "Rnd-pUB1gny", "reviewer": "Reviewer_djSP", "summary": "This paper provides faster algorithms for solving the $\\ell_2$ - Kronecker product regression problem, which is a least squares regression problem where the design matrix is a Kronecker product of several smaller matrices with dimensions $I_n \\times R_n, n \\in [N]$. They also extend their techniques to provide faster algorithms for Kronecker ridge regression, and speed up the Alternating Least Squares (ALS) algorithm for computing the Tucker decomposition of a tensor. Prior work (Diao et.al, 2019) for Kronecker product regression cared about the case where $R_n \\gg I_n$ for all $n \\in [N]$ and this work improves the runtime when $R_n$ might be comparable to $I_n$ and justify that this case is important through their experiments.\n\nReference:\nHuaian Diao, Rajesh Jayaram, Zhao Song, Wen Sun, and David Woodruff. Optimal sketching for kronecker product regression and low rank approximation. Advances in neural information processing systems, 32, 2019.", "questions": "I briefly went over the Diao et.al [16] paper and their runtimes seem to have a factor of $1/\\epsilon^2$ (see footnote 9 on page 6) rather than $1/\\epsilon$ as claimed on Line 46 of this paper?\n\nEditing suggestions:  \nIt would be better if you define $\\omega$ (along with a reference and an approximate value) before using it on page 2.\n\nLine 93, 260 - “pseduoinverse” -> “pseudoinverse”\n\n", "limitations": "Primarily theoretical work.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Originality: I am not an expert in this area so I am not entirely sure about other related work and how novel this work's contributions are compared to other papers in the area.\n\nQuality: The submission is technically sound. All claims are well-supported with proofs.\n\nClarity: The submission is clearly written and well-organized.\n\nSignificance: Yes, the problems being studied are significant. Kronecker product regression is a well-studied problem with several applications.\n", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept