PaperlayerView exact revision on arXiv

Source-linked AI summary

Subspace Iteration Randomization and Singular Value Problems

Ming Gu

arXiv:1408.2208v1math.NA

TL;DR

The paper addresses how to compute reliable fixed-rank approximations and singular values without the high cost of a full SVD. It analyzes randomized algorithms through subspace iteration and shows that, with very high probability, they yield accurate, rank-revealing approximations and singular values, especially for rapidly decaying spectra.

  • Problem

    Full SVD computation is costly, while existing randomized analyses do not establish convergence of computed rank-k approximations to the truncated SVD or jointly characterize approximation and singular-value accuracy.

  • Method

    The paper develops a novel error analysis of randomized algorithms within the subspace iteration framework, using Gaussian start matrices and linear-algebraic and statistical tools.

  • Results

    With failure probability at most 6e^-p, randomized algorithms can provide very good low-rank approximations when σk+1 ≪ ∥A∥2, with rapidly decaying spectra also yielding accurate singular values.

  • Takeaways & Limitations

    The computed low-rank approximations are rank-revealing, and randomized algorithms can efficiently and reliably estimate condition numbers and matrix 2-norms in the supported cases.

  • Takeaways & Limitations

    The analysis focuses on the 2-norm and Frobenius norm and notes that subspace iteration can converge slowly without over-sampling or favorable starting matrices.

Abstract

from arXiv · show

A classical problem in matrix computations is the efficient and reliable approximation of a given matrix by a matrix of lower rank. The truncated singular value decomposition (SVD) is known to provide the best such approximation for any given fixed rank. However, the SVD is also known to be very costly to compute. Among the different approaches in the literature for computing low-rank approximations, randomized algorithms have attracted researchers' recent attention due to their surprising reliability and computational efficiency in different application areas. Typically, such algorithms are shown to compute with very high probability low-rank approximations that are within a constant factor from optimal, and are known to perform even better in many practical situations. In this paper, we present a novel error analysis that considers randomized algorithms within the subspace iteration framework and show with very high probability that highly accurate low-rank approximations as well as singular values can indeed be computed quickly for matrices with rapidly decaying singular values. Such matrices appear frequently in diverse application areas such as data analysis, fast structured matrix computations and fast direct methods for large sparse linear systems of equations and are the driving motivation for randomized methods. Furthermore, we show that the low-rank approximations computed by these randomized algorithms are actually rank-revealing approximations, and the special case of a rank-1 approximation can also be used to correctly estimate matrix 2-norms with very high probability. Our numerical experiments are in full support of our conclusions.

1. Introduction.

The paper develops randomized low-rank approximation within the subspace iteration framework, targeting accurate approximations and singular values at lower computational cost than a full SVD. Its analysis shows high-probability accuracy, rank-revealing behavior, and rapid convergence when singular values decay quickly.

  • Motivation: The truncated SVD is optimal for fixed-rank approximation, but computing a full SVD costs O(mn^2) flops and can be prohibitively expensive for large sparse matrices.The paper focuses on the 2-norm and Frobenius norm, where the truncated SVD attains the smallest possible errors.
  • Existing methods: Subspace iteration is efficient on modern architectures but typically needs strong singular-value separation and a good starting matrix for fast convergence.Its main computations use matrix-matrix products and QR factorizations, while the start matrix can critically affect convergence.
  • Approach: Gaussian starting matrices support the analysis through rotational invariance and established singular-value probability-density results.A rotated standard Gaussian matrix remains standard Gaussian, simplifying the probabilistic analysis.
  • Approach: The paper analyzes randomized algorithms within the subspace iteration framework to obtain both low-rank approximation-error bounds and singular-value accuracy results.This combines theories and technical machinery from randomized algorithms and subspace iteration.
  • Results: With high probability, computed leading singular values remain a good fraction of the true values and converge quickly under rapid decay, while the resulting approximations are rank-revealing.The paper also identifies an alternative starting matrix and reports that its numerical experiments support these conclusions.

2. Algorithms.

The paper develops classical and randomized subspace iteration algorithms for computing low-rank approximations, emphasizing how the starting matrix affects convergence. Randomization is intended to avoid poor starts and can make a single iteration effective, while SVD truncation trades accuracy and cost against rank control.

  • Basic Subspace Iteration: Algorithm 2.1 iteratively builds a subspace from A(A^T A)^qΩ, compresses A to B = Q^T A, truncates B's SVD, and returns QB_k.The input uses an n × ℓ start matrix Ω and produces a rank-k approximation.
  • Basic Subspace Iteration: Its per-iteration cost depends linearly on the target subspace dimension and matrix dimensions, while larger ℓ can improve convergence enough to offset extra work.The trade-off is governed by the singular-value ratio involving σ_{ℓ+1}.
  • Basic Subspace Iteration: A poor deterministic start matrix can be nearly orthogonal to leading singular vectors, causing a large iteration count or non-convergence.The paper analyzes these effects through relative convergence bounds.
  • Randomized Subspace Iteration: Randomized Subspace Iteration differs only by drawing Ω randomly, but this change removes bad-start behavior with near certainty and often makes q = 0 sufficient.The randomized algorithm is a special case of the basic method, so the deterministic analysis also applies.
  • SVD Truncation: SVD truncation can worsen 2-norm and Frobenius-norm accuracy and adds computation, but it can retain reliable leading singular values and a good rank-k approximation.Whether to truncate depends on the desired rank and approximation quality.

3. Setup.

This section establishes matrix-analysis tools for assessing low-rank approximations and singular-value errors. It shows that restricted SVD truncation is optimal in Frobenius norm, while 2-norm analysis is more limited.

  • Matrix-Analysis Tools: Orthogonal projection cannot increase singular values: σ_j(A) ≥ σ_j(Q^T A), and the same principle applies to submatrices.These interlacing relations provide baseline limitations for compressed approximations.
  • Approximation Bounds: The Reverse Eckart and Young theorem relates rank-k approximation error to the approximation errors in the leading singular values.When the residual parameter η is much smaller than σ_{k+1}, the resulting bounds can be substantially sharper.
  • Approximation Bounds: A small η guarantees good approximations to all k leading singular values, with η = 0 forcing equality of the leading singular values.The paper later obtains stronger singular-value results through its Section 4 analysis.
  • Restricted SVD: The rank-k truncated SVD of Q^T A is optimal for the restricted Frobenius-norm problem.This generalizes the optimality property of the truncated SVD under the constraint that the approximation lies in range(Q).
  • Restricted SVD: The corresponding 2-norm restricted problem is harder because B_k need not solve it, making the paper's 2-norm bounds potentially weaker.The analysis uses a mixed 2-norm/Frobenius-norm relation.

F + 2trace

The paper uses orthogonality of the truncated SVD components to establish a Frobenius-norm relation underlying its approximation analysis.

  • F + 2trace: The cross term A_k(A − A_k)^T vanishes, yielding relation (3.6) after combining the resulting identity with equation (3.8).This orthogonality supports the Frobenius-norm analysis.
  • F + 2trace: The deterministic convergence section derives a relative singular-value lower bound and a matrix approximation-error upper bound, then uses them for probabilistic analysis.The framework also covers Algorithm 1.1 as the q = 0 special case of randomized subspace iteration.

4. Deterministic Analysis.

The deterministic analysis identifies how spectral gaps and the randomized start influence convergence of singular values and low-rank approximations. It derives lower and upper bounds governed by the auxiliary matrix H1 and related start-matrix factors.

  • Deterministic Framework: The output QB_k is the rank-k truncated SVD of Q Q^T A because Q has orthonormal columns.Thus the algorithm's compressed SVD directly represents an orthogonal-projection approximation.
  • Convergence Factors: Choosing p controls a trade-off between oversampling for reliability and oversampling for faster convergence.The partition Σ = diag(Σ1, Σ2, Σ3) isolates the leading, intermediate, and trailing singular-value blocks.
  • Convergence Factors: A spectral gap with σ_{ℓ−p+1} much smaller than σ_k makes H1 converge rapidly as q increases, provided the randomized start factor is not too large.When p = ℓ − k, this gap disappears and the bounds align more closely with the earlier theorem.
  • Analysis Structure: The section develops deterministic lower bounds for singular values and upper bounds for rank-k approximations before the statistical analysis of Section 5.The analysis assumes full row rank for Ω̂_1 in the stated lemma and extends some relations by continuity.
  • Singular-Value Convergence: Theorem 4.3 gives lower bounds for σ_j(B_k), and σ_k(B_k) is a good approximation to σ_k when ||H1||_2 is small.Interlacing supplies the matching upper bound σ_k(B_k) ≤ σ_k.
  • Singular-Value Convergence: Fast singular-value convergence requires σ_{ℓ−p+1} ≪ σ_j, whereas a large ||Ω̂_2||_2 ||Ω̂_1^†||_2 can slow convergence despite rapid singular-value decay.Randomization is intended to keep this start-matrix factor reasonably sized with near certainty.
  • Low-Rank Approximation Error: Theorem 4.4 supplies an upper bound for the rank-k approximation error, governed by the same two key factors as singular-value convergence.The Gaussian case permits average-case analysis for all p through these bounds.

5. Statistical Analysis.

The statistical analysis establishes average-case and large-deviation bounds for randomized subspace iteration, showing that oversampling and singular-value decay strongly affect convergence and accuracy.

  • Statistical tools: Gaussian-matrix propositions provide the statistical tools used to analyze singular-value and low-rank approximation error bounds.The analysis uses rank, pseudo-inverse, Lipschitz, and tail bounds for Gaussian matrices.
  • Average-case bounds: Theorem 5.6 bounds the approximated leading singular values of a rank-k approximation computed by Algorithm 2.2.The theorem analyzes all j = 1, · · ·, k and separates cases according to the oversampling parameter p.
  • Parameter choices: Oversampling generally improves singular-value convergence, while the choice p = 0 can be slower than p > 0.The paper recommends avoiding p = 0 in general, except in settings such as a large singular-value gap or very limited oversampling.
  • Average-case bounds: For rapidly decaying singular values, Algorithm 2.2 can converge rapidly even with q = 0, making little computational work sufficient for an excellent low-rank approximation.The convergence becomes faster as the singular values decay faster.
  • Large-deviation bounds: Large-deviation guarantees allow an arbitrarily small but nonzero exception probability, although small p can require substantially more iterations.For p ≤1, the large-deviation bound is described as particularly severe; even p = 2 may require extra iterations for tiny ∆.

6. Rank-revealing Factorizations.

This section defines rank-revealing approximations and shows that Algorithm 2.2 produces them with high probability, while contrasting them with strong RRQR factorization.

  • Definition: A rank-revealing rank-k approximation has rank at most k and must preserve all k leading singular values within dimension-dependent polynomial factors.This extra singular-value condition makes it stronger than an ordinary rank-k approximation.
  • Definition: Preserving leading singular values matters because an ordinary low-rank approximation may have poor singular values even when its approximation error is small.The paper illustrates this issue using matrices whose leading k + 1 singular values are identical.
  • Randomized approximation: Algorithm 2.2 produces a rank-revealing rank-k approximation with probability at least 1−∆.The section compares this randomized approximation with the strong RRQR factorization of Gu and Eisenstat.
  • Comparison with RRQR: The randomized approximation satisfies properties compatible with the strong RRQR inequalities, but it does not include the permutation selecting k linearly independent columns.The missing permutation means it does not provide the same column-selection information, although matrix-matrix operations may make it cheaper.

7. Condition Number Estimation.

The paper applies randomized rank-1 approximation to matrix norm estimation, obtaining a reliable 2-norm estimate and improving initialization for a randomized version of Hager’s method.

  • Motivation: Condition estimation typically requires separate estimates of ∥A∥ and ∥A−1∥, ideally without direct entry access and at lower cost than factorization or inversion.Matrix norm estimators access matrices through matrix-matrix or matrix-vector products.
  • Randomized Hager’s Method: Randomized Hager’s Method computes a rank-1 approximation, initializes Hager’s method with its right singular vector, and returns a 1-norm estimate.The iteration count can be restricted to as few as 1 or 2 because local maximization is unnecessary for this use.
  • Guarantee: Algorithm 7.2 provides a reliable ∥A∥1 estimate with exception probability at most ∆.The guarantee is stated in Corollary 7.1 for any 0 < ∆≪1.
  • Numerical evidence: The experiments report that ℓ = 5 worked well for matrix norm estimation and that randomized methods can require different effort for different singular-value distributions.The paper frames the experiments as support for its probabilistic analysis.
  • Numerical example: In the numerical example, Algorithm 7.2 estimated ∥A∥1 as approximately 2.46×10^11, whereas Hager’s method returned approximately 4.72×10^1 for a true value near 8.35×10^11.The example demonstrates robustness against a starting vector that causes the classical method to fail badly.

8. Numerical Experiments.

The experiments evaluate improved randomized subspace methods across norm estimation, low-rank approximation, structured matrix preconditioning, and face recognition. They show strong performance for rapidly decaying singular values, with application-specific trade-offs involving oversampling and approximation accuracy.

  • 8.1. Improved Randomized Power Iteration.: Algorithm 8.1 combines a rank-ℓ2 approximation from Algorithm 1.1 with an approximate right singular vector start for Algorithm 2.1.The procedure first uses ℓ1>ℓ2≥k, then runs subspace iteration with the improved start matrix.
  • 8.1. Improved Randomized Power Iteration.: Algorithm 8.1 reached ∥A∥2 in about 250 steps versus about 350 for Algorithm A.3 when σ2/σ1 was large.For a small σ2/σ1 ratio, both algorithms performed similarly, while the extra cost of computing the initial Ω offset Algorithm 8.1’s modest speed advantage.
  • 8.2. low-rank approximation.: For matrices with relatively rapid singular-value decay, Algorithm 1.1 outperformed Algorithm 2.2 for q>0 and also outperformed MATLAB’s svds.The comparison used k=50 and counted matrix-vector multiplies; the paper attributes expected advantages partly to matrix-matrix products on modern architectures.
  • 8.3. Structured Matrix Computations.: All tested oversampling choices substantially reduced PCG iterations, but larger p values usually produced only small additional reductions.Because higher p increases HSS preconditioner construction cost, the experiment found small p typically sufficient; p between 10 and 20 is cited as a rule of thumb.
  • 8.4. Eigenfaces.: The face-recognition experiment used 200 training images from 40 subjects and the remaining images for classification.The images were 92×112 grayscale pixels, with five training images per subject.
  • 8.4. Eigenfaces.: In face recognition, p=40 produced results very similar to truncated SVD despite some singular values being accurate to only 1 to 2 digits.Smaller p values performed worse than truncated SVD, demonstrating that constant-factor or limited singular-value accuracy can suffice for this application.

9. Conclusions and Future Work.

The paper concludes that randomized subspace-iteration algorithms can provide accurate low-rank approximations, singular values, rank-revealing factorizations, and reliable condition-number estimates. It also identifies convergence analysis for singular vectors and parameter selection as future work.

  • 9. Conclusions and Future Work.: Future work includes convergence analysis for singular vectors and extensions to subspace iteration methods for symmetric and nonsymmetric eigenvalue problems.The paper expects singular-vector results to be compatible with its singular-value analysis.
  • 9. Conclusions and Future Work.: Choosing q and ℓ to minimize the overall cost for a fixed rank k remains an open issue.The conclusion identifies parameter selection as an immediate direction for future research.
  • 9. Conclusions and Future Work.: The resulting low-rank approximations are suggested to be rank-revealing factorizations and can reliably estimate condition numbers.The rank-1 case is connected to estimating matrix 2-norms.

Appendix S2. Further Convergence Considerations.

This appendix examines how singular-value decay affects randomized subspace-iteration parameter choices and convergence rates. It focuses on matrix-vector-product costs and identifies the expression controlling singular-value convergence.

  • Appendix S2. Further Convergence Considerations.: The appendix identifies a singular-value distribution where Algorithm 1.1 is likely to outperform Algorithm 2.2 and another where choosing ℓ can speed convergence.The convergence analysis is organized around how singular-value decay influences parameter choices.
  • Appendix S2. Further Convergence Considerations.: The analysis counts only Au and A^T v products, each costing O(mn) flops, because truncated SVD costs are smaller when ℓ≪n.Under this accounting, total matrix-vector-product cost is O(mnℓ).
  • Appendix S2. Further Convergence Considerations.: The expression analyzed in this section is treated as the key factor controlling singular-value convergence.The supplied passage does not reproduce the expression itself.

S2.1. Rapidly decaying singular value distributions.

For exponentially or faster decaying singular values, the appendix compares Algorithm 1.1 and subspace iteration at roughly equal matrix-vector-product cost. The comparison suggests little convergence advantage for Algorithm 2.2 over Algorithm 1.1.

  • S2.1. Rapidly decaying singular value distributions.: The decay model applies when singular values, except possibly the first few, decrease exponentially or faster.It is parameterized by a constant α>0 and applies for s,t>1.
  • S2.1. Rapidly decaying singular value distributions.: The comparison equalizes work by allowing both algorithms the same number of matrix-vector products, bℓ.Algorithm 2.2 uses bℓ=(2q+1)ℓ matrix-vector products in the stated setup.
  • S2.1. Rapidly decaying singular value distributions.: Algorithm 1.1 has a tighter upper bound on the compared convergence ratio than Algorithm 2.2 under the rapid-decay model.The appendix uses this bound comparison to assess singular-value convergence at comparable cost.
  • S2.1. Rapidly decaying singular value distributions.: The comparison suggests little convergence advantage for Algorithm 2.2 over Algorithm 1.1 when singular values decay exponentially or faster.This is the section’s stated conclusion for the rapid-decay model.

S2.2. Slowly decaying singular value distributions.

For slowly decaying singular values, the optimal sampling parameter ℓ_opt is characterized analytically and yields exponential convergence. A naive parameter choice is slower, while automatically selecting parameters remains an open challenge.

  • S2.2. Slowly decaying singular value distributions.: The optimal ℓ, denoted ℓ_opt, is the unique solution of g(ℓ) = 0.The supplied passage states that ℓ_opt satisfies this defining equation.
  • S2.2. Slowly decaying singular value distributions.: Choosing ℓ = ℓ_opt makes Algorithm 2.2 converge exponentially when bℓ is much larger than ℓ_opt.The resulting convergence rate is described as exponential under the stated sampling regime.
  • S2.2. Slowly decaying singular value distributions.: The optimized convergence rate is e−O(bℓ/k), whereas the naive choice ℓ = k + p yields e−O(bℓ/k^2).The naive choice is slower by a factor of O(k) in the exponent.
  • S2.2. Slowly decaying singular value distributions.: The best choice ℓ_opt can significantly accelerate the overall convergence rate.This is the paper’s stated consequence of optimizing ℓ.
  • S2.2. Slowly decaying singular value distributions.: Designing an efficient algorithm that automatically selects parameters for different singular value distributions remains necessary.The section identifies automatic parameter selection as an unresolved research need.

S2.3. Adaptive Randomized Algorithms.

The section presents an adaptive randomized subspace iteration method for computing rank-k approximations to a tolerance, then evaluates it on latent semantic indexing data. The experiments report greater agreement with the truncated SVD when randomized approximations are more accurate, with q = 2 outperforming q = 0 in this example.

  • S2.3. Adaptive Randomized Algorithms.: The adaptive method targets rank-k approximations for a supplied accuracy tolerance τ and failure tolerance Δ.Its inputs include matrix dimensions, tolerances, sampling and iteration parameters, and a sampling limit.
  • S2.3. Adaptive Randomized Algorithms.: Algorithm S2.1 draws a random test matrix, forms Y, orthogonalizes its columns, compresses A to B = Q^T A, and returns QB_k.The core procedure computes a truncated SVD of the compressed matrix B.
  • S2.3. Adaptive Randomized Algorithms.: The adaptive algorithm may add sampled columns iteratively and stops when the sampling size exceeds the limit for the given tolerance.The implementation introduces parameters c and b to control initial convergence and minimum columns sampled per iteration.
  • S2.3. Adaptive Randomized Algorithms.: In the LSI experiment, the TDT2 subset contains 9,394 documents and over 36,000 terms across 96 semantic categories.The corpus was collected from six news, radio, and television sources.
  • S2.3. Adaptive Randomized Algorithms.: More accurate randomized approximations produce more agreement with the truncated SVD in query matches, and q = 2 performs significantly better than q = 0.The authors caution that agreement does not always mean a better match, but suggest a small positive q for this LSI example.
Loading 1408.2208v1…