PaperlayerView exact revision on arXiv

Source-linked AI summary

Randomized Block Krylov Methods for Stronger and Faster Approximate Singular Value Decomposition

Cameron Musco, Christopher Musco

arXiv:1504.05477v4cs.DScs.LGmath.NA

TL;DR

The paper addresses the cost and limitations of randomized SVD methods, including gap-dependent guarantees, inadequate spectral-norm behavior, and the distinction between low-rank approximation and PCA quality. It introduces randomized block Krylov methods and proves stronger accuracy with fewer iterations, while analyzing practical runtime improvements.

  • Problem

    Randomized SVD methods need nearly optimal low-rank approximation and PCA without expensive computation, singular-value-gap dependence, or spectral-norm and principal-component limitations.

  • Method

    The paper analyzes Simultaneous Iteration and a randomized block Krylov method, with guarantees based on subspaces that align with leading and intermediate singular vectors.

  • Results

    Both algorithms achieve low-rank approximation and PCA guarantees, while Block Krylov Iteration reduces the iteration requirement to q = O(log d/√ε) and both run in O(nnz(A)kq) time.

  • Takeaways & Limitations

    The results support using block Krylov and modified Simultaneous Iteration as gap-independent alternatives to sketch-and-solve methods for approximation and PCA.

  • Takeaways & Limitations

    Small-block methods remain constrained by duplicate singular values and large singular-value clusters, which introduce dependence on gaps between adjacent top-k singular values.

Abstract

from arXiv · show

Since being analyzed by Rokhlin, Szlam, and Tygert and popularized by Halko, Martinsson, and Tropp, randomized Simultaneous Power Iteration has become the method of choice for approximate singular value decomposition. It is more accurate than simpler sketching algorithms, yet still converges quickly for any matrix, independently of singular value gaps. After $\tilde{O}(1/ε)$ iterations, it gives a low-rank approximation within $(1+ε)$ of optimal for spectral norm error. We give the first provable runtime improvement on Simultaneous Iteration: a simple randomized block Krylov method, closely related to the classic Block Lanczos algorithm, gives the same guarantees in just $\tilde{O}(1/\sqrtε)$ iterations and performs substantially better experimentally. Despite their long history, our analysis is the first of a Krylov subspace method that does not depend on singular value gaps, which are unreliable in practice. Furthermore, while it is a simple accuracy benchmark, even $(1+ε)$ error for spectral norm low-rank approximation does not imply that an algorithm returns high quality principal components, a major issue for data applications. We address this problem for the first time by showing that both Block Krylov Iteration and a minor modification of Simultaneous Iteration give nearly optimal PCA for any matrix. This result further justifies their strength over non-iterative sketching methods. Finally, we give insight beyond the worst case, justifying why both algorithms can run much faster in practice than predicted. We clarify how simple techniques can take advantage of common matrix properties to significantly improve runtime.

1 Introduction

The paper motivates randomized SVD methods that avoid singular-value-gap dependence while targeting spectral-norm approximation and PCA quality. It contrasts iterative methods with faster sketching approaches whose guarantees can be inadequate for heavy-tailed data.

  • Traditional SVD algorithms typically require O(nd^2) time, motivating randomized techniques for low-rank approximation and PCA.
  • Singular-value-gap-dependent runtime bounds become ineffective when nearby singular values are common, because classical analyses track individual singular-vector convergence.
  • For low-rank approximation and PCA, capturing nearly the variance of the top singular vectors avoids needing to distinguish closely spaced singular values.
  • Heavy singular-value tails can make Frobenius-norm error insufficient, allowing good multiplicative error without alignment to large singular vectors.
  • Spectral-norm guarantees are stronger than Frobenius-norm guarantees because they control recovery up to the tail-noise threshold.
  • Random-start Simultaneous Power Iteration achieves spectral-norm guarantees after approximately ˜O(1/ε) iterations and has become a standard randomized SVD method.

2 Our Results

The paper establishes gap-independent guarantees for Block Krylov Iteration, improving the iteration complexity of Simultaneous Iteration while also providing stronger per-vector PCA guarantees. Experiments show Block Krylov Iteration performs substantially better in practice, and gap-dependent analyses can further accelerate both methods on matrices with decaying singular values.

  • Faster Algorithm: The analysis gives the first gap-independent bound for a Krylov subspace method.Earlier theoretical results for Krylov and Lanczos methods were substantially more limited or depended on singular value gaps.
  • Stronger Guarantees: (1+ε) spectral-norm low-rank approximation can fail to ensure accurate principal components, so the paper targets per-vector variance capture instead.For some matrices, even spectral and Frobenius low-rank approximation errors can remain good while the returned approximate principal components have significantly lower quality.
  • Stronger Guarantees: On SNAP/AMAZON0302, sketch-and-solve achieves errors below 1.001 times optimal in Frobenius norm and 1.038 times optimal in spectral norm, yet returns poor per-vector error.This illustrates why low-rank approximation metrics alone may be insufficient for evaluating principal components.
  • Stronger Guarantees: Both Block Krylov Iteration and a modified Simultaneous Iteration achieve the stronger per-vector PCA guarantee in gap-independent runtimes.The guarantee requires each approximate singular vector to capture nearly as much variance as the corresponding true singular vector and does not require convergence to that vector when gaps are small.
  • Faster Algorithm: O(log d/√ε) iterations give Block Krylov Iteration the same low-rank approximation guarantees as Simultaneous Iteration, which requires O(log d/ε) iterations.Both algorithms run in time O(nnz(A)kq), excluding lower-order terms.
  • Beyond Worst Case: Gap-dependent analysis can be tighter when ε is below the relative singular-value gap, while gap-independent bounds better predict iteration counts for high-dimensional data problems.The algorithms can also exploit decaying singular values and use a wider initialization block to accelerate convergence.

3 Background and Intuition

Randomized low-rank approximation methods trade expensive full SVD computation for sketches or iterative subspace construction, while spectral-norm accuracy requires suppressing the singular-value tail. Block Krylov Iteration uses lower-degree polynomials to achieve comparable denoising with fewer iterations, although its analysis must address how to compute a suitable approximation within the Krylov subspace.

  • Randomized Approximation: Randomized sketch-and-solve methods compress A into a small matrix, enabling fast, parallel, and pass-efficient low-rank approximation.Sparse Johnson–Lindenstrauss methods can reduce Frobenius-error approximation to O(nnz(A)+n poly(k/ε)) time.
  • Spectral-Norm Accuracy: Spectral norm error is harder than Frobenius error because lower singular values corrupt sketches when the singular-value tail is large.Simultaneous Iteration addresses this by powering A, which preserves singular vectors while separating singular values.
  • Spectral-Norm Accuracy: Powering makes singular values above (1 + ε)σk+1 much larger than values at or below σk+1, effectively denoising the spectrum for accurate approximation.The resulting basis aligns well with singular vectors whose singular values exceed the threshold.
  • Krylov Acceleration: Block Krylov Iteration constructs polynomial transforms of A from a Krylov subspace, allowing lower-degree polynomials than Aq to suppress tail singular values.Chebyshev polynomials are used because they can push the tail near zero with lower long-run growth.
  • Krylov Acceleration: The analysis shows that the best Frobenius-error rank-k approximation within the Krylov span suffices to obtain a near-optimal spectral-norm approximation.This resolves the obstacle that the best spectral-norm approximation within a fixed subspace is not known to be efficiently computable.
  • Per-Vector Guarantees: Per-vector guarantees are strongest when σk ≥ (1 + ε)σk+1, whereas small gaps can leave intermediate singular vectors outside the computed basis.For Frobenius approximation, the loss is bounded by ε∥A − Ak∥2 and depends on singular values near the threshold.

4 Preliminaries

The preliminaries define the SVD-based optimal rank-k approximation and the subspace-restricted approximations used by the algorithms. They also introduce randomized sketches and Chebyshev polynomials as the main analytical tools.

  • SVD and Low-Rank Approximation: The SVD decomposes A into orthonormal singular vectors and nonnegative singular values, with partial SVD algorithms returning the top k vectors.The rank-k optimum is formed by retaining the largest k singular values.
  • SVD and Low-Rank Approximation: The optimal rank-k approximation Ak minimizes both Frobenius and spectral norm error among all rank-k matrices.The remainder A − Ak represents the singular-value tail beyond rank k.
  • Subspace Approximations: Given an orthonormal basis Q, Frobenius-optimal approximation within its span is obtained by projecting A to Q and taking a rank-k approximation.The computation uses the small matrix M = QT(AAT)Q and its eigendecomposition or SVD.
  • Linear Algebra Tools: An orthonormal basis for span(M) is a matrix Q whose projection QQT fully recovers the columns of M.QQT is an orthogonal projection because it is symmetric and idempotent.
  • Randomized Sketches: A Gaussian sketch AΠ yields a basis Z for its column span with high probability, providing the randomized low-rank approximation primitive used in the proofs.This block-sketch result replaces earlier single-vector random-initialization analyses for block methods.
  • Polynomial Approximation: Chebyshev minimizing polynomials provide degree-q transforms tailored to a threshold α and gap γ for suppressing unwanted singular values.For odd q, the polynomial contains only odd-powered monomials.

5 Implementation and Runtimes

The implementations build Krylov or powered subspaces from randomized starting matrices, stabilize them through orthonormalization, and extract singular directions from a small projected matrix. Block Krylov uses a larger subspace and achieves its runtime with q = Θ(log d/√ε), compared with Simultaneous Iteration’s q = Θ(log d/ε).

  • Implementation: Algorithm 1 permits alternative random starting matrices and extra columns p > k to improve approximation accuracy.Its output basis can then be post-processed for the desired approximation or PCA guarantee.
  • Implementation: For Block Krylov Iteration, setting Z = QŪk produces the best rank-l Frobenius approximation within the Krylov span for every l ≤ k.If only near-optimal low-rank approximation is needed, Z can instead be set to Q.
  • Numerical Stability: Orthonormalizing Krylov blocks after each iteration preserves the exact-arithmetic column span while substantially improving conditioning.This stabilizes computations because powering can make K poorly conditioned.
  • Simultaneous Iteration: Simultaneous Iteration computes K using repeated multiplications by A and AT, then orthonormalizes K and performs a small projected SVD.Its stated iteration choice is q = Θ(log d/ε).
  • Block Krylov Iteration: Block Krylov Iteration can compute Q through Block Lanczos recurrence or explicit Krylov construction followed by QR decomposition.The explicit approach avoids some stability issues, while the recurrence yields a block-tridiagonal projected matrix.
  • Block Krylov Iteration: Block Krylov Iteration uses a kq-dimensional subspace and sets q = Θ(log d/√ε), with additional costs from its larger projected matrix.The runtime includes terms scaling with nnz(A)kq, n(kq)^2, and (kq)^3.

6 Error Bounds

The main approximation lemma establishes Frobenius, spectral-norm, and per-vector guarantees for both algorithms with probability 99/100. Its three properties explain accuracy for dominant, intermediate, and all relevant singular-vector ranges.

  • Main Approximation Lemma: With probability 99/100, Algorithms 1 and 2 return a basis satisfying the paper’s low-rank approximation and per-vector guarantees.The guarantees are formalized through the main approximation lemma and subsequent theorems.
  • Property 1: Property 1 shows that the returned subspace aligns closely with the top m singular vectors, where σ_i ≥ (1 + ε/2)σ_{k+1}.The polynomial separates singular values above this threshold from those below σ_{k+1}.
  • Property 3: Property 3 bounds total error using w, the number of singular values in the intermediate range near σ_k.This property is critical for near-optimal Frobenius-norm low-rank approximation.
  • Property 2: A rank-k subspace within the Krylov space contains vectors that align nearly as well with A as the kth singular vector.This construction ensures good performance even for singular vectors in the intermediate range.
  • Algorithmic Guarantees: Theorems 10–12 transfer the lemma’s properties to near-optimal spectral-norm, Frobenius-norm, and per-vector guarantees.The same conclusions apply to both Simultaneous Iteration and Block Krylov Iteration.

7 Improved Convergence With Spectral Decay

When the spectrum decays rapidly, initializing with more than k random vectors yields substantially faster convergence. The resulting bounds depend on spectral decay rather than adjacent singular-value gaps.

  • Practical Acceleration: The modified initialization can significantly accelerate both methods when singular values decay rapidly.The analysis focuses on Block Krylov Iteration, while the arguments extend to Simultaneous Iteration.
  • Spectral Decay: For rapidly decaying spectra with σ_{p+1} ≤ c · σ_k, the ε dependence becomes logarithmic rather than polynomial.Here p ≥ k need not be much larger than k.
  • Practical Acceleration: Larger intermediate subspaces can have limited additional runtime because matrix multiplication dominates post-processing and memory costs.The paper notes this can hold even for 2k or 10k starting vectors.
  • Practical Acceleration: The resulting convergence varies with σ_k/σ_{p+1} rather than the gap between adjacent singular values.This retains the approximation-focused advantage of avoiding precise singular-vector convergence.

8 Experiments

Experiments on three datasets confirm that Block Krylov Iteration converges faster than Simultaneous Iteration for spectral-norm and per-vector accuracy. Both methods achieve very accurate Frobenius error quickly and can outperform worst-case predictions.

  • Frobenius Error: Both algorithms obtain very accurate relative Frobenius-norm error with very few iterations.Heavy singular-value tails allow performance beyond the worst-case analysis.
  • Convergence Results: Block Krylov Iteration converges much more rapidly than Simultaneous Iteration for spectral-norm low-rank approximation and per-vector error.This matches the theoretical analysis.
  • Convergence Results: Nearly optimal error is often reached in fewer than 8 Block Krylov iterations, whereas Simultaneous Iteration can take much longer to reach 1% error.The comparison concerns spectral-norm and per-vector metrics.
  • Runtime: On the 20 NEWSGROUPS dataset, Block Krylov Iteration outperforms Simultaneous Iteration for small ε.Its additional memory and post-processing costs are small relative to each iteration’s large matrix multiplication.
  • Spectral Gaps: Observed singular-value gaps are small: σ_k/σ_{k+1} − 1 equals .004 for AMAZON0302 and .011 for 20 NEWSGROUPS, versus .042 for EMAIL-ENRON.These values challenge treating singular-value gaps as constant, even for small ε.

Frobenius Norm Low-Rank Approximation

The Frobenius-norm analysis begins with a deterministic subspace lemma and applies Gaussian sketching to show that the resulting basis supports relative-error low-rank approximation.

  • Deterministic Lemma: A deterministic lemma analyzes an orthonormal basis for the column span of AS when S has rank k.This lemma is used as the foundation for the Frobenius-norm result.
  • Gaussian Sketching: For Gaussian Π ∈ R^{d×k}, an orthonormal basis Z for span(AΠ) satisfies the Frobenius-norm low-rank approximation guarantee with probability at least 99/100.The probability can be increased by changing a constant, affecting only logarithmic terms in applications.
  • Conclusion: The argument establishes a relative Frobenius-norm approximation from the Gaussian sketch’s captured subspace.The proof combines the sketching lemma with norm inequalities.
  • Proof Strategy: The proof applies the deterministic lemma with S = Π and uses Gaussian rotational invariance and matrix concentration.Spectral submultiplicativity controls the resulting norm terms.

Chebyshev Polynomials

The lemma constructs a degree-q polynomial by scaling a Chebyshev polynomial, then establishes its derivative and boundedness properties using standard Chebyshev identities. Choosing odd q additionally ensures that the polynomial contains only odd-degree terms.

  • Construction: A required degree-q polynomial is constructed from the standard Chebyshev polynomial T_q(x), defined by a three-term recurrence.The construction is well defined because T_q(x) is positive for arguments greater than 1.
  • Chebyshev properties: T_q(x) is bounded by 1 on [-1,1], while its closed-form behavior above 1 supports the polynomial’s growth analysis.These standard properties provide the basis for proving the lemma’s boundedness and derivative conditions.
  • Polynomial properties: The scaled polynomial satisfies property 1, and proving property 2 reduces to showing that its derivative is at least 1 for x ≥ (1 + γ)α.The derivative is analyzed by the chain rule and standard formulas for T′_q.
  • Derivative analysis: The derivative formula expresses T′_q as a positive combination of lower-degree Chebyshev polynomials, with an additional q term when q is odd.Since the relevant Chebyshev terms are positive above 1, the derivative remains nonnegative and the required inequality extends beyond the threshold.
  • Polynomial properties: For x > (1 + γ), the derivative-related expression does not decrease, completing the proof of property 2; boundedness on [0, α] is handled using T_q(x/α) ≤ 1.The proof invokes the standard bound for Chebyshev polynomials on [-1,1].
  • Parity: When q is odd, the scaled Chebyshev polynomial contains only odd-powered monomials.This follows from the parity structure of odd-degree Chebyshev polynomials.

Additive Frobenius Norm Error Implies Additive Spectral Norm Error

This lemma bounds the spectral-norm error of a rank-k approximation using the optimal rank-k approximation error of A. The proof applies Weyl’s monotonicity theorem after decomposing A into its approximation error and rank-k component.

  • Statement: For any matrix A and rank-k matrix B, the lemma relates ||A − B||_2 to the optimal rank-k approximation error of A.The result is stated for B satisfying the lemma’s rank and error condition.
  • Proof: Weyl’s monotonicity theorem is applied to the decomposition A = (A − B) + B.The resulting singular-value inequality compares singular values of A with those of the residual and B.
  • Proof: Because B has rank k, its (k+1)-st singular value is zero.This rank property removes the corresponding term from the Weyl-based inequality.
  • Dimensional case: The argument also covers n < d by applying the same reasoning to A^T and B^T.Transposition puts the matrices into the dimensional setting used by the theorem.
  • Conclusion: The squared top singular value of A − A_k equals ||A − A_k||_2^2, yielding the lemma’s final spectral-norm bound.This identity converts the singular-value inequality into the stated approximation-error expression.
Loading 1504.05477v4…