Source-linked AI summary
An algorithm for the principal component analysis of large data sets
Nathan Halko, Per-Gunnar Martinsson, Yoel Shkolnisky, Mark Tygert
TL;DR
Very large data sets may not fit in RAM, limiting direct PCA computation. This paper adapts a randomized block Lanczos method for out-of-core PCA, and reports nearly optimal accuracy, including on disk-resident data sets with less than one hundredth fitting in RAM.
Problem
The paper addresses how to compute PCA for data sets too large to be stored in a typical computer’s RAM.
Method
The paper adapts a randomized block Lanczos method using a numerically stable rank-k approximation procedure for out-of-core data.
Results
The algorithm produces nearly optimal spectral-norm accuracy with high probability, and all reported disk-resident examples succeeded when less than a hundredth of each data set fit in RAM.
Takeaways & Limitations
The procedure enables PCA of very large data sets stored on disk, with core steps that parallelize easily.
Takeaways & Limitations
The algorithm can underflow or overflow when floating-point exponent range cannot represent both ||A||2 and ||A||2i+1.
Abstract
from arXiv · showhide
Recently popularized randomized methods for principal component analysis (PCA) efficiently and reliably produce nearly optimal accuracy --- even on parallel processors --- unlike the classical (deterministic) alternatives. We adapt one of these randomized methods for use with data sets that are too large to be stored in random-access memory (RAM). (The traditional terminology is that our procedure works efficiently "out-of-core.") We illustrate the performance of the algorithm via several numerical examples. For example, we report on the PCA of a data set stored on disk that is so large that less than a hundredth of it can fit in our computer's RAM.
1. Introduction.
The paper adapts a randomized block Lanczos method to compute PCA for data sets too large for RAM, targeting nearly optimal accuracy while reducing data-access costs.
- PCA is widely used in machine learning, statistics, data analysis, data mining, and information retrieval.
- The paper computes PCA through a randomized version of the block Lanczos method.
- The method typically needs only a couple of iterations to achieve nearly optimal accuracy with overwhelmingly high probability.The stated probability is typically 1 −10−15 or greater and is independent of the data analyzed.
- The algorithm is adapted for data sets that cannot be stored in a typical computer’s random-access memory.
- The procedure minimizes how many times each matrix entry is accessed.A related strategy instead minimizes disk seeks and maximizes approximation dimensions for a given RAM size.
2. Informal description of the algorithm.
The algorithm first uses randomized matrix products and orthogonalization to identify most of a matrix’s range, then applies compact linear algebra to obtain a low-rank approximation.
- The method constructs a low-rank approximation UΣV⊤ to a real matrix A.
- Random Gaussian vectors are multiplied by A and repeated products involving A(A⊤A)^i suppress smaller singular values.
- A pivoted QR decomposition forms Q with orthonormal columns that approximate most of the range of A.
- The algorithm’s efficiency comes from using an approximation basis with substantially fewer columns than A.
- An SVD of a reduced matrix is used to recover U, Σ, and V, whose leading blocks form the final approximation.
3. Summary of the algorithm.
The paper specifies a numerically stable randomized rank-k approximation whose spectral-norm error is nearly optimal with high probability, while identifying parameter choices and numerical caveats.
- The algorithm constructs a rank-k approximation UΣV⊤ to a real matrix A using randomized range finding and subsequent linear-algebraic manipulations.
- For most applications, i = 1 or i = 2 is sufficient, while increasing i can improve accuracy but reduces efficiency.
- The approximation achieves nearly optimal spectral-norm accuracy with high probability independent of A.The benchmark is the best rank-k error σk+1; failure probability is small for C = 10 and negligible for C = 100.
- The implementation uses a Gaussian test matrix, pivoted QR decomposition, reduced products, an SVD, and extraction of the leading k blocks.
- FFT variants offer negligible acceleration for out-of-core computations, so the paper favors simpler techniques.
- The method can underflow or overflow when floating-point exponent range cannot represent both ||A||2 and ||A||2i+1.Processing A/||A||2 or using the alternative algorithm in [8] are stated remedies.
4. Out-of-core computations.
The out-of-core implementation keeps smaller intermediate matrices in RAM while handling the full data matrix through entry-wise computation or disk streaming.
- Steps 2, 4, 5, and 6 require only O((i+1)l(m+n)) stored entries and can run entirely in RAM.
- Steps 1 and 3 involve the full mn-entry matrix A and therefore require specialized out-of-core handling.
- When entries, rows, or columns can be evaluated individually, the algorithm computes the needed products immediately using only RAM.
- When A is stored on disk in row-major format, the implementation retrieves batches of rows, forms products with G, stores results in H(0), and repeats.
5. Computational costs.
The paper quantifies computational, memory, disk-access, and data-transfer costs for its out-of-core PCA implementations. For most applications, small iteration counts and l approximately equal to k keep the randomized method practical for large matrices.
- For most applications, i ≤ 2 suffices, whereas classical Lanczos generally requires many iterations, making out-of-core computation prohibitive.The paper identifies iteration count as a central cost distinction between the randomized and classical methods.
- Costs with on-the-fly evaluation of matrix entries: On-the-fly evaluation costs 2(i + 1)C_A + O(ilN_A + i^2l^2m) flops.C_A measures the cost of evaluating all nonzero entries, and N_A is the number of nonzero entries.
- Costs with storage on disk: With A stored on disk, forming the main matrix products requires O(lmn) flops, O(mn/j) disk accesses or seeks, and O(mn) transferred words.Here j denotes the number of floating-point words of RAM available to the algorithm.
- Costs with storage on disk: The full disk-based algorithm transfers approximately 2(i + 1)mn floating-point words.This transfer estimate is reported after summing the costs of the six algorithmic steps.
6. Numerical examples.
Numerical tests evaluate the algorithm on synthetic and face-image data, using both on-disk storage and on-the-fly generation. The experiments report accuracy behavior, runtime scaling, and a 150 GB PCA computed with only 1.5 GB of RAM.
- Experimental setup: The experiments used a 1.5 GB RAM laptop connected to an external hard drive, with examples storing matrices on disk or generating them on-the-fly.The processor was a single-core 32-bit 2-GHz Intel Pentium M, and experiments ran in Matlab 7.4.0.
- Experimental setup: Pseudorandom-number quality and distribution had little effect on the algorithm’s accuracy across the numerical experiments.The reported results used Matlab’s built-in pseudorandom number generator.
- Synthetic data: The synthetic-data tests found correlations with the leading singular vectors tending to 1 as m increased, while on-the-fly runtime was roughly proportional to m.The correlations corresponded to the three greatest singular values, and the runtime behavior agreed with the paper’s running-time model.
- Measured data: For face images, the algorithm computed a rank-50 PCA of a 393,216 × 102,042 matrix containing 150 GB of data in 12.3 hours.The data consisted of duplicated and partially randomized FERET face images, processed from disk on the 1.5 GB RAM laptop.
- Measured data: The face-image experiment produced computed singular values and eigenfaces corresponding to the five greatest singular values, but brute-force low-rank approximation alone is inadequate for face recognition.The authors present the computation as useful for processing or preprocessing very large data sets, while specialized image-processing methods remain needed.
7. An application.
The paper applies its out-of-core PCA algorithm to millions of noisy cryo-electron microscopy images of the E. coli 50S ribosomal subunit. The computed singular vectors support satisfactory denoising by projection onto the leading 150 vectors.
- Data set: The data set contains 3,600,000 images, each 129 pixels wide and high, formed from noisy projections of the E. coli 50S ribosomal subunit.The images represent uniformly random orientations, white Gaussian noise 32 times larger than the original images, and rotations from 0 through 359 degrees.
- Results: The algorithm computed 250 singular values and 25 corresponding right singular vectors, while denoising used projections onto the vectors associated with the 150 greatest computed singular values.The resulting denoised images were described as clearly satisfactory.
8. Conclusion.
The paper develops PCA techniques for data sets too large for RAM and evaluates them on synthetic, simulated, and physical data. The scheme succeeded on data sets for which less than a hundredth fit in RAM, while achieving nearly optimal spectral-norm accuracy and offering parallelizable core steps.
- Contribution: The paper describes PCA techniques for data sets too large to be stored in random-access memory and illustrates them on standard tests, simulations, and physical measurements.
- Results: Less than a hundredth of several on-disk data sets could fit in the computer’s RAM, yet the scheme always succeeded.
- Accuracy and deployment: Theorems, proofs, and numerical validations demonstrate nearly optimal spectral-norm accuracy, with similar results available for the Frobenius/Hilbert-Schmidt norm.The core procedure steps also parallelize easily.
Appendix.
The appendix describes a randomized power-method procedure for estimating a matrix’s spectral norm without explicitly forming D⊤D, and reports its accuracy in numerical examples.
- The procedure estimates the spectral norm of D to check whether a matrix approximation is accurate.It uses D as the difference between the original matrix and its approximation.
- The estimator applies D and D⊤ successively to vectors rather than forming D⊤D explicitly.
- The estimate p_j,k(D) is produced by j power-method steps initialized with k independent Gaussian random vectors.The vectors have independent, identically distributed entries with zero mean and unit variance.
- The estimate approaches the true spectral norm rapidly as j increases and is accurate within a factor of two with very high probability for fairly small j.
- The numerical examples use j = 6 and set k equal to the rank of the constructed approximation.The procedure is applied with D = A − UΣV⊤.