PaperlayerView exact revision on arXiv

Source-linked AI summary

Randomized Matrix Decompositions using R

N. Benjamin Erichson, Sergey Voronin, Steven L. Brunton, J. Nathan Kutz

arXiv:1608.02148v5stat.COcs.MSstat.ME

TL;DR

Massive datasets make traditional low-rank matrix decompositions computationally challenging. The paper presents the R package rsvd, which uses randomized reductions to support several approximate decompositions. The reported computational advantage is substantial when the target rank is relatively small and the matrix shape is favorable, with accuracy controlled through oversampling and power iterations.

  • Problem

    Massive datasets make traditional low-rank matrix algorithms difficult to run, while matrix decompositions remain central to dimensionality reduction, data processing, and compression.

  • Method

    The paper presents rsvd, an R package implementing randomized SVD, PCA, robust PCA, interpolative decomposition, and CUR decomposition.

  • Results

    Randomized SVD can provide substantial computational advantages over other R SVD routines when the target rank k is relatively small, with approximation quality controlled by oversampling and power iterations.

  • Takeaways & Limitations

    Randomized methods provide scalable, tractable low-rank computations, trading approximation accuracy against computational cost.

  • Takeaways & Limitations

    Significant speedups depend on matrix shape and target rank, are less impressive for tall and thin matrices, and the input matrix must fit into fast memory.

Abstract

from arXiv · show

Matrix decompositions are fundamental tools in the area of applied mathematics, statistical computing, and machine learning. In particular, low-rank matrix decompositions are vital, and widely used for data analysis, dimensionality reduction, and data compression. Massive datasets, however, pose a computational challenge for traditional algorithms, placing significant constraints on both memory and processing power. Recently, the powerful concept of randomness has been introduced as a strategy to ease the computational load. The essential idea of probabilistic algorithms is to employ some amount of randomness in order to derive a smaller matrix from a high-dimensional data matrix. The smaller matrix is then used to compute the desired low-rank approximation. Such algorithms are shown to be computationally efficient for approximating matrices with low-rank structure. We present the \proglang{R} package rsvd, and provide a tutorial introduction to randomized matrix decompositions. Specifically, randomized routines for the singular value decomposition, (robust) principal component analysis, interpolative decomposition, and CUR decomposition are discussed. Several examples demonstrate the routines, and show the computational advantage over other methods implemented in R.

1. Introduction

Large datasets make traditional matrix decompositions increasingly difficult to compute, motivating randomized methods that reduce high-dimensional matrices before approximation. The rsvd package brings randomized SVD, PCA, RPCA, interpolative, and CUR decompositions to R, with substantial speedups when matrix shape and target rank are favorable.

  • Motivation: Large-scale datasets challenge traditional matrix algorithms because data growth is outpacing computational performance.SVD is especially computationally expensive despite its widespread use in dimensionality reduction, processing, and compression.
  • Randomized framework: Randomized matrix algorithms derive smaller matrices that preserve essential information before computing approximate decompositions.This probabilistic framework has been used for SVD, PCA, pivoted QR and LU, dynamic mode decomposition, and related tasks.
  • Package contribution: The rsvd package provides randomized routines for SVD, PCA, robust PCA, interpolative decomposition, and CUR decomposition.The package is available through CRAN and GitHub for use in R.
  • Performance and scope: Randomized SVD produces significant speedups over existing R SVD routines when computing dominant singular values and vectors, particularly for rapidly decaying spectra.Its computational benefits also transfer to PCA and can accelerate RPCA, ID, and CUR decomposition.

2. Probabilistic framework for low-rank approximations

The probabilistic framework replaces direct low-rank decomposition of a massive matrix with randomized subspace construction followed by deterministic compression. Oversampling and power iterations improve the sampled basis when singular values decay slowly, trading accuracy against additional computation.

  • Framework: Low-rank approximation factors a large matrix A into smaller matrices E and F, reducing storage from mn words to mr + nr words.The columns of E span A’s column space, while the rows of F span its row space.
  • Two-stage framework: The framework constructs an orthonormal basis Q so that A ≈ QQ⊤A, then forms B := Q⊤A and decomposes the smaller matrix.Randomness is used in constructing Q; the subsequent compression stage is deterministic.
  • Random projections: Random projections form Y from independent randomly weighted combinations of A’s columns, reducing the number of columns from n to k while approximately preserving Euclidean distances.The projections use yi := Aωi, with random vectors typically drawn from a standard normal distribution.
  • Oversampling: Oversampling uses l := k + p random projections, with p = 5 or 10 often sufficient to obtain a basis comparable to the best possible basis.Additional projections reduce variation in the random test matrix and improve sketch quality.
  • Power iteration scheme: Power iterations accelerate singular-value decay, helping the sketch emphasize dominant information and suppress unwanted information.For slowly decaying spectra, q = 1, 2, or 3 iterations can considerably improve approximation accuracy, but require q additional passes over the input matrix.

3. Randomized singular value decompositions

The SVD provides stable low-rank representations and supports core computational tasks, but computing a truncated SVD directly becomes costly for massive matrices. Its economy and reduced forms retain only the components needed for the target rank.

  • Role of the SVD: The SVD supports low-rank approximation, pseudoinverse computation, least-squares solutions, and machine-learning methods including PCA and robust PCA.It is described as a numerically stable matrix decomposition and a workhorse algorithm for these applications.
  • Economy SVD: The economy SVD computes only the singular vectors and values associated with the n right singular vectors and is more compact when n ≪ m.The economy form is the default form of R’s base svd() function.
  • Reduced SVD: A reduced rank-k SVD represents A approximately as a sum of k rank-one matrices, with k chosen according to reconstruction accuracy or desired dimensionality.Choosing k near the effective rank favors reconstruction, while a much smaller k favors a low-dimensional representation.
  • Optimality: The Eckart–Young theorem identifies the truncated SVD as the optimal rank-k reconstruction in the least-square sense.The corresponding reconstruction error is reported in both spectral and Frobenius norms.
  • Computational cost: Computing the full SVD of an m × n matrix costs O(mn^2), after which the first k components can be extracted.This cost makes direct truncated SVD computation challenging for massive datasets.

3.3. Randomized algorithm

The randomized SVD embeds the decomposition in a probabilistic low-dimensional projection framework, computing a basis Q, decomposing a smaller matrix B, and reconstructing the requested factors. The approach offers lower communication costs, parallelism, and broad applicability.

  • Algorithmic rationale: Randomized algorithms compute an approximate rank-k SVD at cost O(mnk), substantially below full-SVD truncation when matrix dimensions are large.The algorithms target low-rank structure and are controlled by the target rank k.
  • Pipeline: The randomized SVD first constructs Q ∈ R^{m×l}, projects A to B := Q⊤A ∈ R^{l×n}, computes B’s deterministic SVD, and recovers the left singular vectors using Q.Oversampling l = k + p and power iterations q can improve the basis before decomposing B.
  • Implementation: The implementation exposes oversampling and power-iteration parameters, then returns the first k left singular vectors, singular values, and right singular vectors.The target rank must satisfy k < min{m, n}; good performance is expected when k is much smaller than the ambient dimensions.
  • Practical advantages: The randomized SVD requires few, at least two, passes over the input matrix, reducing communication costs that can dominate big-data computation.Sequential reads from storage into fast memory may cost more than theoretical arithmetic complexity suggests.
  • Practical advantages: Randomized algorithms are highly parallelizable because their expensive steps are matrix–matrix operations suited to multithreading and distributed computing.Their approximation applies to matrices with arbitrary singular-value decay, approaching optimal truncated SVD performance for rapid decay with high probability.

3.4. Theoretical performance

The randomized approximation has an expected error bound governed by oversampling and power iterations. Increasing oversampling moves the bound toward the optimal σk+1(A), while power iterations improve performance by accelerating singular-value decay.

  • Expected error: The randomized low-rank approximation has an expected error bound relative to the optimal truncated SVD.The expectation is taken over a Gaussian test matrix, with oversampling parameter p ≥ 2 assumed.
  • Oversampling: Increasing the oversampling parameter p drives the bound’s second and third terms toward zero, approaching the theoretically optimal value σk+1(A).The result connects larger oversampling with tighter approximation-error control.
  • Power iterations: Power iterations accelerate singular-value decay while preserving eigenvectors, yielding better performance for matrices with slowly decaying spectra.The power parameter q therefore provides a second mechanism for controlling approximation quality.

3.5. Existing functionality for SVD in R

R provides accurate full SVD routines and several partial SVD implementations. Partial methods reduce computation by targeting dominant components, but their performance depends on matrix structure and algorithmic settings.

  • The base svd() function interfaces with LAPACK routines known for numerical stability and full double-precision accuracy.
  • Computing a truncated SVD by first forming the full SVD can require tremendous time for large-scale data.
  • Krylov subspace methods efficiently approximate dominant singular vectors and values using operations such as matrix-vector multiplication.
  • Approximating left singular vectors through A V_k Σ^-1 is possible, but eigendecomposition of A^⊤A squares the condition number.
  • R packages including svd, RSpectra, and irlba provide competitive partial SVD routines based on PROPACK, ARPACK-inspired methods, and implicitly restarted Lanczos bidiagonalization.

3.6. The rsvd() function

The rsvd() function provides an R-compatible randomized low-rank SVD interface. Its target rank, oversampling, power iterations, and random test-matrix distribution expose the main speed–accuracy controls.

  • The rsvd package provides an efficient low-rank SVD routine with an interface similar to base svd().
  • The input matrix A and target rank k are mandatory, with speedups expected when k < min{m, n}/4.
  • Oversampling p and power iterations q control approximation accuracy, with defaults p = 10 and q = 2.
  • The function supports normal, uniform, and Rademacher random test matrices through the sdist argument.
  • The returned model contains singular values d and approximate left and right singular-vector matrices u and v.

3.7. SVD example: Image compression

The image-compression example uses dominant singular components to reconstruct a grayscale image. Randomized SVD achieves reconstruction close to truncated deterministic SVD while reducing runtime, especially when power iterations are retained.

  • A 1600 × 1200 grayscale image is represented using a compact low-rank basis for compression.
  • The truncated SVD retains the dominant k = 100 singular values and vectors to reconstruct the image.
  • 12.1% reconstruction error is achieved using only k = 100 singular values and vectors.
  • The randomized reconstruction error is about 0.122, close to the optimal truncated SVD.
  • Randomized SVD with power iterations is nearly as accurate visually as deterministic SVD, whereas omitting subspace iterations reduces quality.
  • Across 20 runs, randomized routines attain substantial speedups while maintaining similar reconstruction errors when q ≥1.

3.8. Computational performance

The performance evaluation compares randomized and existing R SVD routines on low-rank and sparse matrices. Randomized SVD gains become more pronounced with increasing dimensions, while sparse-specialized methods perform better on sparse inputs.

  • The benchmark uses Microsoft R Open with multithreaded LAPACK on a quad-core Intel Core i7-7700K system.
  • The evaluation benchmarks SVD routines on generated low-rank matrices with varying dimensions and intrinsic rank r = 200.
  • For target rank k = 20, randomized SVD gains become pronounced as matrix dimensions increase, after RSpectra and irlba perform best at small sizes.
  • The rsvd() function achieves substantial speedups, while additional power iterations improve approximation accuracy and expose a time–accuracy trade-off.
  • Sparse-matrix benchmarks use matrices with about 5% nonzero elements.
  • RSpectra and propack show considerably better performance on sparse and structured matrices because they are specifically designed for those inputs.

4. Randomized principal component analysis

Randomized PCA uses randomized SVD to efficiently approximate dominant principal components when large data matrices have low-rank structure. The section explains PCA’s variance-based formulation, eigenvector interpretation, visualization, and examples showing comparable accuracy with faster computation.

  • Conceptual overview: PCA constructs uncorrelated principal components that sequentially capture the greatest possible remaining variation in the data.The first component captures most variation, followed by subsequent components capturing the remaining variation.
  • Conceptual overview: PCA can be formulated equivalently as variance maximization or least-squares minimization of unexplained variation.The equivalence follows because total variation equals explained plus unexplained variation.
  • PCA formulation: The principal directions are orthonormal eigenvectors of the covariance or correlation matrix, with eigenvalues measuring explained variation.The eigenvectors form a rotation matrix and project observations into principal-component coordinates.
  • Choosing the target rank: Retaining k components is application-dependent: too many adds little, while too few produce a poor approximation.The eigenvalues quantify the variance captured by the retained components.
  • Randomized PCA: Randomized PCA uses randomized SVD to efficiently approximate the dominant k principal components, implemented by rpca().The approach is especially relevant when large-scale data can be approximated by a few principal components.
  • Examples: The handwritten-digit example reports a 32.8% relative reconstruction error and near-optimal PCA results, with randomized PCA much faster than prcomp().Both randomized and deterministic PCA achieve about 99.18% classification accuracy, while reconstructed digits show virtually no noticeable difference.

5. Randomized robust principal component analysis

Randomized robust PCA separates corrupted data into low-rank and sparse components while reducing the computational burden of iterative optimization. The rsvd implementation performs competitively in simulations, but its advantage depends on matrix shape and target rank.

  • Robust PCA formulation: Robust PCA models a data matrix as low-rank structure plus sparse corruption, optionally including a noise component.The sparse component represents corrupted entries or outliers.
  • Robust PCA formulation: Principal component pursuit can exactly separate low-rank and sparse components under broad assumptions, but robust PCA remains computationally challenging because optimization is iterative.The inexact augmented Lagrange multiplier method is presented as a computationally favorable alternative.
  • Randomized algorithm: The randomized RPCA algorithm replaces the costly SVD inside IALM with randomized SVD, then updates the low-rank and sparse components by rank prediction and soft thresholding.The algorithm alternates updates of the low-rank matrix, sparse matrix, and Lagrange multiplier.
  • Randomized algorithm: The randomized RPCA implementation automatically switches to deterministic SVD when the predicted target rank exceeds min{m, n}/4.This rule avoids using randomized SVD for relatively large target ranks.
  • Handwritten-digit example: On corrupted handwritten digits, the method separates low-rank structure from sparse salt-and-pepper corruption, retaining comparable deterministic and randomized reconstructions.The corruption retains about 10% nonzero elements, and the figures display the recovered low-rank and sparse components.
  • Performance evaluation: The randomized routine has no computational advantage for the tall, thin digit matrix because accurate approximation requires a target rank of 232.Randomized methods become more advantageous for higher-dimensional problems with low-rank structure.
  • Performance evaluation: Across 50 simulation runs, both rsvd randomized and deterministic RPCA outperform the rpca package when the maximum iteration count is 50.The comparison also reports rpca-package performance after 100 iterations.

6. Additional functionality

The section presents CUR and interpolative decompositions in rsvd, emphasizing interpretable row and column subsets and randomized acceleration for low-rank approximation.

  • CUR decomposition: CUR approximates A using C, U, and R, where C and R contain selected columns and rows from the input matrix.The middle factor U is chosen to make the Frobenius-norm approximation error small.
  • CUR decomposition: CUR factors preserve the original data structure, making their rows and columns more interpretable than dense SVD factors.This can also preserve sparsity, which supports compression when sparse SVD factors would reduce storage benefits.
  • CUR decomposition: The rsvd package implements deterministic and randomized CUR decompositions, using interpolative decomposition to construct C and R.The rcur() routine supports both modes, with randomized computation selected by default.
  • Interpolative decomposition: Interpolative decomposition represents A approximately as C Z, with C formed from a small set of columns and Z serving as the interpolation matrix.The selected columns are determined from dominant pivots in a pivoted QR factorization.
  • Interpolative decomposition: Randomized interpolative decomposition first computes a randomized QB decomposition, then applies ID to the smaller matrix B.The paper reports this implementation as more accurate but slightly more computationally demanding than the referenced alternative.
  • Interpolative decomposition: The ID routines support column or row decompositions, with row ID obtained by applying the procedure to the transposed input matrix.Oversampling and power-scheme parameters p and q control the randomized computation.

7. Conclusion

The conclusion positions randomized low-rank methods as scalable tools for big-data computations, trading approximation accuracy for lower computational cost. It also identifies matrix shape, in-memory size, and future sparse-matrix support as important boundaries.

  • Conclusion: Randomized linear-algebra methods offer scalable, tractable computation for low-rank problems that may otherwise be infeasible.Their trade-off is between approximation accuracy and computational cost.
  • Conclusion: Randomized SVD provides strong theoretical error bounds, with approximation quality controlled through oversampling and power iterations.Its computational advantage over other R SVD routines is substantial when the target rank k is relatively small.
  • Conclusion: The performance advantage of randomized methods depends on matrix shape and is less impressive for fat matrices.The conclusion therefore cautions that matrix geometry affects the computational benefit.
  • Conclusion: The rsvd routines target mid-sized problems whose input matrices fit into fast memory.The authors recommend Microsoft R Open to use available computational resources more fully.
  • Conclusion: Planned package developments include randomized routines for additional statistical methods and better support for large-scale sparse matrices.Listed directions include linear discriminant analysis, principal component regression, canonical correlation analysis, and matrix completion.
Loading 1608.02148v5…