PaperlayerView exact revision on arXiv

Source-linked AI summary

Low Rank Approximation and Regression in Input Sparsity Time

Kenneth L. Clarkson, David P. Woodruff

arXiv:1207.6365v4cs.DS

TL;DR

The paper addresses the cost of subspace embeddings and the resulting numerical linear algebra algorithms. It introduces sparse embedding matrices based on hashing and random signs, achieving subspace embeddings computable in input-sparsity time and enabling near-linear-time algorithms across several tasks. The analysis relies on subspace structure rather than norm preservation for arbitrary large vector sets.

  • Problem

    The paper seeks subspace embeddings with small dimension and fast matrix multiplication, because prior approaches were not optimal for sparse inputs.

  • Method

    The paper uses sparse embedding matrices whose hashing structure and random signs exploit the special structure of vectors in a fixed subspace.

  • Results

    The resulting algorithms have leading time O(nnz(A)) or O(nnz(A) log n) for least-squares regression, low-rank approximation, leverage scores, and ℓp-regression.

  • Takeaways & Limitations

    Sparse subspace embeddings provide a common route to input-sparsity-time algorithms for several randomized numerical linear algebra problems.

  • Takeaways & Limitations

    The sketching dimension improves when leverage scores are roughly equal, an assumption used in applications such as matrix completion.

Abstract

from arXiv · show

We design a new distribution over $\poly(r \eps^{-1}) \times n$ matrices $S$ so that for any fixed $n \times d$ matrix $A$ of rank $r$, with probability at least 9/10, $\norm{SAx}_2 = (1 \pm \eps)\norm{Ax}_2$ simultaneously for all $x \in \mathbb{R}^d$. Such a matrix $S$ is called a \emph{subspace embedding}. Furthermore, $SA$ can be computed in $\nnz(A) + \poly(d \eps^{-1})$ time, where $\nnz(A)$ is the number of non-zero entries of $A$. This improves over all previous subspace embeddings, which required at least $Ω(nd \log d)$ time to achieve this property. We call our matrices $S$ \emph{sparse embedding matrices}. Using our sparse embedding matrices, we obtain the fastest known algorithms for $(1+\eps)$-approximation for overconstrained least-squares regression, low-rank approximation, approximating all leverage scores, and $\ell_p$-regression. The leading order term in the time complexity of our algorithms is $O(\nnz(A))$ or $O(\nnz(A)\log n)$. We optimize the low-order $\poly(d/\eps)$ terms in our running times (or for rank-$k$ approximation, the $n*\poly(k/eps)$ term), and show various tradeoffs. For instance, we also use our methods to design new preconditioners that improve the dependence on $\eps$ in least squares regression to $\log 1/\eps$. Finally, we provide preliminary experimental results which suggest that our algorithms are competitive in practice.

1 Introduction

The paper targets randomized numerical linear algebra algorithms whose running times are not input-sparsity optimal. It introduces sparse subspace embeddings and uses them to obtain near-linear-time algorithms for regression, low-rank approximation, leverage scores, and ℓp-regression.

  • Motivation: Existing randomized methods address least-squares regression, low-rank approximation, and leverage-score computation, but often require more than input-sparsity time.These problems support applications including data mining, recommendation systems, information retrieval, web search, clustering, and mixture learning.
  • Motivation: O(nnz(A)) + poly(d/ε) time is the desired regression complexity, while prior low-rank methods required O(nnz(A)(k/ε + k log k) + n·poly(k/ε)).The gap is especially relevant when nnz(A) is much smaller than nd.
  • Results: The paper achieves algorithms whose leading time term is O(nnz(A)), sometimes with a logarithmic factor, for several numerical linear algebra problems.The results cover least-squares regression, low-rank approximation, approximate leverage scores, and ℓp-regression.
  • Results: Least-squares algorithms include running times O(nnz(A) log(n/ε) + r^3 log^2 r + r^2 log(1/ε)) and O(nnz(A)) + ˜O(d^3ε^-2).The paper also gives generalized regression in O(nnz(A) log n + r^2((r + d′)ε^-1 + rd′ + r log^2 r + log n)).
  • Results: Low-rank approximation runs in O(nnz(A)) + n·poly(k log(n)/ε), while all leverage scores and ℓp-regression use O(nnz(A) log n) plus polynomial lower-order terms.The leverage-score result simultaneously (1+ε)-approximates all n scores for constant ε.
  • Techniques: A subspace embedding preserves the norms of every vector in the column space after applying a randomized matrix S.The paper designs such embeddings with t = poly(d/ε) and computes SA in input-sparsity time, then plugs them into existing algorithms.

2 Sparse Embedding Matrices

A sparse embedding matrix is constructed by hashing input coordinates into t buckets and randomly signing them. The resulting matrix has one nonzero entry per input column.

  • Construction: For parameter t, the construction defines a random map h:[n]→[t] that assigns each input coordinate uniformly to a bucket.The map is used to determine where each coordinate contributes in the sketch.
  • Construction: The binary matrix Φ has Φ_h(i),i = 1 and zeros elsewhere, so every column contains exactly one nonzero entry.Thus Φ aggregates coordinates according to the hash map.
  • Construction: The diagonal matrix D independently assigns each coordinate a +1 or −1 sign with equal probability.The signs randomize contributions that collide in the same bucket.
  • Construction: The product ΦD is called a sparse embedding matrix.This is the paper’s basic sketching matrix for preserving norms on a subspace.

3 Analysis

The analysis separates large and small coordinates, controls hashing collisions and random-sign cross terms, and then lifts fixed-vector norm preservation to the entire subspace using a net argument.

  • Setup: The rows of an orthonormal basis U are ordered by leverage magnitude, and thresholds define the heavy-coordinate index s and related tail parameters.These parameters organize the analysis according to the concentration of leverage scores.
  • Tail analysis: For the tail coordinates, Bernstein-type bounds and bucket concentration control the hashed mass in each bucket with high probability.The event E_h is obtained by applying a concentration bound to each bucket and taking a union bound.
  • Tail analysis: Random signs and the Hanson–Wright inequality bound deviations of the hashed tail norm for a fixed unit vector in the column space.The hashed norm is represented as a quadratic form in the independent signs.
  • Heavy coordinates: The first s coordinates are handled directly by requiring that they do not collide, an event whose failure probability is at most s^2/t.This isolates the potentially large coordinates from the collision analysis.
  • Cross terms: Cross terms between heavy and light coordinates are bounded separately, completing fixed-vector norm preservation under the required events.The analysis combines bounds for large entries, small entries, and their cross terms.
  • Subspace extension: A net argument converts fixed-vector preservation into preservation for every vector in the r-dimensional subspace.The argument bounds the quadratic form on a finite net and extends the bound to the whole unit sphere.
  • Main theorem: Theorem 11 gives t = O((r/ε)^4 log^2(r/ε)) and success probability at least 9/10, while applying ΦD to A takes O(nnz(A)) time.The theorem also gives a smaller dimension under additional leverage-score conditions.

4 Partitioning Leverage Scores

The analysis partitions leverage scores into constant-factor groups and bounds errors from within-group collisions and cross-group collisions. These bounds yield a subspace embedding with dimension scaling as O(P(r/ε)^2 log^6(r/ε)).

  • Leverage-score partitioning: The method groups leverage scores above T into q = O(log(r/ε)) constant-factor classes and analyzes collisions within each class.The analysis leaves poly(log(r/ε)) factors unoptimized.
  • Within-group collisions: With probability 1 − o(1), each leverage-score group’s colliding rows have squared spectral norm O(q(β_j + 1/n_j + r/t)).This is established for all groups simultaneously.
  • Within-group errors: Sampling sufficiently many rows from each group preserves every component’s norm with error O(ε/q), giving total within-group error O(ε).The argument uses a net over the collision-supported subspaces.
  • Within-group errors: Choosing t = O(P(r/ε)^2 log^6(r/ε)) controls the within-group contribution with failure probability at most δ1 + O(1/log r).The theorem allows error proportional to ε/(Pδ1).
  • Cross terms: Cross-group collisions are bounded by controlling bins containing items from multiple groups, with Pr[E] ≥ 1 − δ1.A separate balls-and-bins event bounds the number of items per bin with probability at least 1 − 1/r.
  • Combined guarantee: Combining within-group and cross-term bounds gives ∥Sy∥ = 1 ± C′ε/(Pδ1) for all unit y ∈ C(A), using t = O(P(r/ε)^2 log^6(r/ε)).The failure probability is at most δ1 + O(1/log r).

5 Generalized Sparse Embedding Matrices

The generalized sparse embedding combines random hashing with independent block embeddings to handle matrices whose leverage scores are uneven. It achieves a subspace embedding using t = O(rε^-4 log(r/εδ)(r + log(1/εδ))) rows while applying S in input-sparsity time.

  • Construction: Rows are randomly hashed into q groups, independently transformed by block matrices B^(i), permuted, and combined into the generalized embedding S.The map is defined as a block-diagonal matrix composed with the permutation.
  • Implementation: S · A can be computed in O(nnz(A)(log(r/εδ))/ε) time because each nonzero is processed by one sparse block.The permutation costs O(nnz(A)), and each block multiplication costs O(a).
  • Small-entry vectors: Each column-position block N_j has the distribution of a sparse embedding matrix with qv = t/a rows and n columns.It has one random ±1 nonzero per column, with locations induced by the hash and permutation.
  • Small-entry vectors: For the low-leverage portion, the construction preserves all restricted subspace norms with failure probability at most δ/r.This follows by applying the low-dimensional subspace guarantee to each hash group and taking a union bound.
  • Final guarantee: Theorem 28 gives an embedding dimension t = O(rε^-4 log(r/εδ)(r + log(1/εδ))) and application time O(nnz(A)ε^-1 log(r/δ)).With probability at least 1 − δ, all y ∈ C(A) satisfy ∥Sy∥2 = (1 ± ε)∥y∥2.

6 Approximating Leverage Scores

The section develops a faster leverage-score approximation procedure using the paper’s sparse embedding, achieving simultaneous approximation with near input-sparsity running time.

  • For general ε, the running time becomes O(nnz(A) log n) + poly(rε^-1 log n).
  • For constant ε, all n leverage scores are simultaneously (1+ε)-approximated in O(nnz(A) log n + r^3 log^2 r + r^2 log n) time.The algorithm succeeds with probability at least 2/3.
  • The algorithm first identifies the rank and independent columns, then computes a change of basis and randomized row-norm estimates.The change-of-basis construction makes the row norms correspond to leverage scores up to relative error.
  • The procedure replaces the earlier embedding with Π1 = F·S, combining the sparse embedding S with a fast Johnson–Lindenstrauss transform F.S·A costs O(nnz(A) log r), while F·(SA) costs O(r^3 log^2 r).
  • The embedding composition preserves norms on the column space with probability at least 7/10, and the overall leverage-score algorithm succeeds with probability at least 2/3.

7 Least Squares Regression

The section applies sparse, sampling-based, and composed embeddings to least-squares regression, obtaining relative-error algorithms whose leading costs are linear or nearly linear in the input sparsity.

  • ℓ2 regression achieves a (1+ε)-factor approximation with probability at least 2/3 in O(nnz(A) + poly(dε^-1)) time.The stated bound is O(nnz(A) + O(d^3ε^-2 log^7(d/ε)).
  • 7.4 Generalized Regression: Algorithm: Generalized regression with B having multiple right-hand sides is solved to relative error ε in O(nnz(A) log n + r^2(rε^-1 + rd′ + r log^2 r + d′ε^-1 + log n)) time.
  • 7.4 Generalized Regression: Algorithm: The generalized regression algorithm also yields a coreset of size O(r(ε^-1 + log r)).
  • 7.1 Generalized Regression and Affine Embeddings: Sparse embeddings, leverage-score sampling, and SRHTs provide different affine-embedding dimensions for regression analyses.The sparse embedding requires t = O(ε^-2r^2 log^6(r/ε)), while leverage-score sampling requires t = O(ε^-2r log r) for a weak affine embedding.
  • 7.7 Iterative Methods for Regression: Iterative preconditioning reduces the dependence on ε to logarithmic: O(log(1/ε)) iterations suffice for relative error ε.The resulting work is O((nnz(A) + r^2) log(1/ε)) after reducing the matrix to rank r.
  • 7.7 Iterative Methods for Regression: The iterative method is comparable to conjugate-gradient methods when d′ < r, but is asymptotically worse than normal equations when d′ > r.

8 Low Rank Approximation

The section reduces low-rank approximation to generalized regression and uses composed embeddings to obtain a relative-error rank-k approximation with input-sparsity-dependent running time.

  • The target is a low-rank matrix whose Frobenius distance from A is within a factor 1+ε of the best rank-k error Δ_k.
  • The main embedding composition is justified because products of subspace embeddings and matrix-product approximations preserve the corresponding properties up to O(ε) error.
  • For A ∈ R^n×n, the algorithm finds rank-k factors L, D, W with orthonormal columns or diagonal structure such that ||A−LDW^T|| ≤ (1+ε)Δ_k.The failure probability is 1/10.
  • The construction computes a basis from AR^T, sketches U and A with a sparse embedding composed with an SRHT, and solves a sketched rank-k problem.
  • The algorithm’s dominant costs include O(nnz(A)) plus terms such as Õ(nk^2ε^-4), Õ(k^3ε^-5), and Õ(nk^2ε^-2).
  • The resulting factorization satisfies ||LDW^T−A|| ≤ (1+ε)^2Δ_k.

9 ℓp-Regression for any 1 ≤p < ∞

The section extends the embedding-based approach to ℓp regression for every constant p ∈ [1,∞), using well-conditioned bases and sampling to construct relative-error coresets and solutions.

  • For every constant p ∈ [1,∞), the algorithm constructs a coreset and solution achieving relative error 1+ε with probability at least 1/2.The entire construction runs with probability 1−o(1) within the stated running-time guarantee.
  • The construction divides A into row blocks and applies generalized sparse embeddings to the blocks before forming the sampling structure.
  • A well-conditioned basis is obtained by applying a sparse embedding and changing the basis of the resulting column space.The basis construction uses Condition(A), including a preliminary projection and basis transformation.
  • The basis-construction time is O(nnz(A) log n) + poly(rε^-1), replacing a leading O(nr^5 log n) term by an input-sparsity-dependent term.
  • The method solves ℓp regression by sampling rows according to norms of the corresponding well-conditioned basis.A second projection estimates these row norms for efficient sampling.
  • The overall complexity is O(nnz(A) log n) + poly(rε^-1) when r < log n.

10 Preliminary Experiments

Preliminary experiments evaluated a simplified low-rank approximation method on 1,155 sparse matrices, finding promising practical behavior relative to the paper’s general bounds. The experiments varied the sketch sizes and assessed approximation error and conditioning.

  • Results: The preliminary experiments suggest that the simplified low-rank approximation technique may perform much better in practice than the general theoretical bounds.This is presented as an initial empirical indication rather than a definitive guarantee.
  • Experimental setup: The tested method skipped the randomized Hadamard transform and combined a sparse embedding with leverage score sampling.Its Frobenius error was compared with that of the best rank-k approximation.
  • Experimental setup: The experiments varied t_R, the number of columns of R-hat, while restricting t_R ≤ d/5.The tested values followed floor(1.6^z − 0.5) for integer z ≥ 1.
  • Experimental setup: The number of sampled rows t_S was selected so that cond(SU) was at most 1.2.Because U has orthogonal columns, sufficiently large leverage-score samples achieve this conditioning target.
  • Experimental setup: 1,155 sparse matrices from 70 application-oriented sub-collections were tested, with n up to about 7000 and at most 10^5 nonzeros.The matrices came from the University of Florida Sparse Matrix Collection.
  • Results: The Figure 1 Pareto curve summarizes error across tested sketch sizes, allowing at most 1% of points to violate each plotted tradeoff.Figure 2 similarly indicates that t_R/t_S need be no smaller than about 1/110 for cond(SU) ≤ 1.2.

A Deferred proofs

The deferred proofs establish norm and approximation bounds through moment inequalities, concentration, subspace-embedding properties, and norm comparisons. They also derive the stated runtime for constructing a well-conditioned basis in the ℓ_p-regression analysis.

  • Deferred proofs: A subspace embedding implies that A^T S^T S A − I has operator norm at most ε_0^2, which bounds the approximation error term.The proof uses this spectral bound together with Frobenius-norm subadditivity and ∥CD∥_F ≤ ∥C∥_2∥D∥_F.
  • Deferred proofs: The sparse embedding norm analysis uses a hash-function representation and the second-moment method to bound ∥SA∥_F^2.The construction writes S as ΦD and groups entries according to the hash function.
  • Deferred proofs: Chebyshev’s inequality yields the lemma’s bound when t = Ω(ε^-2).The argument follows bounds on the relevant terms in the preceding moment calculation.
  • Deferred proofs: Hoeffding’s inequality gives failure probability at most 2 exp(−2ε^2t/α^2), so t = Θ(ε^-2(log n)^2) suffices for failure probability at most 1/10.The argument assumes the relevant row bound with α = O(log n).
  • Deferred proofs: For ℓ_p regression, the proof applies a block decomposition, relates 2-norm and p-norm quantities, and uses a well-conditioned basis to obtain simultaneous bounds.The resulting basis computation takes O(nnz(A) log n) time.
Loading 1207.6365v4…