Efficient Leverage Score Sampling for Tensor Train Decomposition
Authors: Vivek Bharadwaj
Beheshteh T. Rakhshan
Osman Asif Malik
Guillaume Rabusseau
Venue: NeurIPS 2024 poster (2024)
Keywords: leverage score sampling
tensor train decomposition
alternating least square
Primary area: optimization
TL;DR: We propose an efficient leverage score based randomized alternating least squares algorithm for tensor train decomposition.

Abstract: Tensor Train~(TT) decomposition is widely used in the machine learning and quantum physics communities as a popular tool to efficiently compress high-dimensional tensor data. In this paper, we propose an efficient algorithm to accelerate computing the TT decomposition with the Alternating Least Squares (ALS) algorithm relying on exact leverage scores sampling. For this purpose, we propose a data structure that allows us to efficiently sample from the tensor with time complexity logarithmic in the product of the tensor dimensions. Our contribution specifically leverages the canonical form of the TT decomposition. By maintaining the canonical form through each iteration of ALS, we can efficiently compute (and sample from) the leverage scores, thus achieving significant speed-up in solving each sketched least-square problem. Experiments on synthetic and real data on dense and sparse tensors demonstrate that our method outperforms SVD-based and ALS-based algorithms.

Efficient Leverage Score Sampling for Tensor Train Decomposition
Vivek Bharadwaj Beheshteh Rakhshan Osman Asif Malik", Guillaume Rabusseau
EElectrical Engineering and Computer Science Department, UC Berkeley *Encube Technologies
Mila & DIRO, Université de Montréal *CIFAR AI Chair
Abstract
Tensor Train (TT) decomposition is widely used in the machine learning and quan- tum physics communities as popular tool to efficiently compress high-dimensional tensor data. In this paper, we propose an efficient algorithm to accelerate comput- ing the TT decomposition with the Alternating Least Squares (ALS) algorithm relying on exact leverage scores sampling. For this purpose, we propose a data structure that allows us efficiently sample from the tensor with time complexity logarithmic the lensor size. Our contribution specifically leverages the canonical form of the TT decomposition. By maintaining the canonical form through each iteration of ALS, we can efficiently compute (and sample from) the leverage scores, thus achieving significant speed-up in solving each sketched least-square problem. Experiments on synthetic and real data on dense and sparse tensors demonstrate that our method outperforms baseb and ALS-based algorithms.
1 Introduction
Tensor decomposition methods have recently found numerous applications in machine learning. Their ability to perform operations efficiently on very high- nigh-dimensional tensors makes them suitable for data science and machine learning problems. For example, they have been used for neuro- ceuro-imaging. and signal processing [Zhou al., (2013] (Sidiropoulos al. [201 Cichocki and Phan 2009), super- vised learning [Novikov et al. 2016] Stouudeemeee and Schwab, 2016], feature extraction [Bengua jtt al. 2015) and scaling up Gaussian processes [zmailov et al. | 2018]. The most popular decom- positions are the CANDECOMPPPRAAAL (CP), Tucker [Hitchcock]] 1927 Tucker 1966), Tensor (011T and Tensor Ring (TR) Zhao al. [2016] decompositions. However, finding a rank-r CP decomposition is NP-hard Kolda and Bader 2009 Hillar and Lim| [013] and in
Train (TT) Oseledets
(TR)
the Tucker decomposition, the number of parameters grows exponentially with the order of the tensor. While the number of parameters is linear in the tensor's 's order for both TT and TR decompositions, TR is known to have numerical stability issues. The TT decomposition, in contrast, can be used to represent tensor in compressed format where the number of parameters scales linearly with the order of a tensor. Additionally, finding a good approximation to the best rank-r TT decomposition is feasible which lets the user perform various linear algebra operations efficiently and compute the TT decomposition with stable algorithms.
Due to the high-dimensional nature of tensors, designing efficient algorithms for computing the TT decomposition is crucial. A popular method for computing the TT decomposition of an N- dimensional tensor X is the TT-SVD algorithm [Oseledets] 2011 which uses sequence of singular values decompositions on the tensor unfoldings to produce the TT representation in a single pass.
Equal contribution
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

Since TT-SVD requires performing SVDs of unfoldings of X' its cost is exponential in N. Alternating Least Square (ALS) another popular approach [Holtz et al. 2012] find the TT approximation. Starting with crude guess, each iteration of ALS involves solving sequence of least squares problems. While ALS is the workhorse algorithm in many tensor decomposition problems, the computational cost is still exponential in the order of a tensor (N), since each iteration requires solving least squares problems involving unfoldings of X. These issues have led to the search for alternatives based on randomization and sampling techniques. A cheaper alternative to the TT-SVD with strong accuracy guarantees can be implemented by replacing the exact singular value decomposition (SVD) with well-studied randomized counterpart [Halko et al. 2011 Huber et al., Randomized variants of the TT-ALS approach have received little attention. In Chen et al. 2023 the authors propose randomized ALS algorithm that uses TensorSketch [Pham and Pagh 2013] in each iteration. In this work, we also propose a novel randomized variant of the TT-ALS algorithm that relies on exact leverage score sampling. Notably, the sketch size in TensorSketch TT-ALS[Chen et al. [2023] has an exponential dependence on the tensor dimension / whereas our algorithm avoids any dependence of the sketch size on
2017].
Our Contributions. In this paper, we propose a new sampling-based ALS approach to compute the TT decomposition: rTT-ALS. By using exact leverage score sampling, we are able to significantly reduce the size each ALS least squares problem while providing strong guarantees on the approxi- mation error. At the core of rTT-ALS, we leverage the TT canonical form to efficiently compute the exact leverage scores and speed up the solutions of least square problems in each iteration of ALS. To the best of our knowledge, rTT-ALS is the first efficient TT decomposition by the ALS algorithm which relies on leverage scores sampling. We provide experiments on synthetic and real massive sparse and dense tensors showing that rTT-ALS can achieve up 26x speed-up compared to its non-randomized counterpart with little to no loss in accuracy.
Our core contribution is the following theorem, which shows that can efficiently compute subspace embedding of a left-orthogonal chain of TT tensor cores by efficiently sampling according to their squared row norms (for technical definitions and details see subsection 3.I.
Theorem 1.1 (Row- (Row-norm-ss norm-squared sampling for 3D core chains). Let A1, A₁ be sequence of 3D tensors, A E R x Ra (with Ro = 1). Assume that the left- eeft-matricization of each core is orthogonal. Let Acj be the IT I x R3 matrix obtained by unfolding the contraction of the tensor chain AI.....A; Then there exists data structure to randomly sample rows from Agj according to the distribution of its squared row norms with the following properties:
1. The data structure has construction time O (E RR.. t
When R
R
R2 and I I structure linear in the sizes of the input cores.
1y- the runtime is O(jRR)). The space overhead of the data
2. The data structure produces a single row sample from Agj according to the distribution of its squared row norms time EE log (Tx R Rk) RR)). When all ranks Rg and physical dimensions Ig are equal, this complexity is O(j log I).
We highlight that the runtime required to construct the data structure asymptotically identical to the runtime required to compute the canonical form the tensor train subchain, i.c., Açj. by successive QR decompositions. This implies that the data structure construction and subsequent updates of the data structure do not increase the asymptotic complexity of our method.
2 Related work
Randomized algorithms and leverage score sampling-based methods [Mahoney et al. 2011 Woodruff et al. 2014 Drineas et al. 2006a] have been used widely in a large body of research particularly
in tensor decomposition problems
the past two decades (Malik and Becker 2021 Larsen and
Kolda 2022 Fahrbach et al. 2022 just to name a few.
Cheng et al. 2016 propose SPALS, the first ALS-based algorithm relying on leverage score sampling for the CP decomposition. Their proposed method reduces the size of the least squares problem in each iteration of ALS with sub-linear cost per iteration in the number of entries of the input

---

tensor. Larsen and Kolda [222] extends this method by combining repeated sampled rows in deterministic and random sampling fashion. However. both of these methods use leverage score approximations and therefore require a number of samples which is exponential in the number of tensor modes in order to achieve relative-error performance guarantees. [Malik] [2022] proposes a method which avoids this exponential dependency on the number of tensor modes by using higher- quality leverage score estimates for the CP decomposition. The method is further improved by MMalik et al. [022] to use exact rather than approximate leverage scores which is applicable for arbitrary tensor decompositions. Recently, [Bharadwaj et al. 2023] provided a novel data structure to efficiently sample from the exact distribution of the factor matrices' leverage scores in the Khatri-Rao product with time complexity logarithmic in the tensor size, leading to further improvements on the work in [Malik al. 2022]. Moreover, MMalik and Becker 2021] proposed an ALS-based algorithm for finding the TR decomposition using leverage scores approximation. However, the runtime of the order of a tensor. The sampler we propose in
their method has exponential dependency this paper is built on the work by |Bharadwaj et al. (2023). extending it to the TT decomposition and leveraging the canonical form for further speed-up.
There are also a variety of non-ALS-basee randomized algorithms for computing the TT decomposi- tion. (Huber al. 2017] leverages randomized SVD for the TT decomposition which accelerates the classical TT-S VD algorithm proposed by (Oseledets 201T]. To handle situations where the exact TT rank is unknown. [Che and Wei 2019] propose adaptive randomized algorithm that can achieve near optimal TT approximation. Yu et al. 2023] present a method leveraging randomized block Krylov subspace iteration for computing TT approximations. Most of the algorithms for TT decomposition are based on the randomized SVD for matrices introduced by (Halko et al.. 2011] In the quantum physics community, the ALS algorithm is widely used for finding TT decomposition and often yields more accurate results than TT-SVD. The randomized TT-SV method struggles to scale for high-order tensors as it requires generating a random Gaussian matrix at each step that can only handle small-order tensors. More closely related to our work are those using sketching and sampling in each iteration of ALS approximate the TT decomposition. Recently, [Chen et al.. 2023) introduced an algorithm that employs TensorSketch [Pham and Pagh, 2013 in each iteration of regularized ALS approach for TT decomposition. However, the sketch size has an exponential dependency to tensor dimension. By contrast, our proposed algorithm requires the sketch size with no dependence on the tensor dimension and it depends only on the column size of the design matrix and and 8 parameters.
3 Preliminaries
We use capital letters A to denote matrices and script characters AA to denote multidimensional arrays. We use Matlab notation for slices of matrices and tensors. We use the tuple notation to indicate the position of entries of arrays. For example, A(i1, 22, d3) indicates the (i1, i2, i3)-th element of A. A [i, and A [: refer to the i-th row and column of A, respectively; for three-dimensional tensor
A E R Ra Ji x Rs the matrix AA [., i E IR Ri x Ra is the i-th lateral slice of A. For a positive integer n, we use [n] to denote the set of integers from 1 to n. For i [4],.... . AN E [IN] the notation
:]
1) 1 1 1j will be helpful for tensor unfoldings. We use Ø and O to denote the Kronecker and Khatri-Rao products, respectively (see definitions in Appendix A}. We use Ia to denote the d identity matrix, A for the transpose of A, + for the pseudo-inverse of A, I|e for the Frobenius norm and I1- for the Euclidean norm of vector, We use Ô to indicate the presence of multiplicative terms polylogarithmic in R and 1/5.
in 1 + Mn
3.1 Tensor Train Decomposition
Let X E R be an N-dimensional array. A rank (Ri...., Rx-1) tensor train (TT) de- composition of a tensor X E 1-1xxxxxx factorizes it into the product of N third-order tensors An E R RRR--xL, xR, for n E [N] (with Ro Rx 1):
(iN)= X(1,, ((N))
Σ
II
(n(n--1,
To. n=1
3

---

for all 11 E [41] ... AN E [Ix], where each ranges from to Rx. A tensor network representation of a TT decomposition is shown in Figure We call A A2, AN core tensors and we use TT((A,)==)) to denote TT tensor with factors AI, ... An.
R4
A
Figure 1: Tensor Train decomposition of 5-dimensional tensor in tensor network notation.
Definition 3.1. The mode-n unfolding of tensor X E eR4x 1R1xxII is the matrix ((n) x2lxxxlzp 1
defined dlement-wise by ((n) (ins in ... in in-14 i. = X(i1,...,ix)
def
As special case, we denote the left (resp. right) matricization of a 3-dimensional tensor A E 1xxxxx by AL (A) E ERllaass and A A(1) R1illlss.
Given TT decomposition TT((A,) and an index j. we will often use the left-chain Acj E RIK[]] Rx and right-chain Azj E R,, xIIx unfoldings obtained by matricizing the contraction of all cores the left and on the right side of the j-th core. Formally,
j-1 II (r(--1, and Asi (ri,izs)
Acj ((<s)
Σ
Σ
II
Au Ar(
1yt1..,NN k=j+1
where icj = i1 ..-1 and
رمg iy+1...dy. We also use A#1 def
Acj & Ab3 E
Rllez 1x R to denote the unfolding of the contraction of all cores except the j-th one.
We conclude by introducing the canonical form of the TT decomposition Holtz et al. 2012 Evenbly 2018 2022] which will be central to the design of our algorithm.
Definition 3.2. A TT decomposition (((A)))
R xxxxxxIxx in a canonical format with respect to fixed index j E [N] if A A Tn. for all < j. and A A j (see Figure 2
T Ta-- for all
R1
Ro
Ro
R1
Ao
Figure 2: Orthonormal TT decomposition. The cores at the left side of A3 are left-orthonormal and the cores at the right ottht-orthonormal.
Note that any TT decomposition can efficiently be converted to canonical form w.r.t. any index j [N] by performing series of QR decompositions the core tensors HHlltz al. 2012 Evenbly 2018
3.2
Alternating Least Squares with Tensor Train Structure.
The decomposition problem consists in finding a low-rank approximation (AA()))) of given tensor X: argmin
".TAA,...,,,)))). where X' is the target tensor with this a non-convex optimization problem, the popular alternating least-squares (ALS) approach can be used to find an approximate solution (Kolda and Bader 2209). Fixing all cores except the j-th one, the low rank approximation problem can be reformulated as linear least squares problem:
A....n || dimensions I .... IN. AN.....AA
argmin A, ((Ac) & AT,) (A)) XTI
The ALS approach finds an approximate solution by keeping all cores fixed and solving for the j-th one. Then repeat this procedure multiple times for each j [N] until some convergence

---

criteria is met. While ALS is the workhorse algorithm in many tensor decomposition problems, the computational cost is still O(1*) ) for combine ALS with core orthogonalization to efficiently compute the exact leverage scores. This will also lead to a stable algorithm for computing TT. To compute the orthogonalized TT approximation, we start with crude TT decomposition in canonical form (see Definition (.2) where all cores except the first one are right-orthonormal. After optimizing the first core, a QR decomposition is performed and the non-orthonormal part is merged into the next core. This procedure repeats until reaching the right side of the decomposition. The same procedure is then repeated from the right until reaching the left side (see the tensor network illustration in Appendix A.I). Even after computing the canonical form, which makes the design matrix orthonormal for each linear least squares problem, we still need to multiply the matricized tensor by the chain of TT cores in canonical form. Given that the tensor may have hundreds of millions of nonzeros, this multiplication is a significant computational bottleneck: without sketching, the runtime for this matrix multiplication scales as ((nnz(X)) R²), where nnz(2') is the number of nonzeros in the tensor. Sketching allows us to select only a subset of rows from the design matrix and the corresponding subset of rows from the matricized tensor, reducing the cost to O(N R²).. This approach also leads to providing computational benefits for computing the leverage scores and an efficient sampling scheme which will be discussed in Section
IN I. In this work, to reduce this cost, we will
3.3 Sketching and Leverage Score Sampling
There exists vast literature on randomized algorithms Mahoneyy et al. 2011 Woodruff et al. 2014 to solve the over-determined least squares problem min, || - bilr where A R >> R. Regardless of the structure of both A and b, solving this least-squares problem costs O To reduce this cost, we can randomly select rows of A and by proposing sketching operator S with
< I. Therefore, instead of solving the original least squares problem, we consider solving the downsampled version of the form min, |SAr Sbll r. where S E RJ and reduce the cost to O()R²) The goal is to find a "good" sketch S to approximate the solution of the least squares problem at each step of the ALS algorithm. When each entry of S is selected according the rows of A leverage scores, strong guarantees can be obtained for the solution of the downsampled problem.
Definition 3.3. (Leverage scores) Suppose A E RIXR with I > R. The i-th leverage score of the matrix A is defined
4(A) = Alt. (AT A) : for i [1].
(2)
Definition 3.4. (Leverage score sampling) Let I and E 10, be a probability distribution [((A)
vector with entries Pt- rank(A]' where rank(A) = 4(A. Assume $1....,; $5 are drawn i.i.d according to the probabilities "1, .P. The random matrix E RJxx defined element-wise by S(j,t) JJp, 8 and 0 otherwise is called leverage score sampling matrix for A.
The following result is well-known and appeared in several works; see, e.g., DDrineas et al. 2006b), [Drineas et al. 2008). [Drineas et al. 2011]. [Larsen and Kolda 2022). We borrow the form presented in Maalk 2022].
Theorem 3.5. (Guarantees for Leverage Score Sampling) Suppose A E RIXR, Let S R be the leverage score sampling matrix defined in 3.4 For any 5 E (0, 1). if = O(R2/E6), then a) min, || bllg with probability 8 8.
i min, ||Ar SbIl, satisfies ||Ar* bl12
Computing leverage scores in Definition costs (11²²) and as costly as directly solving the original least squares problem. In the following section, we will show that the leverage scores can be computed much more efficiently when A is the matrix appearing in the TT-ALS algorithm in canonical form..
3.3 requires computing the pseudo-inverse of A, which
4 Sampling-hased Tensor Train Decomposition
In this section, we show how to efficiently sample rows of A A Ac, SAT, and Xt)) in Equation i according to the exact leverage scores distribution. In doing so, we will also present the sketch of the proof of Theorem [.] (which closely mirrors that of Bharadwaj et al. 2023] with key modifications required to adapt the procedure to tensor core chain).
5

---

For each row ر7 .......-11j: ...NN of A Equation 2 gives
1ip (A+1) = / ( 11A A
Computing +b def (A A is the main computational bottleneck in finding the leverage scores of A#3_ Malik et al. [2022] proposed an algorithm to compute in time O(NTR2 + R²). In this paper, we leverage the fact that when the TT tensor is in canonical form w.r.t mode j. A# is orthogonal, and thus Φ Ra. Therefore, computing + is free of cost. By maintaining the canonical form of the TT tensor throughout the ALS algorithm, we can sketch the least square problems from the leverage score distributions with almost computational overhead. We now explain how to efficiently sample rows of A#1 from the leverage scores distribution.
4.1 Efficient Core Chain Leverage Score Sampling
As discussed above, when the TT tensor is in canonical form, the leverage score of row ji is given by less (A+)) (i (;
Leveraging the Kronecker structure of A*) = Acj A3j-
one can easily show that (, (A*j) = (.,((A))) (Acs) 1iz, (A5,). Sampling from the leverage scores distributions thus boils down to sampling rows of Acj and A with probability proportional to their squared row norms (due the orthogonality of Acj and Azj inherited from the canonical form). Without loss of generality, we detail the sampling procedure for Ac; (the difference between Ag; and Acj amounts to reindexing). The sampling procedure for A>5 will be the same and straightforward.
Let 81 E [4]], 8, E [ be random variables such that the multi-index 85j = $1... follows the leverage score distribution of Asj. Since TT((A,,) ((A))) i) is in canonical form w.r.t. ++ 1, Agi an orthonormal matrix, hence &j. the corresponding row of Agz:
&j is selected with probability proportional to the squared norm of
1 1 (Asylst....sy.
p = $1 83 s;) = R,
A -.glsllAsss........... (...s1...s...1))
(4)
Our sampling procedure will draw a lateral slice from each core starting from Aj and ending with A.. corresponding to a single row Agj. Suppose we have drawn sk+1..... 88, for some < j. To sam-
p(se..)))) ple the k-th index, we need to compute the conditional probability p(sN .......... 85) p(as 15))
The following lemma shows that this can be done efficiently by leveraging the underlying TT structure.
Lemma 4.1 (Conditional distribution for 8x). Consider the events 83 S..... 8 8841, which we abbreviate as 8> 8>k- Then
p(šk = 8k 8 sya) x Tr [T A li sk. A sky HyH
where H = Act1 [t, 8k+1, Aí s1.
The proof is given in Appendix Intuitively, H>* acts as a "history matrix" conditioning on 8>k, while the trace operation corresponds marginalization over 8ck. Unfortunately, updating H through matrix multiplication as each index is selected still requires time O(R²) (assuming R1 = R2 R). order further improve the runtime and reach the quadratic complexity in R claimed in Theorem 1.1 we make the following observation: let y E R be the probability vector for the leverage score distribution of Agj. Then Equation can be rewritten in vector form as q = R2 (Agy [:, 1)2 ... + Agi [: R,}²) Here, the square of each column vector is
elementwise operation. Observe that each Ac, [:, r is probability vector (positive entries summing one) due the orthonormality of Acj- Hence q is mixture distribution. To sample from q, thus suffices to select a single column â of Ag; uniformly at random and restrict the sampling procedure to Ac₁ [;, 2 More formally, let be uniformly distributed over [Ry] and let i i,,... follow the conditional distributions defined by
plie to 16+11 t4+1 .,, -,, = ||AA [;, t. hs
(5)
where fos ARt1 (1 .... ...... Ag r.) We have the following result.

---

Lemma 4.2. For any choice of ",,.... fix s, =t,, 85-1 = Sk Ey. After marginalizing over it the conditional distribution of is satisfies p(fx- = is [i*k = (sk) = p(s: = 8k مم* = 8>k)-
As a consequence, the joint random variable (f1,..., Eg) follows the desired squared row-norm distribution of As, after marginalizing over f. The proof appears in Appendix B Notice that the "history matrix" " Н> has been replaced by a vector b This vector can updated by matrix- vector multiplication, yielding a reduced sampling complexity with only quadratic dependency on R.
Our final improvement is to show that each sample from the distribution in Equation (5) can be drawn in time sublinear in the dimension l (after appropriate preprocessing). Letting A be the left unfolding of An. one can check that
R R
p(ia- = Ls |5az = byk, f (=))=
Σ 1=0
(AF Rk-1 + 2,:] hys)
The probability of selecting the slice 8% is thus the sum of Rx-1 consecutive entries from the probability vector (A h>*)² )2 As result, we can sample & by first sampling an index in the range given by (Ag-h then performing integer division by R to obtain the corresponding slice index 4. The advantage here lies in an efficient data structure for sampling from the weight vector (At- h>*)², given by the following lemma:
Lemma 4.3 (Bharadwaj et al. 2003]] Adapted). Given a matrix A E RRXRR, there exists data structure with construction time O(IR)) and space usage O(IR) such that, given any vector h E RR single sample from the un-normalized distribution of weights (A hi2 can be drawn in time O (R))
The adaptation of this lemma is given in Ap- pendix [B.3 Lemma 4..3 enables us efficiently draw samples according to the distribution in Equation 0 and therefore gives us procedure to sample from the entire core chain. Construct- ing the data structure above for each matrix A 1 S k S j. costs O (RE) with a linear space overhead in the input core sizes. Draw- ing sample from the k-th data structure re-
Algorithm 1 ConssuuCChhnnSSmpplee... An)
1: for k = 1..N do
2: Z BuildSamplle(AAA)
Algorithm 2 ChainSampleLLeft((..))
1: for d = 1..J do f Uniform-sample(l1..)
2:
quires time O OR, log R Rx.). Summing
3:
h :
up this runtime over 1 < k < j gives the stated complexity in Theorem 1.] [ [nd summarize the procedures to efficiently draw J samples from left-orthogonal core chain. The construction procedure builds a set of data struc- tures Zμ of the form given by Lemma 4.3 on the left-matricization of each tensor core. For each of J rows to draw, the sampling algorithm selects column i uniformly random from the left matricization. It then initializes the history vector / and successively samples indices iy-1,.... according to the conditional distribution, updating the history vector at each step. Appendix B.4 provides a rigorous proof of the correctness of the procedure sketched in this section.
4:
for k = j...1 do
5:
is == RowSample Zx. ht/ R&-1
6:
h h= h- A- [i ta = (fx)as)
7:
8:
return ti....ty
While our procedure shares similarities with the Khatri-Rao product leverage score sampler, signifi- cant adaptations are required to sample from a tensor train core chain. The factors of a Khatri-Rao product can be sampled in any order, since the Khatri-Rao product of several matrices is commutative up to permutation of its rows. By contrast, our sampling procedure requires us to sample from core A, down A1, since Lemma 4.1 exploits the orhhogottyyyy of the each core in its derivation. Starting the sampling procedure at A, leads to a "history matrix" to keep track of prior draws instead of the vector that would arise starting from core Ai. Here, our second innovation of sampling a column uniformly at random is required to bring down the overall sampling complexity. We can now state the following guarantee for Randomizze--T--ALS (rTT-ALS) applying the data structure in Theorem I.I
The proof is given in Appendix B.5
7

---

Corollary 4.4. (rTT-ALS) For any E. E (0. 1) the sampling procedure proposed above guarantees that with J = Ó(R² / /c8) samples per least-square problem, we have
[At)(A))T)) -x11 S +1+
min (Aj)(z)
(A))
with probability (1 6), where A₃ is the solution of the sketched least-squares problem. for all least-squares solve. The efficient sampling procedure of Theorem 1.7 brings the overall complexity to # Rd 2j=1 N Noggg + 11) where "#it" is the number of ALS iterations.
Algorithms I andd 2 refer to procedures "BuildSampler" and "RowSample" that were first used to sample from the Khatri-Rao product. The k-th BuildSampler data structure creates full binary tree that truncated to log 1x Rx-1/ / R& levels, each caching an Rg x R& matrix containing information from AL. To draw a sample, the RowSample procedure executes a random walk from the root to the leaves requiring O (RR) work at each internal node. The sampler performs matrix- multiplication with the cached data at each internal node and compares the output value to a threshold, using the comparison to branch either left or right in the random walk. Assuming R1
R₃ 1,. the storage cost of each sampler is OTT (computed by multiplying the
and 1 = 1 matrix size stored at each node by the maximum node count). By multiplying the tree depth by the matrix- -vector multiplication cost at each node, we get runtime cost O log I) to draw one sample from Af. For the motivation behind the procedure, details of its correctness, and pseudocode, we refer the reader to the original work (Bharadwaj et al. 2023
Experiments
All experiments were conducted on CPU nodes of the NERSC Perlmutter, an HPE Cray EX supercomputer, and the Mila Quebec AI Institute compute cluster. Our code is available at https:///ithub. cwm/ wttp:/raaaaaabb//rrraaaa
this section, we demonstrate the effectiveness of the proposed rTT-ALS on two types of tensors: (i) synthetic and real dense datasets and (ii) real sparse datasets. We use the fit as evaluation metric (higher is better); fit(X², = 1 - | X11F// /||-2/// where X is the TT approximation and X' is the target tensor. The goal of the dense tensor experiments is to show that rTT-ALS has better time complexity than TT-ALS and TT-SVD while matching rTT-SVD in terms of fit. The sparse tensor experiments show that SVD- decompositions cannot handle high-order (sparse) tensors. We compare rTT-ALS with the classical TT-ALS. The runtime improvements are most significant for large sparse tensors. Figure
accuracy (y-axis, higher is better) against ALS iteration time, for rTT-ALS versus non-randomized ALS. The speedup per iteration can be as high as 26x for lower ranks. Partieularly, for the NELL--2 tensor, the plot shows that accuracy within three significant figures of non-randdd ALS was achieved roughly 3-4x faster than an optimized non-randomized ALS baseline.
5.1 Decompasition of Synthetic and Real Dense Datasets
We compare ITT-ALS to three other methods; TT-SVD Oseledets 2011]. Randomized TT-SVD (FTT-VVD) [Hber et al., 2017) and TT-ALS (Holtz et al. 2012]. We use TensorLy (Kossaifi et al., 2019] for SVD-based methods and our own implementation for deterministic TT-ALS. For simplicity, we set Ri =... RN -1 R for all experiments. For all algorithms, we illustrate the quality of performance by fit and runtime.
Synthetic Data Experiments. For the synthetic data experiment, we generate random tensors of size I x I for E {100, 500} and of TT rank R 20 (by drawing each core's components ii.d. from a standard normal distribution). A small Gaussian noise with mean zero and standard deviation of 10-6 6 is added to each entry of the resulting tensor. We then run the four methods to find rank R = approximation of the target tensor. ALS- b-ase methods are initialized using their SVD-based counterpart (TT-ALS with the output of TT-SVD and rTT-ALS with the output of rTT-SVD) and are run for 15 iterations. The sample count for rTT-ALS fixed to = 5000 for all values of 7. The average fit over 5 trials for all four algorithms are reported as function of the dimension in Figure 3 rTT-ALS about 2x faster than TT-ALS and 3x faster than TT-SVD for I = 500. Although r TVD is the fastest method, it achieves poor performance in terms of fit.
8

---

svd svd
Itavd nt-svd
2.09
24
(proposal)
rth al5 (proposal)
0.0
Drsendion Size 21
Dimordinn Size
Figure 3: Fit (left) for J 5000 and running time (right) averaged over 5 trials for the synthetic data experiment.
Table Decomposition results for real datasets with J 2000 and the target rank R Time is in seconds.
Pavia Uni.
Tabby Cat
MNIST
DC Mall
Method
Fit
Time
Fit
Time
Fit
Time
Fit
Time
TT-ALS
0.61
4.16
0.65
44.570
0.46
8.29
0.59
21.86
TTT-ALS (proposal)
0.60
0.82
0.65
7.360
0.45
2.20
0.59
2.81
TT-SVD
0.61
6.65
0.65
136.189
0.46
17.19
0.59
41.45
TTTSVD
0.61
0.33
0.65
4.285
0.46
0.65
0.59
0.46
Real Data Experiments. For the real data experiment, we consider four real images and video datasets (more details about datasets are given in Appendix C): (i) Pavia University is a hyper-spectral image dataset size (610 x 340 x 103). (ii) DC Mall is also dataset of hyper-spectral images of size (1280 x 307 x 191). Both datasets are the image height and width, and the third dimension the number of spectral bands, (iii) the MNIST dataset of size (60000 x 28 28), and iv) Tabby Cat is the three- diree-diinsiinall tensor of size (720 x 1280 X 286) which contains grayscale videos of a man sitting on a park bench and a cat, respectively. The first two dimensions are frame height and width, and the third dimension is the number of frames. For all datasets, the preprocessing step is done by tensorizing data tensors into digheo-dimeeniooal tensors. Table illlsttrate the results for a single trial when R = 5. For all datasets keep the sample count fixed to J 2000. Similarly to the synthetic data experiments, rTT-ALS is faster than TT-ALS and TT-SVD (up to 10x faster than TT-ALS).
three-dimensional tensors where the first two dimensions
5.2 Approximate Sparse Tensor Train Decomposition
We next apply rTT-ALS to three large sparse tensors from FROSTT [Smith al. 2017]. Table 2 gives the fits achieved by our method to decompose these tensors. The largest of these tensors, NELL-2, has around 77 million nonzero entries with mode sizes in the tens of thousands. Fits for sparse tensor decomposition are typically low, but the factors of the resulting decomposition have
Eoron
NELL-2
0.18
0.07
0.000
0.025
0.16
0.12
0.05
2.020
0.10
2.013
0.05
0.07
2.000
0.02
Randen LS5
Emart
Easet Bandom L55
2.005
Emant Randen 1S5
0.00
0.00
102
25c
350
400
2.000
300 122
Time D1
Time 30
Time (1
Figure 4: Fit as a function of time for three FROSTT tensors, R 6, J 216 for ITT-ALS. Thick lines are averages of 5 fit-time traces, shown by thin dotted lines.

---

Table 2: Average fits and speedup, 216 for randomized algorithms, 40 ALS iterations. The speedup is the average per-iteration runtime for single exact ALS sweep divided by the average time for single randomized sweep.
Uber
Enron
NELL-2
R
rTT-ALS
TT-ALS
Speedup
rTT-ALS
TT-ALS
Speedup
ITT-ALS
TT-ALS
Speedup
4
0.1332
0.1334
4.0x
0.0498
0.0507
17.8x
0.0213
0.0214
26.0x
6
0.1505
0.1510
3.5x
0.0594
0.0611
12.4x
0.0265
0.0269
22.8x
8
0,1646
0.1654
3.0x
0,0669
0.0711
10.5x
0.0311
0.0317
22.2x
10
0.1747
0.1760
2.4x
0.0728
0.0771
8.5x
0.0350
0.0359
20.5x
12
0.1828
0.1846
1.5x
0.0810
0.0856
7.4x
0.0382
0.0394
15.8x
6.40
31966
31344
Senplan))
Sunglee
Tuntir
Don
Figure 5: Final fit of sparse tensor decomposition for varying sample counts. Each boxplot reports statistics for 5 trials. The blue dashed lines show the fit for non-randomized ALS.
successfully been mined for patterns [Larsen and Kolda 2022 For these experiments, we chose all decomposition ranks equal with R1
Ry R. and tested over range of values for R.
The fits produced by rTT-ALS match those produced by the non-randomized ALS method up to variation in the third significant figure for Uber and NELL-2, with slightly higher errors on the Enron tensor. We kept the sample count for our randomized algorithms fixed J 20 throughout this experiment. As result, the gap between the fit of the randomized and exact methods grows as the target rank increases, which our theory predicts.
Table 2 also reports the average speedup per ALS sweep of rTT-ALS over the exact algorithm. On the NELL-2 sparse tensor with target rank 12, the non-randomized ALS algorithm requires an average of 29.4 seconds per ALS sweep, while rTT-ALS requires only 1.87 seconds. Figure shows that our method makes faster progress than its non randomized counterpart across all three tensors. Because we could not find well-documented, nigh-performance library for sparse tensor train decomposition, we wrote fast multithreaded implementation C++, which serves as the baseline method in these figures and tables (the code will be made publicly available).
Figure [ruure]]show shows the impact of varying the sample count the final fit. We find modest increases in accuracy for both Uber and NELL-2 as the sample count increases by a factor of (starting from J 215). Increasing J has smaller impact for the Enron tensor, which is generally more difficult to decompose beginning with i.i.d. random factor initialization [Larsen and Kolda 2022
Conclusion
We proposed a sampling-basee AL.S method leveraging an efficient data structure to sample from the exact leverage scores. More precisely, show that by exploiting the canonical form of the TT decomposition, leverage scores can be computed efficiently for all the least squares problems of AL.S. We provide strong theoretical guarantees for the proposed data structure. Experiments massive dense and sparse tensors confirm the theoretical results. The sampling algorithm we proposed could be extended to more general tree-based tensor network structures, leveraging canonical forms in similar spirit to ITT-ALS.
10

---

Acknowledgments and Disclaimers
V. Bharadwaj was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Department of Energy Computational Science Graduate Fellowship under Award Number DE-SC0022158. This research used resources of the National Energy Research Scientific Computing Center, DOE Office of Science User Facility, under Contract No. DE-AC02- BT. Rakhshan was supported by an IVADO Fellowship. O. A. Malik contributed this work while he was affiliated with Lawrence Berkeley National Laboratory where he was supported by the Office of Science of the DOE under G. Rabusseau acknowledges the support of the Canadian Institute for Advanced Research (CIFAR AI chair program) and the Natural Sciences and Engineering Research Council of Canada (Discovery program, RGPIN-2011 ME1IN-2019059944 2019-05949).
05CH1 1231 using NERSC award ASCR-ERCAPO0244
Award Number DE-AC02-05CHHE
This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any ageney thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.
References
Johann A Bengua, Ho N Phien, and Hoang D Tuan. Optimal feature extraction and classification of tensors via matrix product state decomposition. In 2015 IEEE International Congress on Big Data, pages 669-672. IEEE, 2015.
Vivek Bharadwaj. Osman Asif Malik, Riley Murray, Laura Grigori, Aydin Buluc, and James Dem- mel. Fast exact leverage score sampling from khatri-rao products with applications to tensor decomposition. Advances in Neural Information Processing Systems, 36, 2023.
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.
Zhongming Chen, Huilin Jiang, Gaohang Yu, and Liqun Qi. Low-rank tensor train decomposition using tensorsketch. arXiv preprint arXiv:2309. 20t.:2339000093,, 2023.
Dehua Cheng, Richard Peng, Yan Liu, and Ioakeim Perros. Spals: Fast alternating least squares via implicit leverage scores sampling. Advances in neural information processing systems, 29, 2016.
Andrzej Cichocki and Anh-Huy Phan. Fast local algorithms for large scale nonnegative matrix and tensor factorizations. IEICE transactions on fundamentals of electronics, communications and computer sciences, 92(3):708- 721, 2009.
Petros Drineas, Ravi Kannan, and Michael W Mahoney. Fast monte carlo algorithms for matrices i: Approximating matrix multiplication. SIAM Journal on Computing, 36(1):132-157 2006a.
Petros Drineas, Michael W. Mahoney, and Shan Muthukrishnan. Sampling algorithms for & re- gression and applications. In Proceedings of the Seventeenth Annual ACM-SIAM Symposium on Discrete Algorithmm, pages 1127-1136, 2006b.
Petros Drineas, Michael W. Mahoney, and S. Muthukrishnan. Relative-error CUR matrix decomposi- tions. SIAM Journal on Matrix Analysis and Applications, 30(2):844-881, 2008
Petros Drineas, Michael W Mahoncy, Shan Muthukrishnan, and Tamás Sarlós. Faster least squares approximation. Numerische mathematik, 1 117(2):19-2449, 2011.
Glen Evenbly. Gauge fixing, canonical forms, and optimal truncations in tensor networks with closed loops. Physical Review B, 98(8):085155, 2018.

---

Glen Evenbly. A practical guide to the numerical implementation of tensor networks i: Contractions, decompositions and gauge freedom. arXiv preprint arX1v:2202.02238, 2022.
Matthew Fahrbach, Gang Fu, and Mehrdad Ghadiri. Subquadratic kronecker regression with ap- plications to tensor decomposition. Advances Neural Information Processing Systems, 35: 28776-28789, 2022.
Nathan Halko, Per- Per-Gunnar Martinsson, and Joel A Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM review, 53(2): 217-288, 2011.
Christopher Hillar and Lek-Heng Lim. Most tensor problems are np-hard. Journal of the ACM (JACM), 60(6):1-33) 2013.
Frank L Hitchcock. The expression of tensor or polyadic as sum of products. Journal of Mathematics and Physics, 6(1-4):164-18, 1927.
Sebastian Holtz, Thorsten Rohwedder, and Reinhold Schneider. The alternating linear scheme for tensor optimization in the tensor train format. SIAM Journal on Scientific Computing, 34(2): A683-A713, 2012.
Benjamin Huber, Reinhold Schneider, and Sebastian Wolf. A randomized tensor train singular value decomposition. In Compressed Sensing and its Applications: Second International MATHEON Conference 2015, pages 261-290. Springer, 2017.
Pavel Izmailov, Alexander Novikov, and Dmitry Kropotov. Scalable gaussian processes with billions of inducing inputs via tensor train decomposition. In International Conference on Artificial Intelligence and Statistics, pages 726-735. PMLR, 2018.
Tamara G Kolda and Brett W Bader. Tensor decompositions and applications. SIAM review, 51(3): 455-500, 2009.
Jean Kossaifi, Yannis Panagakis, Anima Anandkumar, and Maja Pantic. Tensorly: Tensor learning in python. Journal of Machine Learning Research (JMLR), 20(26), 2019.
Brett W Larsen and Tamara G Kolda. Practical leverage-based based sampling for low-rank tensor decom- position. SIAM Journal on Matrix Analysis and Applications, 43(3):1488 1121148881111 1517,, 2022.
Michael W Mahoney al. Randomized algorithms for matrices and data. Foundations and Trends@ in Machine Learning, 3(2):123-224, 201
Osman Asif Malik. More efficient sampling for tensor decomposition with worst-case guarantees. In International Conference on Machine Learning, pages 14887-14917. PMLR, 2022.
Osman Asif Malik and Stephen Becker. A sampling basedd method for tensor ring decomposition. In International Conference on Machine Learning, pages 7400-7411. PMLR, 2021.
Osman Asif Malik, Vivek Bharadwaj. and Riley Murray. Sampling-bar based decomposition algorithms for arbitrary tensor networks, 2022.
Alexander Novikov, Mikhail Trofimov, and Ivan Oscledets. Exponential machines. arXiv preprint 1.05.06605 2016.
Ivan Oseledets. Tensor-train decomposition. SIAM Journal on Scientific Computing, 33(5): 2295-2317, 2011.
Ninh Pham and Rasmus Pagh. Fast and scalable polynomial kemelss via explicit feature maps. In Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 239-247, 2013.
Nicholas D Sidiropoulos, Lieven De Lathauwer, Xiao Fu, Kejun Huang, Evangelos E Papalexakis, and Christos Faloutsos. Tensor decomposition for signal processing and machine learning. IEEE Transactions signal processing, 65(13):355133522 2017.
12

---

Shaden Smith, Jee W. Choi, Jiajia Li, Richard Vuduc, Jongsoo Park, Xing Liu, and George Karypis.
FROSTT: The Formidable Repository of Open Sparse Tensors and Tools, 2017. URL http //frostt.io/
Edwin Stoudenmire and David J Schwab. Supervised learning with tensor networks. Advances in neural information processing systems, 29, 2016.
Ledyard R Tucker. Some mathematical notes on three-mode factor analysis. Psychometrika, 31()): 279-311, 1966.
David P Woodruff al. Sketching as tool for numerical linear algebra. Foundations and Trends in Theoretical Computer Science, 10(1-2)11117, 2014.
Gaohang Yu, Jinhong Feng. Zhongming Chen, Xiaohao Cai, and Liqun Qi. A randomized block krylov method for tensor train approximation. arXiv preprint vrXiv:2308.01440 2023.
Qibin Zhao, Guoxu Zhou, Shengli Xie, Liqing Zhang, and Andrzej Cichocki. Tensor ring decomposi- tion. arXiv preprint
2016.
H. Zhou, L. Li, and H. Zhu. Tensor regression with applications in neuroimaging data analysis. Journal of the American Statistical Association, 108(502):540-5555 552, 2013.
13

---

7 NeurIPS paper checklist
1. Claims
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Question:
Answer: [Yes]
Justification: The claims supported by theoretical results and experiments.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
[ns
We did not find any particular limitation of our approach that we thought was relevant to discuss the paper.
Justification:
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: The main theorems and lemmas are in section 4 with proofs in appendix.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: In the experiment section we addressed all details about the datasets and the metrics evaluate our method.
5. Open access to data and code
Question: Does the paper provide open access the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
The GitHub repository link is provided in the experiments section.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: All the details all in the experiment section.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer:
[Yes]
Justification: We reported all the results for several runs in plots and tables in the experiment section. However, we did not report the error bars they too small to report.
8. Experiments Compute Resources
Question: For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time execution) needed to reproduce the experiments?
Answer: [Yes]
14

---

Justification: We reported the time each algorithms takes in the experiment section. Details are given in appendix.
9. Code Of Ethies
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics
Answer: [Yes]
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
This paper focuses theoretical aspects. We did not feel relevant to dis- cuss the potential positive and negative societal impacts democratizing access to tensor decomposition methods.
Justification:
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer:[NA]]
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
15

---

A Additional Notation
Definition A.1. Let E Tgmxx and RRXX then the Kronecker product, B Rmpxxnq defined by
011
Oin B
Q21 B
a22 B
a2n B
A
LmmB
0m B
Omn B
Definition A.2. Let A RmxR R and B E R R then the Khatri-Rau product, A B E Runx defined by
A B
01 2 bi
& b2 an & bn
the columns of A. bi...., bR E R are the columns of B and the columns of A B the subset of the Kronecker product. In the corresponding tensor network diagram, the copy tensor captures the fact that the second indices the same.
where ai...., aR E IR
A.1 Details about Orthogonalization of the TT Decomposition
Figure 6 illustrates the single-site TT-ALS method, which begins with a TT decomposition in canonical form initialized by a crude guess. Core Ai of the decomposition sweeps from left-to-right and right-to-left, the algorithm holds all but one core constant and solves the optimal value for the remaining core. After updating each core, by QR decomposition non-orthonormal part is merged the left right (depending the direction of the sweep), step which is called core orthogonalization
non-orthogonal; in
Ri
R1
Ra
Rs
A
step:
F2
R
As
QR
1
R
R
R5
R1
step: 2
8
R
R
QR
in
R
R
Rs
R.
step:
Figure 6: Half-sweep of TT-ALS. In each non-QR step the fully colored core is optimized and in each QR step the non-orthogonal component (depicted by black circle) is absorbed the next core. This procedure repeats until reaching the right side of the decomposition then the same procedure is repeated from right until reaching to the left side (not demonstrated in this figure.)
16

---

B Proofs
B.1 Proof of Lemma 4.1
Noting that Ag, 11...11...5
s3. is row vector, write
p(š 8k
83k 82k)
Σ
p(81 81 A A...A.j -.. =z)
E
1
A 81 11....88 89.:1 Asy(s1....
R2
Σ
-r
A< 81 1
Agi 11..
R
81
A li sg ... A (1.1d
Ai (1s1 ... A, s3.
R
$1
EA A 11...1 .... ... (i
Tr
A₂ 1.s.s.1
Rg
R Σ Tr [[(](188,)))
A :2 T A2 [; 82, A, A1
In the expressions above, the summation over each variable 85, 1 1 S k. is taken over the range [T,]- The first step follows by marginalizing over random variables $,,..., - The second step follows from Equation (4). The third step rewrites an inner product two vectors as the trace of their outer product. The fourth step follows from the definition of Agj. The fifth step follows from the linearity of the trace by moving the summation over $1 into the product expression. The last step follows from the definition of the left-orthonormality property on A that is, Σ., Ai [t. s1, -A [: s1,:] = ArT Af I. By successively moving summation operators into the product expression to repeat the last step (exploiting the lehhoorthoyyymality of each core in the process), we find
p(šx
sp.
8>* s>k)
Tr
A, 8j A [ [1, sk.:] AA [s, 8k,:] ..... Ai 83,8
R
(8)
1
Tr [MTx A J sl As 1., sk Hs
R2
where the last line follows from the definition of H
B.2 Proof of Lemma 4.2
We write
1
p(g 3h 18
syk)
Tr
H
Au s]:""
As [1, sk. HsK
Rg
Rg
1
e
HJRh A sk As [1, sk, H er)
R
r=1
B2
(9)
[1Tx As [t, sk:
As surt hz8)
R2
r=1
R
1
p(fs to bs t =r))
R
The first step follows from Lemma 4.1 The second step follows from the definition of the trace. The third step follows from the definitions of hy* and H The fourth step follows from the definition of the variables ..,, f,. Now observe that p(f r) 1/R, for 1 < Rí, write
17

---

R,
p(8k = 8k
85n = 85k)
=
Σ
p(1x = Lx 12a = r)p(f = r)
(10)
r=1
=p(Ox 4. 1isx= = bzx);
which completes the proof.
B.3 Efficient Sampling Data Structure
Lemma 4.3 first appeared as Lemma 3.2 in the original work by BBaradwajj et al. 2023]. state a condensed form of the original claim below:
Lemma B.1 Bharadwaj et al. 2023], Original), Given U E RXR Y E RRXR with Y p.s.d., there exists a data structure parameterized by positive integer that requires OMRR)) time to construct and additional space space O [M] F/ After construction, the data structure draw sample from the distribution defined elementwise by
UuYU.Y [s] C-- [s, hhh) u [s,:
When Y is a rank-1 matrix, the runtime drops to in time O(R log log[M/F M/P + FRR). O log] M/ + FR).
In the statement above, C is appropriate normalization constant. To prove our adapted lemma, take Y = [1], a matrix of all ones that is rank-1, and set F = R. Then
Uu,U.Y C- [s, (hhT)
[s.]] = C- hh)
This is the target probability distribution of Lemma 4.3 and the runtime to draw each sample is O(R² log(M) RR) R²) O log( M/ R.). The choice F= R also induces space usage O (M RR. linear in the size of the input. Our modified claim follows.
B.4 Proof of Theorem
We provide short end-to-end proof that shows that Algorithms ] [add] correctly draw samples from Agı (the matricization of the left-orthogonal core chain) according to the distribution its squared row norms while meeting the runtime and space guarantees of Theorem 1.1
Construction Complexity: The cost of Algorithm follows from 4.3 with M IR the row count of A for < k S Using this lemma, construction of each sampling data structure Za requires time O - R2). The space required by sampler Zμ R summing over all indices l gives the construction claim in Theorem
The complexity to draw samples in Algorithm 2 is dominated by calls to the RowSample procedure, which as discussed in Section 4 O log R Rg)) Summing the complexity over indices 1 < k < j yields the cost claimed by Theorem 1.1 to draw single sample. The complexity of calling the RowSample procedure repeatedly dominates the complexity to
Sampling Complexity:
update the history vector h over all loop iterations, which is O | R
for each sample.
Our task is to show that Algorithm 2 each sample td 1 1 d < J. is multi-index that follows the squared row norm distribution on the rows of Agj. To do this, we rely on lemmas proven earlier. For each sample, the variable is uniform random draw from [Ry], and is initialized the corresponding basis vector. By Equation (6) and Lemma 4.3] Line 5 from Algorithm 2 draws each index l correctly according to the probability distribution specified by Equation (5). The history vector updated by Line 6 of the algorithm so that subsequent draws past iteration & of the loop are also drawn correctly according to Equation (5). Lemma 4.2 (relying on Lemma 4.1] shows that the multi-index f ... f, drawn according to Equation (5) follows the same distribution as $1 $.. which was defined follow the squared norm distribution on the rows of Agj- This completes the proof.
Correctness:
D
18

---

B.5 Proof of Corollary 4.4
Since A#j RT R and X() Rl TexI, we draw ((1F/Ed) samples to achieve the error bound (1 + e) with probability (1 - 8) for each least squares solve in the down-sampled problem 3.5). By Theorem I.1. the complexity of drawing J samples with our data structure is
Σ kサj
log 1&11? -
Σ Rt ((50) log I Kgン
where we suppose that R1 R₂ .... RN-1 and I
IN. The cost of sampling corresponding subset of ()) is (JJ1) = -((((())))) Solving the downsampled least squares problem also costs O ( = Ô (1, R² ((0)) Summing them all together for 1 S S N gives
N
Ô 1/58
Σ
Σ kgj
R log Is
+ IRI1)
N
Rt /ed
Σ
(N 1) log 1 +
j=1
=
R /e0
Σ
N log 1, + +
where we wrote the last equation considering the fact that N dominates (N 1).
Details about Datasets & Experiments
C.1
Datasets
For the real dense datasets experiment, we truncated and reshaped the original data tensors in to the fourth order tensors as follows.
Pavia University dataset: The original has dimensions (610,340, 103). We truncate it to (600, 320, 100), permute the modes to dimensions (100, 320, 600) tensor and reshape it into tensor of dimensions (100, 320, 24, 25). It is available
attp://lesun.weebe www.wwwwlloooooo.. wweebly..com/yyperspettr... datasaltdddattttttttaaaa
Tabby Cat dataset is permuted to (286, 720, 1280) and reshaped to tensor of size (286, 720, 40, 32). The video is color and converted to grayscale by averaging the three color channels. It is available at
The MNIST dataset was reshaped into a tensor of size (200, 600, 10) and is available
at https: http:///agg..a..aaaaaaaaaaa//aaa
The Washington DC Mall dataset was truncated to dimensions (1280, 306, 190) before
reshaping into a tensor size (1280, 306, 10, 19). It is available at
https://engineett purdue edu /aau Embiehiiiii MaaisiiehhMMMlleeee hwwpeehhhhhhhthpaeeeeee /ttpa/aat...eeepppppppcaaaaaaaaa httml.
The sparse tensors Uber, Enron, and NELL-2 were downloaded from the FROSTT collection (Smith et al., 2017]. The dimensions of these tensors were unchanged from the versions available online. Consistent with established practice [Larsen and Kolda 2022], we computed the logarithm of the tensor values in the Enron and NELL-2 datasets before performing our experiments.
C.2 Computing Resources
The dense data experiments were conducted on MILA cluster nodes with 4 CPUs and 16GB of RAM each. Sparse tensor decomposition experiments were conducted on NERSC Perlmutter nodes with CPUs and 512 GB of RAM each.
19


Reviews:
[{"review_id": "aQIL0Vfodl", "reviewer": "Reviewer_yjyo", "summary": "This paper considers the problem computing the tensor train TT decomposition of large tensors and proposes a novel randomized approach for efficiently solving the problem. In particular, the Alternating Least Squares (TT-ALS) algorithm is considered and exact leverage scores sampling approach is proposed to accelerate the approach. A data structure based approach is devised to efficiently compute the leverage scores of the appropriate matricization of the tensors and sample rows from them. Theoretical analysis show the sampling complexity. Numerical experiment results are presented which to illustrate the performance of the proposed algorithm.", "questions": "The paper presents an interesting approach to compute the TT decomposition of large tensors. \nHowever, I have the following comments that can further improve the paper:\n\n1. r-TT-SVD seems better - \n\nBased on numerical results presented, it appears  the randomized SVD based approach  (r-TT-SVD) takes significantly lower time, yet achieves comparable performance (fit) as other methods. Are TT-ALS and rTT-ALS popular and relevant for TT decomposition. How do these methods differ in the quality of the decompositions computed?\n\nIn some cases, rTT-SVD might have a slightly lower fit. But, since the runtime is so low, perhaps a larger sketch size  would yield similar fit at lower cost. Details about the sketch sizes used are missing. \n\n\n2. Presentation:\nCertain aspects of the presentation can be improved, and there are few minor typos.\n\ni. Are TT decompositions popular in tensor applications? Overall, they seem to have higher computational cost and are not very interpretable. Few alternate tensor decompositions achieve better compression with lower cost. A discussion on the motivation for the use of TT decomposition will benefit the paper (readers can better appreciate the results).\n\nii. Introduction has tensor jargons which might not be known to general AI audience, such as 3D core chains, left-matricization, contraction of tensor, etc. It is better to avoid these or define them before use. \n\niii. The computational  cost of solving eq (1) can be added.\n\niv. Theorem 1.1, point 1 has a $j$ missing in the runtime cost. \n\nv. Not sure what is the subscript (2) in Corollary 4.4.\n\nMany of the results in the paper is based on another Arxiv paper (which is not peer-reviewed). So, correctness of these results is not established.", "limitations": "See above comments\n\n\n------\nPost Rebuttal:\n\nI have read other reviews and authors' responses to all reviews. The responses have adequately addressed my concerns. The work has many merits. I have raised my score.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The strengths of the paper  are:\n1. A novel randomized approach is presented for computing the TT decomposition of large tensors.\n2. New data structure and sampling approach are proposed. Theoretical results are presented.\n3. Numerical results show that the proposed method is efficient.", "weaknesses": "The weaknesses of the paper are:\n1. The randomized SVD approach seems to be more efficient in terms cost compared the the prosed method.\n2. Certain aspects of the paper can be improved (see below).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their comments and feedback, please find responses below:\n\n**Weaknesses:**\n> The randomized SVD approach seems to be more efficient in terms cost compared the the prosed method.\n\nRandomized SVD cannot scale to large tensors. Indeed, to decompose a tensor using randomized SVD at each step we need to produce a Gaussian tensor of size $s_j\\times I^{(N-1)}$, where N is the order of the tensor, $I$ is the physical dimension of each mode and $s_j$ is the size of the sketch. This approach works well when $N$ is small. But randomized svd cannot handle very high-order tensors. In the experiment part, we only use tt-svd and randomized tt-svd for the dense tensors with small $N$, but these two algorithms cannot scale to the case of large sparse tensors in the second part of the experiments.\n\n**Questions:**\n> The paper presents an interesting approach to compute the TT decomposition of large tensors. However, I have the following comments that can further improve the paper:\nr-TT-SVD seems better. Based on numerical results presented, it appears the randomized SVD based approach (r-TT-SVD) takes significantly lower time, yet achieves comparable performance (fit) as other methods. Are TT-ALS and rTT-ALS popular and relevant for TT decomposition. \n\nyes : TT-ALS is a very popular algorithm, notably in the quantum physics community, where it is closely related to the density renormalization group (DMRG) algorithm. \n> How do these methods differ in the quality of the decompositions computed?\n\nRegarding the quality of the decomposition, TT-ALS usually finds better solutions than TT-SVD (TT-SVD can actually be used as an initialization for TT-ALS). Note however that the focus of our paper is on improving the TT-ALS algorithm using randomization. Indeed, TT-ALS is a popular algorithm, but, as mentioned previously, the cost of solving the least-squares problems in ALS is exponential in the order of the tensor which is not efficient for decomposing very large tensors. The main goal of our paper is to make the ALS efficient with randomization and doable for high-order tensors. Therefore, the improvement of a TT-ALS algorithm is the main focus of our paper.\n\n> In some cases, rTT-SVD might have a slightly lower fit. But, since the runtime is so low, perhaps a larger sketch size would yield similar fit at lower cost. \n\nThat’s correct, by increasing the sketch size we can have a better quality of sketch for rTT-SVD. However, as mentioned earlier, the SVD-based approach is not suitable for decomposing high-order tensors. Even if by increasing the sketch size rTT-SVD still suffers from the curse of dimensionality when N is very large. Therefore, the quality of the sketch of rTT-SVD cannot be compared to rTT-ALS in the case of sparse high-order tensors (section 5.2 second part of the experiments.)\n\n> Details about the sketch sizes used are missing.\n\nIn Section 5.1, the sketch sizes are set to $J=5000$ and $J=2000$ for the synthetic and real data, respectively (we mention it in the text but we will add this information in the caption of Figure 3 and Table 1 as well). \n\n> Presentation: Certain aspects of the presentation can be improved, and there are few minor typos.\ni. Are TT decompositions popular in tensor applications? Overall, they seem to have higher computational cost and are not very interpretable. Few alternate tensor decompositions achieve better compression with lower cost. A discussion on the motivation for the use of TT decomposition will benefit the paper (readers can better appreciate the results).\n\nTT decomposition enjoys stable and numerical stabilities while finding rank-r CP decomposition is NP-hard and the number of  parameters of a Tucker decomposition grows exponentially with the order of the tensor. In contrast, the number of parameters is linear in the order of the tensor for both TT and TR decompositions, but TR decomposition is known to suffer numerical stability issues.  For these reasons, TT is more popular in tensors and quantum physics communities (where it is known as Matrix Product States / MPS). To address your concern about motivating the TT format, we will include a more comprehensive motivation in the introduction.\n\n>  ii. Introduction has tensor jargons which might not be known to general AI audience, such as 3D core chains, left-matricization, contraction of tensor, etc. It is better to avoid these or define them before use.\n\nIf we understand correctly, your concern is mostly about the main theorem in the introduction which may cause difficulties for other audiences. We will clarify this by clearly stating before the theorem that all relevant definitions are given in Section 3.1.\n\n> iii. The computational cost of solving eq (1) can be added.\n\nSolving Eq (1) exactly (without randomization) has a computational cost of $O(I^N)$. It is mentioned in the introduction but we will emphasize it in the revision in section 3.2.\n\n> iv. Theorem 1.1, point 1 has a 𝑗 missing in the runtime cost.\n\nCorrect. We will add j in the final version. Thank you for catching this typo. \n\n> v. Not sure what is the subscript (2) in Corollary 4.4.\n\nThis denotes the second mode unfolding of the tensor which is defined in Definition 3.1 in the paper.\n\n> Many of the results in the paper is based on another Arxiv paper (which is not peer-reviewed). So, correctness of these results is not established.\n\nThe paper “Fast exact leverage score sampling from Khatri-rao products with applications to tensor decomposition.” by Bharadwaj et el. has been published in Neurips 2023, we will correct the citation in the camera-ready version.", "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": "", "extra_scores": {}}, {"review_id": "oLgOnSFeto", "reviewer": "Reviewer_BTzr", "summary": "This paper presents an efficient algorithm to use leverage-score sampling to solve least squares problems arising as subproblems in a larger alternating least squares (ALS) algorithm for building an approximate Tensor Train (TT) decomposition. The paper reports empirical evaluation of the proposed algorithm, showing runtime improvement compare to the baseline un-sampled problems.", "questions": "- Lines 96-97: Sketching is a good way to circumvent the need to compute leverage scores, and TensorSketch allows us to sketch matrices with Kronecker structure. And indeed, the authors point out that [Chen et al. 2023] do this. How do their results compare to the present paper? Why is what you are proposing a better idea than using TensorSketch?\n- Lines 64-66: What is the significance of this observation regarding the complexity of the algorithm?\n- Lines 125 - 126: One the canonical form has been computed, can we solve the ALS subproblem efficiently?\n- Did you compare empirically your method to the method of [Chen et al. 2023].\n- Is your algorithm computing the canonical form anew in each ALS iteration, or are you somehow updating it?\n- Does your algorithm maintain/computes the leverage scores of each row (the number of such scores is exponential), or do you present only an efficient way to sample from it without forming it?", "limitations": "Nothing to add.", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- Empirical evaluation show a runtime improvement with respect to unsketched TT-ALS, without reducing fit.\n- The paper presents an efficient algorithm for leverage score sampling of ALS subproblems in TT-ALS.", "weaknesses": "Major comments (affecting recommendation):\n- The proposed method is compared empirically only to a small number of baselines. Seems that the most likely competatitor (Chen et al. 2023) is not compared against.\n- Runtime improvement in the experiments is modest. This is a sampling based method, so if sampling is very efficient (due to efficient leverage score approximation) I would expect a big improvement in runtime. This might not be the case if either more iterations are now needed for convergence or computing the sampling probabilities is very expensive. Unclear from then experiments which is the case.\n- The paper does not present an analysis end-to-end. It only analyzes how to subsample the ALS subproblems, and thus accelerate them. Not analyzed is how the solution of the subsampled system in lieu of the full subproblem affect ALS convergence.\n- Many details required for implementation of the algorithm are omitted. Many points regarding the algorithm are not clear. \n\nMinor comments (do not affect the recommendation):\n- Lines 125 - 126: add citation.\n- Figure 3 left: is this fit in the y-axis or misfit?\n- Line 108-109: I think you mean *omission* of multiplicative terms.\n- Lines 58-60: shouldn't runtime be O(jIR^3)?\n- Line 54: The term \"left matricization\" is used without defining it yet.\n- Line 55: I think it should be R_j\n- Eq (2): I think on the right it should be A[i,:]^T. Similar issue with Eq (3).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their comments and respond to the feedback below:\n\n**Weaknesses:**\n> [W1][Comparison against (Chen et al. 2023)]\n\nThe TensorSketch approach by Chen et. al. requires an **exponential** sample count in the tensor dimension q to achieve the $(\\epsilon, \\delta)$-guarantee on the accuracy of each least-squares solve. See Theorem 4.4 of their work and note the 3^q term in the sample complexity, where (in their notation) q is the tensor dimension. By contrast, our algorithm requires a sample count that has **no dependence** on the tensor dimension q; the sample count depends only on the column count of the design matrix and the parameters $\\epsilon$nd $\\delta$. The time to draw these samples in our method is linearly proportional to q, the tensoring dimension, a modest cost. The efficiency of our method stems from the fact that we sample from the exact distribution of squared row norms of the design matrix. By contrast, TensorSketch applied without tree structure modifications, such as those proposed by Ahle et al., suffers from an exponential output sketch size. We will clarify this point at the end of Section 2 and in the Introduction.\n\n> [W2][Runtime improvement]\n\nThe runtime improvements are most significant for large sparse tensors. Figure 4 shows the accuracy (y-axis, higher is better) vs. ALS iteration time for our method vs. non-randomized ALS. The speedup per iteration can be as high as 26x for lower ranks. For the NELL-2 tensor, the plot shows that accuracy within three significant figures of non-randomized ALS was achieved roughly 3-4x faster than an optimized non-randomized ALS baseline. We will note this at the beginning of Section 5 (Experiments) and will highlight these facts further in Section 5.2.\n\n> [W3][End-to-end analysis]\n\nCorrect. Our analysis provides guarantees on the accuracy of individual least-squares problems, in line by prior works by Cheng et. al, Larsen and Kolda, and Malik. Under additional assumptions on the tensor, convergence guarantees can be derived for sketched ALS: see, for example, https://proceedings.mlr.press/v119/gittens20a/gittens20a.pdf. The global guarantees derived in this paper by Aggour, Gittens, and Yener depend on sketching guarantees established for each linear least-squares problem, which we provide in our work.\n\n> [W4][Implementation details]\n\nThe revised version of our draft will include (in section 4.1) detailed descriptions of the BuildSampler and RowSample functions, the data structure described in the appendix to draw samples from each orthonormal core flattening. Our Git repository, which is public at https://anonymous.4open.science/r/fast_tensor_leverage-EB01, includes a simple, slow reference implementation written in Python for this data structure, providing sufficient detail to understand and fully replicate our method. \n\n> Minor comments\n\nThanks for catching these typos, we will correct and address them all. For Figure 3, the y-axis is the fit. \n\n**Questions:**\n>[Q1] [Our approach vs TensorSketch]\n\nWe believe we addressed this point in the comment above on the exponential complexity of TensorSketch without appropriate modifications. Our method offers sub-exponential sample complexity and worst-case runtime complexity to achieve guarantees on the solution to each least-squares problem. The approach by Chen. et. al. requires a worst-case sketch size exponential in the tensor dimension.\n\n> [Q2][Complexity of the algorithm] \n\nThe significance here is that data structure construction (as well as subsequent updates to the data structure) does not increase the asymptotic complexity of ALS tensor train decomposition.\n\n> [Q3][proceeding only with the canonical form]\n\nEven after computing the canonical form (which makes the design matrix orthonormal for each linear least squares problem orthonormal), we must multiply the matricized tensor against the chain of TT cores placed in canonical form. **This is the major computational bottleneck** - without sketching, the runtime cost for this matrix multiplication scales as  $O(nnz(T) N R^2)$, where $nnz(T)$ is the number of nonzeros in the tensor, $N$ is the order of the tensor, and $R$ is the tensor train rank. The tensor may have hundreds of millions of nonzero values. Sketching allows us to select only a subset of rows from the design matrix and the corresponding subset of rows from the matricized tensor, reducing the cost to $O(N R^2)$ * (number of selected nonzeros extracted by the sketch). We can update our draft to clarify this point in Section 1.1.\n\n> [Q4][Emprical comparison with Chen et al. 2023]\n\nWe did not compare empirically with Chen et. al results. To the best of our knowledge, there is no public code available to replicate their results and compare with them. Also, as explained previously, from the theoretical point of view, the method proposed by Chen et. al cannot scale to the size of the tensors that we consider in the large sparse tensors experiments. \n\n> [Q5][Updating canonical form in each iteration]\n\nWe update the canonical form without recomputing it entirely at each substep of ALS. After solving for the optimal value for each core, we compute the QR decomposition of its left or right flattening and replace the newly-computed core with the appropriately-reshaped Q factor. In this way, only one core needs to be updated after each least-squares solve to maintain the canonical form.\n\n> [Q6][Computing the leverage scores]\n\nWe do not compute the leverage scores of all rows, since, as you pointed out, there are exponentially-many rows in the tensor dimension. To design a computationally-efficient sampling algorithm, we build an efficient data structure to **sample** from this distribution without materializing it.\n\nWe hope we have addressed all your concerns and answered your questions (in which case we kindly ask you to consider increasing your score). We are happy to clarify any additional point during the discussion period.", "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": "", "extra_scores": {}}, {"review_id": "rDFog0SrEu", "reviewer": "Reviewer_oeCx", "summary": "The authors proposed a leverage score sampling-based TT-ALS method to reduce the computational complexity of the traditional TT-ALS. Experimental results verify the performance of the proposed method.", "questions": "1. My main concern about this work is its relation to the leverage score sampling-based TR-ALS [Malik and Becker, 2021]. As TT can be seen as a special case of TR, what is the main contribution of the current work compared with leverage score sampling-based TR-ALS? Please clarify that it is necessary to develop a new method beyond [Malik and Becker, 2021] (i.e. has a better theoretical guarantee or has better performance than leverage score sampling-based TR-ALS).\n\n2. For real data experiments, how about the performance compared with other tensor decomposition structures, such as TR? As you can set the first and last TR-rank as 1 thus it reduces to a TT-like structure.", "limitations": "N/A", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper is well written with good theoretical analysis and desired experimental performance. The method using leveraging score sampling is technically sound for TT decomposition.", "weaknesses": "The contribution and novelty are not clearly stated compared with [Malik and Becker, 2021]. In [Malik and Becker, 2021], the leverage sampling was applied to TR. As TT can be treated as a special case of TR, the author does not clearly state that it is necessary to develop a new method for TT, or that there are some new theories or findings that are distinct from the previous TR structure.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their comments and feedback, please find responses below:\n\n**Weaknesses:**\n> The contribution and novelty are not clearly stated compared with [Malik and Becker, 2021]. In [Malik and Becker, 2021], the leverage sampling was applied to TR. As TT can be treated as a special case of TR, the author does not clearly state that it is necessary to develop a new method for TT, or that there are some new theories or findings that are distinct from the previous TR structure.\n\nIn [Malik et al, 2021], the leverage scores were approximated using a product of simpler distributions. They sample according to the leverage scores of each core. While in our paper, we proposed a novel data structure that can compute the exact leverage scores. Moreover, the time complexity to solve one least-square of ALS with the method of Malik et al.  Malik et al approach is $NIR^4$ +#iter. $NIR^{(2N+2)}/(\\epsilon\\delta)$, which is still exponential in the order of the tensor $N$, while our result does not have any exponential dependency to the order of the tensor. We will add and emphasize in the camera-ready version that our approach differs from theirs in two main aspects: \n- We did not use approximation for finding the leverage scores, instead a novel data structure is proposed for finding **exact** leverage scores.\n- In Malik et al, the runtime for the least-square solve still has an exponential dependency to the order of a tensor while our approach is free of this exponential dependency. \n\n**Questions**\n> My main concern about this work is its relation to the leverage score sampling-based TR-ALS [Malik and Becker, 2021]. As TT can be seen as a special case of TR, what is the main contribution of the current work compared with leverage score sampling-based TR-ALS? Please clarify that it is necessary to develop a new method beyond [Malik and Becker, 2021] (i.e. has a better theoretical guarantee or has better performance than leverage score sampling-based TR-ALS).\n\nIndeed, TT is a special case of the TR decomposition. However, in [Malik and Becker, 2021], leverage score approximation is used for the TR-ALS with time complexities #iter · $NIR^{2N+2}/(\\epsilon\\delta)$ (further improved to and #iter · $N^3R^8 (R + I/\\epsilon)/\\delta $ in [Malik, 2022]), respectively. In contrast, our method has a lower time complexity and provides a better theoretical guarantee (please see corollary 4.4 in our paper). Therefore, our method goes beyond the works of Malik et al. As mentioned in their works, computing leverage scores requires computing the left singular values which has the same cost as solving the original least squares problem. For this reason, they approximate the leverage scores. Instead, we propose a novel data structure to sample from the **exact** leverage scores at a lower cost than their proposed approaches.\n\n> For real data experiments, how about the performance compared with other tensor decomposition structures, such as TR? As you can set the first and last TR-rank as 1 thus it reduces to a TT-like structure.\n\nIt is true that TT is a special case of TR, but note that the main problem we address in our paper is how to scale TT-ALS to very large tensors using randomized techniques. The purpose of our experiments is to demonstrate how we efficiently achieve this goal with exact leverage score sampling. Extending our approach to TR-ALS could be interesting, however it would be very challenging due to the lack of canonical forms for TR decomposition. \n\nWe hope we have addressed all your concerns and answered your questions  (in which case we kindly ask you to consider increasing your score). We are happy to clarify any additional point during the discussion period.", "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": "", "extra_scores": {}}, {"review_id": "loJTyRKbex", "reviewer": "Reviewer_8ErW", "summary": "This paper gives a better randomized alternating least squares algorithm for computing tensor factorizations. It's based on an exact characterization of leverage scores of the matrization of tensors via a suitable intermediate orthonormal representation. This is justified rigorously, and significant empirical gains were demonstrated.", "questions": "Would it be possible to check how the orthogonality conditions are maintained in intermediate steps under the inexact arithmetic caused by round off errors? Aka. are the conditions for the characterizations of leverage scores preserved exactly? (this has been answered)", "limitations": "yes", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The approach taken is natural, and the bounds obtained are quite powerful.\n\nThe experiments considered both synthetic and real data sets, and demonstrated clear gains in the parameter regimes considered.", "weaknesses": "I'm a bit concerned about the setting of the experiments, which seem to be 3-dimensional dense tensors. My understanding is that a lot of the more complicated tensor instances are sparse and in higher dimensions. However, I'm not sure whether those have low rank representations.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their comments and feedback, please find responses below:\n\n**Weaknesses:**\n> I'm a bit concerned about the setting of the experiments, which seem to be 3-dimensional dense tensors. My understanding is that a lot of the more complicated tensor instances are sparse and in higher dimensions. However, I'm not sure whether those have low rank representations.\n\nCorrect. We added dense synthetic and real data experiments to show that our proposed approach has a better time complexity than TT-ALS and TT-SVD and works slightly better than rTT-SVD in terms of fit. The main point is to demonstrate that our new approach works as well as other types of TT decompositions for dense data. As SVD-based methods cannot handle high-dimensional sparse tensors. For the more complicated tensors, we added a sparse tensor data experiment as you can see in the second part, the setting is different and more complicated as SVD-based TT decomposition cannot be used and we just compared our approach with the classical TT-ALS method.\n\n**Questions:**\n> Would it be possible to check how the orthogonality conditions are maintained in intermediate steps under the inexact arithmetics caused by round off errors? Aka. are the conditions for the characterizations of leverage scores preserved exactly?\n\nWe did not observe any numerical instabilities due to a round of errors during the orthogonalization step (which is performed using the stable QR decomposition routine provided by numpy.linalg) and we believe the characterization of leverage scores is preserved exactly (up to machine precision).", "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": "", "extra_scores": {}}]

Author rebuttal:
We thank all the reviewers for their comments and feedback, which we all answered in the individual rebuttals below. 
We summarize some of our answers to the main points raised by the reviewers (see individual rebuttals for details).

- [SVD-based vs ALS-based approach]: Overall TT-ALS is very popular in the quantum physics community and usually finds better solutions than TT-SVD. Randomized TT-SVD cannot scale to high-order tensors. Each step of the randomized TT-SVD decomposition requires generating a random Gaussian matrix (classical random projection) that only can handle small-order tensors. 

- [Comparison with Chen et al, 2023 paper]: The TensorSketch approach proposed by Chen et al, 2023~(https://arxiv.org/pdf/2309.08093) requires **exponential** sketch size in the tensor dimension while our algorithm requires the sketch size that has no dependence on the tensor dimension; the sketch size depends only on the column count of the design matrix and the parameters $\epsilon$ and $\delta$.  

- [Comparison with Malik and Becker, 2021 paper]: The approach proposed by Malik et al, 2021~(https://proceedings.mlr.press/v139/malik21b/malik21b.pdf) requires **approximating** the leverage scores. However, in our paper, we proposed a novel data structure that computes the **exact** leverage scores. Moreover, in Malik et al, the runtime for the least-square solve has an **exponential dependency** to the order of a tensor while our approach is **free** of exponential dependency to the order of a tensor. 

- [Dense and sparse tensors experiments]: For the experiment part, we added two separate sections, one for the dense and one for the sparse tensors. Our purpose for the dense tensor section is to show that our approach has a better time complexity than TT-ALS and TT-SVD and matches rTT-SVD in terms of fit. However, SVD-based decompositions cannot handle high-order (sparse) tensors. 
For the sparse tensors section, we compared our proposed approach with the classical TT-ALS.

We hope we have addressed all your concerns and answered your questions and we are happy to clarify any additional points during the discussion period.

Meta-review: After the rebuttal period and inter-reviewer discussion, the final ratings for this paper are {7, 5, 6, 7} since Reviewer BTzr raised their score from 4 to 5, allowing us to reach a unanimous accept decision.

I read the paper myself and am also fine with it being accepted. It's not as crisp as it should be (theory and experiments), but a paper on leverage score-based sampling for Tensor Train (TT) decompositions is missing from "faster tensor decomposition via leverage scores" line of work (CP decomposition in NeurIPS 2016, Tucker decomposition in NeurIPS 2022).
Decision: Accept (poster)