PaperlayerView exact revision on arXiv

Source-linked AI summary

Practical Leverage-Based Sampling for Low-Rank Tensor Decomposition

Brett W. Larsen, Tamara G. Kolda

arXiv:2006.16438v3math.NAcs.NA

TL;DR

Sparse CP decomposition requires solving least-squares subproblems that can remain costly on very large tensors. The paper sketches these systems using Khatri–Rao structure and leverage-score bounds, with row-combining and hybrid sampling strategies. On an Amazon tensor with 1.7 billion nonzeros, CP-ARLS-LEV was faster and achieved a better solution; the paper also reports order-of-magnitude speedups on large sparse tensors.

  • Problem

    Sparse CP-ALS least-squares solves can be expensive because their cost scales with tensor nonzeros, which can number in the billions.

    Source: §1 · p. 1

  • Method

    CP-ARLS-LEV sketches CP-ALS least-squares systems by sampling rows using leverage-score bounds and practical strategies for repeated or high-probability rows.

    Source: §1 · p. 1, §8 · pp. 22–23

  • Results

    On an Amazon tensor with 1.7 billion nonzeros, CP-ARLS-LEV was faster and achieved a better solution, and the method yielded order-of-magnitude speedups on large-scale sparse tensors.

    Source: §8 · pp. 22–23

  • Takeaways & Limitations

    Combining repeated rows and including high-probability rows deterministically are practical tools for leverage-score matrix sketching beyond tensor decomposition.

    Source: §8 · p. 23

    The sample-count and deterministic thresholds were chosen empirically, and low-variance fit estimates on large tensors required extremely many samples.

    Source: §8 · p. 23

Abstract

from arXiv · show

The low-rank canonical polyadic tensor decomposition is useful in data analysis and can be computed by solving a sequence of overdetermined least squares subproblems. Motivated by consideration of sparse tensors, we propose sketching each subproblem using leverage scores to select a subset of the rows, with probabilistic guarantees on the solution accuracy. We randomly sample rows proportional to leverage score upper bounds that can be efficiently computed using the special Khatri-Rao subproblem structure inherent in tensor decomposition. Crucially, for a $(d+1)$-way tensor, the number of rows in the sketched system is $O(r^d/\epsilon)$ for a decomposition of rank $r$ and $\epsilon$-accuracy in the least squares solve, independent of both the size and the number of nonzeros in the tensor. Along the way, we provide a practical solution to the generic matrix sketching problem of sampling overabundance for high-leverage-score rows, proposing to include such rows deterministically and combine repeated samples in the sketched system; we conjecture that this can lead to improved theoretical bounds. Numerical results on real-world large-scale tensors show the method is significantly faster than deterministic methods at nearly the same level of accuracy.

1. Introduction.

CP decomposition solves a sequence of least-squares problems, but sparse tensors can make even the standard solves expensive. This paper sketches those systems using Khatri–Rao structure and leverage-score bounds, with a practical algorithm and sampling strategies for high-probability rows.

1. Introduction.

  • Sparse CP-ALS solves can be costly even when exploiting sparsity, since their cost is proportional to the number of nonzeros, which can reach billions.The paper targets a different least-squares step from Cheng et al. and presents methods for repeated and high-probability sampled rows.Source: §1 · p. 1

1.1. CP least squares problem.

  • Each CP-ALS iteration solves a sequence of overdetermined least-squares problems, with a Khatri–Rao product design matrix in each subproblem.For sparse tensors, the normal-equations approach can reduce the solve cost to O(nnz(X) r).Source: §1.1 · pp. 1–2

1.2. Related work.

  • Unlike KFJLT sketching [6], which requires FFT preprocessing that destroys sparsity, leverage-score sampling is presented as an alternative for sparse tensors.For the Uber tensor, the sparse representation needs 0.14 GB versus 68 GB for the KFJLT-preprocessed dense version.Source: §1.2 · pp. 2–3

1.3. Our contributions.

  • For suitably small ε, CP-ARLS-LEV uses O(r^d/ε) sampled rows for an ε-accurate solution with high probability, independent of system size and tensor nonzeros.On the large Reddit tensor with 4.7 billion nonzeros, it achieved up to 16× speedup over CP-ALS, reducing compute time from about 4 days to 6 hours.Source: §1.3 · p. 4
  • The paper details CP-ARLS-LEV, proves its sample bound, proposes combining repeated rows and deterministic inclusion of high-probability rows, and reports numerical experiments.For comparison with Cheng et al. [8], the paper sketches the full least-squares problem and uses factor-matrix leverage scores rather than sampling based on rows of the unfolded tensor.Source: §1.2 · p. 3, §1.3 · p. 4

2. Background on least squares problems in CP-ALS and KRPs.

A CP-ALS iteration updates one factor matrix at a time by solving a least-squares subproblem. Its design matrix is a Khatri–Rao product whose rows correspond to tuples of factor-matrix rows.

2. Background on least squares problems in CP-ALS and KRPs.

  • Each of the d + 1 CP-ALS subproblems has the same format; the exposition focuses on the update for the (d + 1)st factor matrix.For the first factor, the unfolding and Khatri–Rao product change to match the other fixed factors.Source: §2 · p. 4
  • The discussion assumes all factor matrices and Z have rank r.The factor matrices are taken to have dimensions n_k × r.Source: §2 · p. 4
  • The design matrix Z is the Khatri–Rao product of the fixed factor matrices, and its rows correspond bijectively to tuples of rows from those matrices.The tuple is the multi-index and the corresponding row number in Z is the linear index.Source: §2 · p. 4

3. Background on sketching for least squares problems.

Row sketching aims to reduce least-squares costs by sampling and weighting rows so the sketched system preserves the original solution quality. Leverage scores guide sampling, while bounds on them make the approach practical for Khatri–Rao systems.

3. Background on sketching for least squares problems.

  • Sketching targets the least-squares cost’s dependence on N, reducing it to O(sr^2 + snr) before further exploiting sparsity and avoiding explicit Khatri–Rao formation.The sampling matrix is restricted to one nonzero per row so that it selects and weights rows.Source: §3 · pp. 4–5

3.1. Weighted sampling.

  • Weighted row sampling makes the sketched squared norm an unbiased estimator of the original squared norm when sampled rows have positive probability wherever the vector is nonzero.The probability distribution affects estimate quality through the leverage scores of Z.Source: §3.1 · p. 5

3.2. Leverage scores and sampling probabilities.

  • Leverage scores measure the relative importance of rows, and coherence is the maximum leverage score; incoherence corresponds to coherence near r/N.Each leverage score is at most 1 and their sum is r.Source: §3.2 · p. 5
  • With sampling probabilities bounded relative to leverage scores by β, the least-squares sample requirement scales as O(r/(βε)) under the stated small-ε conditions.The theorem gives an ε-accurate solution with probability at least 1 − δ.Source: §3.2 · p. 6
  • For Khatri–Rao systems, estimated leverage-score bounds avoid computing exact scores and yield β = 1/r^d; concentrated probabilities motivate combining repeated rows and including high-probability rows deterministically.The paper identifies direct computation of true leverage bounds as too expensive.Source: §3.2 · p. 6

4. Tools for sketching with concentrated sampling probabilities.

The section addresses concentrated sampling probabilities, where a small set of rows can be sampled repeatedly. It presents row combining and hybrid deterministic-random sampling as practical ways to manage that concentration.

4. Tools for sketching with concentrated sampling probabilities.

  • The sampling tools do not assume probabilities are concentrated, which matters because CP-ALS probabilities can shift from near-uniform at initialization to concentrated as factor matrices become structured.The paper reports numerical improvements from the methods in subsections 7.1 and 7.2.Source: §4 · p. 6

4.1. Combine repeated rows.

  • Combining repeated sampled rows yields a smaller system with the same sketched norm, and the experiments use this combination unless otherwise noted.The combined row accounts for the multiplicity of the sampled index.Source: §4.1 · p. 7

4.2. Hybrid deterministic and random sampling.

  • The hybrid scheme includes a subset of high-probability rows deterministically and samples the remainder randomly in proportion to their original probabilities.The deterministic set is presumed to contain the highest-probability indices, though the reweighting analysis does not depend on that choice.Source: §4.2 · p. 7

4.3. Combining rows for the hybrid deterministic and random sam-

  • The random portion of the hybrid sample can also combine repeated rows, with the random weights adjusted for the probability mass assigned to deterministic rows.The random sampling probabilities are rescaled after removing that deterministic mass.Source: §4.2 · p. 7, §4.3 · pp. 7–8

4.4. The nonviable alternative of sampling without replacement.

  • Sampling more rows until enough unique rows appear is not equivalent to sampling without replacement, because dependence makes correct weighting infeasible and can bias the sampled residual.The hybrid method includes high-probability rows first so the remaining random draws have relatively few repeats.Source: §4.3 · p. 8

4. Tools for sketching with concentrated sampling probabilities.

  • For Khatri–Rao least-squares systems, leverage-score upper bounds give a sample requirement of O(r^d/ε), and the sampling and deterministic inclusion are designed to avoid explicitly forming the matrix or all probabilities.The sampling count is independent of the number of rows and nonzeros in the tensor.Source: §4.4 · p. 8

5. Efficient leverage score sampling for KRP matrices.

The section derives leverage-score-based sampling for Khatri–Rao subproblem matrices, with a row-sample bound for accurate least-squares solutions. It also describes implicit sampling and a hybrid procedure for handling high-probability rows.

5.1. Sampling probabilities for KRP matrices and main theorem.

  • Leverage-score upper bounds for Khatri–Rao rows support sampling probabilities with a theoretical bound on β.The authors report that leverage-score sampling also outperforms uniform and two-norm-based sampling in experiments.Source: §5.1 · pp. 8–9, §5 · p. 8
  • s = O(r^d/ϵ) samples yield an ϵ-accurate residual with high probability for sufficiently small ϵ.The sample count is independent of the system’s number of rows and the tensor’s number of nonzeros.Source: §5.1 · p. 8

5.2. Implicit random row sampling for KRP matrices.

  • Sampling each factor mode independently draws Khatri–Rao rows implicitly, avoiding the O(N) cost of explicitly forming the Kronecker product of factor leverage scores.A selected row of the Khatri–Rao matrix can be assembled in O(rd) work from the corresponding factor-matrix rows.Source: §5.1 · p. 9, §5.2 · p. 9

5.3. Implicit computation of high-probability rows for deterministic inclusion.

  • The hybrid procedure identifies rows above a probability threshold for deterministic inclusion, checks a restricted set of candidate combinations, and rejection-samples the remaining rows.The candidate count can be bounded using the threshold; τ = 1/s was used in the reported experiments.Source: §5.3 · pp. 9–10

6. Alternating randomized least squares with leverage score sampling.

CP-ARLS-LEV applies leverage-score sampling to the least-squares problems in CP-ALS and builds reduced systems directly from sampled rows. The implementation also estimates fit using stratified samples to avoid the cost of exact fit calculation.

6.1. Finding indices and weights.

  • SkrpLev combines deterministic high-probability indices with rejection-sampled random indices, adjusts their weights, and combines repeated random rows.Repeated occurrences of a row scale its weight by the square root of its count.Source: §6.1 · p. 11

6.2. Full algorithm.

  • The implementation builds the sampled Khatri–Rao matrix and solves the reduced least-squares problem using the selected rows and weights.Constructing the sampled matrix costs O( s̄dr ); solving costs O( s̄r^2 + s̄nr ), with the latter work reducible using sampled-tensor sparsity.Source: §6.2 · p. 12

6.3. Efficient sampling from sparse tensor.

  • The sparse-tensor implementation stores nonzeros in coordinate format and extracts and reweights the selected unfolding entries to form the sampled right-hand side.It recommends precomputing linearized indices for each mode, requiring (d + 1) nnz(X) additional storage.Source: §6.3 · pp. 12–13

6.2. Full algorithm.

  • CP-ARLS-LEV cycles through tensor modes, sketches each least-squares problem, and groups outer iterations into epochs before checking fit.The stopping rule allows multiple epochs without sufficient fit improvement because randomized steps need not improve fit at every iteration.Source: §6.2 · pp. 11–12, §6.3 · p. 13

6.4. Estimating the fit.

  • The fit estimate uses stratified sampling of zero and nonzero tensor entries because uniform sampling would return predominantly zeros and may be inaccurate.The default sample mix is α = 0.5 nonzero entries, with zeros selected by rejection sampling.Source: §6.4 · p. 13

7. Numerical results.

Experiments compare CP-ARLS-LEV variants with CP-ALS and SPALS on large sparse tensors, examining repeat-row combination, random versus hybrid sampling, runtime and fit, and example factors. On Amazon and Reddit, randomized methods substantially reduce runtime, while fit and speed vary by sampling method and tensor.

7. Numerical results.

  • The experiments compare random and hybrid CP-ARLS-LEV, standard CP-ALS, and SPALS on Uber, Amazon, and Reddit tensors, with ranks and sample counts varied across tests.The factor matrices are initialized with standard Gaussian entries; the tested tensors range from Uber’s 3 million nonzeros to Reddit’s 5 billion.Source: §7 · p. 14

7.1. Combining repeated rows.

  • Combining sampled rows always improves solution time; it reduces sampled right-hand-side nonzeros, with smaller reductions for hybrid sampling because deterministic inclusion yields fewer repeats.The comparison uses Amazon factors of rank 10 and 217 sampled rows, averaged across 10 runs; without combining repeats or hybrid sampling, the mode-2 solve can eclipse fiber extraction.Source: Figure 2 · p. 15, §7.1 · p. 16

7.2. Comparing random and hybrid sampling.

  • At s=219, hybrid sampling achieves 10^-4 accuracy using 1.1% of the rows and about two more accuracy digits than random sampling.At that sample count, fewer than 20% of hybrid rows are deterministically included, yet they account for as much as 90% of the sampling probability; computing pdet adds cost.Source: §7.2 · pp. 16–17, Figure 3 · p. 17

7.3. Equivocal performance on small tensor.

  • For Uber, hybrid sampling improves median final fit over random at each tested sample count but is slower; at s=217, both randomized methods slightly outperform standard fit.At s=217, random is substantially faster than standard, whereas hybrid is somewhat slower.Source: §7.3 · p. 17, Figure 4 · p. 18

7.4. Order of magnitude speed improvement on massive sparse ten-

  • On Amazon, random CP-ARLS-LEV takes less than 25 minutes (13X faster than standard), while Reddit random and hybrid take 6 and 8 hours versus 4 days (>12X speedup).The Amazon comparison also reports that random’s median and best fits are approximately 0.5% worse than standard at the lower sample count; Reddit’s best fits are essentially equivalent across methods.Source: §7.4 · pp. 18, 20, Figure 5 · p. 19, Figure 6 · p. 20

7.5. Example factors for massive tensor.

  • Reddit factors include components focused on non-U.S. news, soccer and sports, and movies and television, illustrating interpretable topics in the decomposition.Component 6’s top subreddit and word factors include worldnews and world; blue bars mark rarer words or subreddits.Source: §7.5 · p. 21, Figure 7 · p. 21

8. Conclusions.

The conclusions present CP-ARLS-LEV as a leverage-score sketch for CP-ALS that applies to sparse tensors, report favorable results against SPALS on Amazon, and identify practical and theoretical directions still open.

8. Conclusions.

  • CP-ARLS-LEV sketches the overdetermined least-squares problem in CP-ALS and can be used on sparse tensors, unlike CP-ARLS-MIX, whose required FFTs destroy sparsity.The paper also provides a numerical comparison on dense tensors in Appendix C.Source: §8 · p. 21
  • On the 1.7-billion-nonzero Amazon tensor, CP-ARLS-LEV was faster and achieved a better solution than SPALS.The authors also report better sample complexity than SPALS.Source: §8 · p. 22
  • Combining repeated sampled rows shrinks the solved system, while deterministic inclusion of high-probability rows increases the number of unique samples.Estimated fit is used instead of exact fit because computing the latter costs almost as much as solving the least-squares problem; with these techniques, the implementation yields order-of-magnitude speedups on large sparse tensors.Source: §8 · pp. 22–23
  • Open questions include choosing sample counts and deterministic thresholds, bounding hybrid sampling's improvement, and finding a more robust stopping condition than estimated fit.On large tensors, low-variance fit estimates required an extremely large number of samples.Source: §8 · p. 23
  • For distributed tensors, the authors describe a possible approach that stores factor matrices on one node and gathers sampled fibers there, while leaving implementation and parallelization for future work.They contrast this with CP-ALS, whose iterations require solving a system involving the entire tensor.Source: §8 · p. 23

A.1. Properties of sketching matrix under structural conditions.

  • Theorem 10 bounds the sketched solution under two structural conditions; a tighter solution bound is available when a constant fraction of B's columns lies in A's column space.The latter assumption is described as typically reasonable when the fit is practically interesting.Source: §A.1 · p. 24

A.2. Proof that sketching matrix meets structural conditions.

The proof establishes the sketching matrix's structural conditions using leverage-score overestimates and randomized product approximation, then uses these conditions to obtain high-probability guarantees.

A.2. Proof that sketching matrix meets structural conditions.

  • With row sampling and rescaling by S using leverage-score overestimates ℓ(A), if s > 144r ln(2r/δ)/(βϵ²), then 1−ϵ ≤ σ_i²(SUA) ≤ 1+ϵ simultaneously for all i with probability at least 1−δ.Source: §A.2 · p. 24
  • For the second structural condition, the proof bounds how well the sampled product U_A^T S^T S B⊥ approximates U_A^T B⊥, which is zero by orthogonality.It applies a randomized matrix-multiplication lemma to control this approximation.Source: §A.2 · p. 25
  • Rows are sampled and rescaled by importance sampling to approximate A^T B; the proof uses an expected-error bound and Markov's inequality to control the chance that the second condition fails.The text derives a sample requirement for failure probability δ from that bound.Source: §A.2 · pp. 25–26

A.3. Main Theorem.

The main theorem combines two structural conditions for a row-sampled least-squares system. Under leverage-score overestimation assumptions, it gives a sample requirement and high-probability bounds for the sketched solution.

A.3. Main Theorem.

  • The theorem considers an overdetermined least-squares problem with a full-rank design matrix A ∈ R^(N×r), r ≪ N, and samples rows according to a probability distribution.The sample count scales with r, the overestimate parameter β, δ, and ϵ.Source: §A.3 · p. 26
  • With leverage-score probabilities bounded below by β times each row's leverage score and s = (r/β) max{C log(r/δ), 1/(δϵ)}, the structural conditions hold jointly with probability at least 1 − δ.The proof allocates failure probability δ/2 to each condition and applies a union bound.Source: §A.3 · p. 26
  • Combining the structural-condition guarantees with Theorem 10 yields high-probability bounds on the sketched solution and residual.The proof states that the joint guarantee holds with probability at least 1 − δ.Source: §A.3 · p. 26

Appendix B. Comparison of Sampling Methods.

The sampling-method comparison finds leverage-score sampling gives higher fit than length-squared sampling, while uniform sampling is not competitive. In the dense synthetic comparison, CP-ARLS-MIX outperforms CP-ARLS-LEV, although CP-ARLS-LEV remains viable.

Appendix B. Comparison of Sampling Methods.

  • Leverage-score sampling achieves higher fit than length-squared sampling, while uniform sampling is not competitive in the Uber-tensor comparison.Figure 10 compares the three sampling methods for rank-25 decomposition of the Uber tensor with 3.3 million nonzeros.Source: Appendix B · p. 26, Figure 10 · p. 27

Appendix C. Comparison on Dense Tensor.

  • On the coherent synthetic dense tensor, CP-ARLS-MIX achieves higher factor-match scores in less time than CP-ARLS-LEV; CP-ARLS-LEV is nevertheless a viable dense-tensor option.CP-ALS scores around 0.814, and CP-ARLS-MIX reaches a score of at least 0.93 in 80%, 70%, and 30% of runs at s = 2^9, 2^10, and 2^11, respectively; CP-ARLS-LEV does so in 10% of runs at s = 2^11.Source: Figure 11 · p. 27, Appendix C · p. 27
  • Figure 11 compares total runtime with the best-run factor-match score across ten shared random initializations on a 500 × 500 × 500 synthetic tensor.The tensor was engineered to have factor matrices with a few rows of concentrated leverage scores; a factor-match score of 1.0 is a perfect match.Source: Appendix C · p. 26, Figure 11 · p. 27

C.1. Data generation.

The experiments specify a synthetic tensor construction and shared run settings, then detail CP-ARLS-LEV preprocessing, per-epoch, and per-iteration costs.

C.1. Data generation.

  • A 500 × 500 × 500 rank-25 tensor uses three sparse, high-leverage columns per factor and 5% Gaussian noise.Each of the first three columns has five nonzeros, producing 15 rows with very high leverage; the planted entries have magnitude 3 with small added noise.Source: §C.1 · p. 28

C.2. Experimental setup.

  • Ten shared random initializations are used per algorithm and sample combination, with randomized runs stopping after five unimproved epochs.CP-ALS uses tolerance 10^-5 and a 250-iteration maximum; randomized runs use epoch size 5 and fit-improvement tolerance 10^-4.Source: §C.2 · p. 28

Appendix D. End-to-end Complexity of Algorithm.

  • With tensor order treated as constant, the complexity analysis tracks rank r, sample count s, fit-sample count sfit, maximum mode size n̄, and tensor nonzeros.For dense tensors, nnz(X) = O(∏k nk).Source: Appendix D · p. 28

D.1. Preprocessing.

  • Preprocessing costs O(nnz(X)) once for unfolding indices and O(r^2 n̄) to compute initial leverage scores.Source: §D.1 · p. 28

D.2. Per-epoch costs.

  • Exact fit computation costs O(r nnz(X)) per epoch, while estimated fit costs O(r sfit) for sparse tensors with billions of nonzeros.Source: §D.2 · p. 28

D.3. Per-iteration costs.

  • Each sampled solve combines index and weight generation, sampled Khatri-Rao and tensor construction, least-squares solution, and updated sampling probabilities.The stated costs include O(s log s) to combine repeated indices, O(sr^2) for QR formation, and O(sr n̄) to apply Q to right-hand sides; sparse right-hand sides can reduce the latter to O(r nnz( X̃ )).Source: §D.3 · pp. 28–29

Appendix E. Detailed Runs on Uber Tensor.

Figure 12 presents individual runs and median trajectories for the Uber tensor comparison.

Appendix E. Detailed Runs on Uber Tensor.

  • Figure 12 gives the full run data for the Uber experiments summarized in Figure 4, including CP-ARLS-LEV and standard CP-ALS.Source: Appendix E · p. 29
  • Each randomized-method marker represents an epoch of five outer iterations, while each CP-ALS marker represents one outer iteration.Source: Appendix E · p. 29
  • Dotted curves show individual runs, and solid curves show medians over interpolated runs.Source: Appendix E · p. 29

Appendix F. Initialization via Randomized Range Finder for Enron Tensor.

On Enron, randomized range finder initialization improves CP-ARLS-LEV fit relative to random initialization, reaching fit comparable to CP-ALS. Its median runtime is also faster than CP-ALS at both reported sample sizes.

Appendix F. Initialization via Randomized Range Finder for Enron Tensor.

  • The comparison uses the 4-way Enron log-count tensor with 54,202,099 nonzeros and hybrid sampling with τ = 1/s.Its entries are log(c + 1), which damps the largest counts while leaving zeros unchanged.Source: Appendix F · p. 29
  • RRF initialization forms factor-matrix columns from random linear combinations of sampled matricized-tensor fibers, tending to decrease the right-hand-side component outside the factor matrix’s column space.The method samples sinit = 105 nonzero fibers and multiplies them by a random Gaussian matrix.Source: Appendix F · pp. 29–30
  • 5.78 times faster at s = 218 and 4.39 times faster at s = 219, RRF-initialized CP-ARLS-LEV has median runtime versus CP-ALS, with fit comparable to CP-ALS.The reported fit improvement over random initialization is shown for the Enron tensor; Figure 13 presents fit and runtime comparisons.Source: Appendix F · pp. 30–31, Figure 13 · p. 31
Loading 2006.16438v3…