Source-linked AI summary
A randomized algorithm for principal component analysis
Vladimir Rokhlin, Arthur Szlam, Mark Tygert
TL;DR
PCA requires low-rank approximation, but existing efficient methods lack good-accuracy guarantees for large matrices in challenging spectral regimes. This paper presents an efficient algorithm for arbitrary-size matrices and shows theoretically and numerically that its accuracy is very close to the best possible, with high probability and stated computational cost.
Problem
PCA reduces to low-rank approximation, yet efficient algorithms lack good-accuracy guarantees when matrices are large and their spectra decay slowly.
Method
The paper describes an efficient randomized algorithm that produces a rank-k approximation to the data matrix for PCA.
Results
Accuracy is very close to the best possible with probability typically 1 −10−15, and the algorithm requires O(nmki) floating-point operations when i > 0.
Takeaways & Limitations
The algorithm provides a practical low-rank approximation approach for PCA on matrices of arbitrary sizes, supported by theoretical analysis and numerical examples.
Takeaways & Limitations
Dense pivoted QR comparisons were limited to m ≤4096 by the available 1 GB of random-access memory.
Abstract
from arXiv · showhide
Principal component analysis (PCA) requires the computation of a low-rank approximation to a matrix containing the data being analyzed. In many applications of PCA, the best possible accuracy of any rank-deficient approximation is at most a few digits (measured in the spectral norm, relative to the spectral norm of the matrix being approximated). In such circumstances, efficient algorithms have not come with guarantees of good accuracy, unless one or both dimensions of the matrix being approximated are small. We describe an efficient algorithm for the low-rank approximation of matrices that produces accuracy very close to the best possible, for matrices of arbitrary sizes. We illustrate our theoretical results via several numerical examples.
1. Introduction.
The paper targets efficient rank-k approximation for PCA, where standard methods can be inaccurate on large matrices with slowly decaying singular values. Its randomized algorithm achieves accuracy close to the best possible with high probability.
- PCA reduces computationally to finding a low-rank approximation of the data matrix.
- O(nm^2) flops may be required by standard SVD algorithms, while pivoted QR uses O(nmk) flops but can lose accuracy.
- When m ≥10,000 and σ1/σk+1 ≤100, pivoted QR typically produces ∥A−B∥∼∥A∥ when noise causes slow singular-value decay.
- The proposed algorithm produces a rank-k approximation whose spectral-norm error is bounded relative to σk+1 with very high probability.
- The algorithm requires O(nmki) operations for i > 0 and O(nmk) flops when i = 0.
- Numerical evidence indicates higher accuracy than classical pivoted QR for matrices with slowly decaying spectra.
2. Preliminaries.
This section collects singular-value facts, approximation results, Gaussian random-matrix bounds, and a monotonicity lemma used later in the analysis.
- Singular values of general matrices: A rectangular block cannot have a greatest singular value exceeding that of the full matrix.
- Singular values of general matrices: An n × l matrix R admits an approximation QS with Q having orthonormal columns and S having k rows.
- Singular values of random matrices: For Gaussian random matrices, the lemmas provide high-probability upper bounds on greatest singular values and lower bounds on least singular values.
- A monotone function: The function f defined on (0, ∞) decreases monotonically for x > α.
3. Mathematical apparatus.
The mathematical apparatus connects randomized Gaussian sketches and power iterations to approximation error bounds for the proposed low-rank algorithm.
- The section develops lemmas used to bound the algorithm’s accuracy.
- If Q has orthonormal columns and QS approximates (G(AAT)^iA)T, then AQQT approximates A under a two-term error bound.
- With very high probability, a matrix F exists with reasonably small norm such that FG(AAT)^iA approximates A.
- Gaussian-sketch lemmas bound singular values of GA and G(AAT)^iA in terms of singular values of A.
4. Description of the algorithm.
The paper develops randomized low-rank approximation algorithms that recover an accurate approximation to A with high probability while avoiding full SVD costs. The basic construction uses random sampling, orthonormalization, a small SVD, and optional power iterations or a Blanczos variant to improve accuracy.
- 4.1. The algorithm.: With very high probability, the resulting factors produce a rank-k approximation whose error is bounded by a constant times σk+1.U and V have orthonormal columns, Σ is nonnegative diagonal, and C is independent of A but depends on algorithm parameters.
- 4.1. The algorithm.: Applying A^T(A A^T)^i to random vectors enhances singular-value decay before linear-algebraic manipulations recover an approximation to A.This power-iteration idea identifies most of the relevant range of A^T and motivates the algorithm’s accuracy.
- 4.1. The algorithm.: The algorithm constructs an approximate SVD of A using randomized matrix products, orthonormal bases, and an SVD of a smaller matrix.The five-step procedure forms a Gaussian test matrix, constructs Q, computes a product involving A, and factors the resulting matrix.
- 4.1. The algorithm.: The proof bounds the projection error by combining two terms involving randomized range capture and the residual after applying a power iteration.The analysis uses the decomposition in (4.17), selects an index based on adjacent singular values, and combines probability bounds Φ and Ψ into Π.
- 4.1. The algorithm.: Choosing l = k + 12, β = 2.57, and γ = 2.43 gives failure probability below 10^-15, although experiments indicate the factor 100l is overly conservative.A stronger bound is available when l ≥ 2k, reducing the right-hand side of the accuracy bound by an additional factor.
- 4.2. Computational costs.: The basic algorithm costs CPCA = (il + k) · CA + (il + l) · CAT + O(l^2 n) floating-point operations and needs only operator applications.It applies A to il + k vectors and A^T to il + l vectors without requiring explicit access to individual matrix entries.
5. Numerical results.
Numerical experiments compare the proposed randomized algorithms with varying iteration counts and with classical pivoted QR, including matrices whose singular values decay slowly. The results indicate improved accuracy for positive iteration counts, favorable scaling, and greater roundoff tolerance for the more expensive blanczos variant.
- The experiments evaluate the proposed five-step algorithm, blanczos, and pivoted QR on matrices with varying sizes and best rank-k approximation accuracies.Tables 1–5 cover the proposed algorithms and blanczos; Table 6 covers pivoted QR-decompositions.
- The algorithm performs much better with i > 0, and its accuracies in Table 1 exceed those in Table 2.Table 1 uses i = 1, whereas Table 2 uses i = 0.
- The proposed algorithm produces higher accuracy than classical pivoted QR-decompositions for matrices whose spectra decay slowly.The comparison is between the accuracies in Tables 1 and 6 for the tested slowly decaying spectra.
- The observed accuracies appear proportional to m^1/(4i+2) σk+1 for Subsection 4.1 and m^1/(4i) σk+1 for Subsection 4.3.These empirical scalings agree with the corresponding theoretical bounds cited by the authors.
- The theoretical accuracy factor 100l appears much greater than necessary in numerical experiments.The authors report this conclusion alongside the high-probability accuracy guarantee.
- Blanczos has accuracy similar to Subsection 4.1 with matching or reduced i, depending on the relationship between (σk)^2i+1 and machine precision.The authors suspect blanczos is a better general-purpose black-box PCA tool because it tolerates roundoff, despite somewhat higher cost.
6. Appendix.
The appendix restates and proves lemmas underlying the randomized low-rank approximation, bounding projection error through intermediate approximations and Gaussian constructions.
- The appendix restates and proves Lemmas 3.1 and 3.2 from Section 3.
- Lemma 6.1 bounds ∥A Q Q^T − A∥ using the approximation errors of F G B and Q S.The bound combines the error in approximating A by F G B with the error in approximating (G (A A^T)^i A)^T by Q S.
- The appendix derives the error bound by applying the triangle inequality, orthonormality of Q, and norm inequalities before combining the resulting estimates.
- The proof constructs F so that F G (A A^T)^i A approximates A with very high probability when G has independent standard Gaussian entries.
- The proof uses an SVD of A and partitions G U and Σ into blocks to construct auxiliary matrices and F.H, R, and Γ are blocks of G U; S, T, and Θ are blocks of Σ, while P and F are then defined from these quantities.