PaperlayerView exact revision on arXiv

Source-linked AI summary

Randomized methods for matrix computations

Per-Gunnar Martinsson

arXiv:1607.01649v3math.NAcs.NA

TL;DR

Low-rank matrix approximation supports storage, computation, and data analysis, while basic randomized SVD does not improve asymptotic cost for dense matrices without fast matrix-vector products. The text presents randomized algorithms and structured projections that reduce computational cost and communication, and can support single-pass processing.

  • Problem

    For dense matrices given as arrays, basic RSVD costs O(mnk) and has no asymptotic edge over classical methods without fast matrix-vector products.

    Source: §6 · p. 12

  • Method

    The text develops randomized projections for low-rank factorization, including structured random matrices for efficient range finding and single-pass variants.

    Source: §5 · p. 9, §6 · p. 12

  • Results

    Randomized methods can attain O(mn log k + k2(m + n)) rank-k approximation cost, require less communication, and enable single-pass factorization.

    Source: §1.3 · p. 5

  • Takeaways & Limitations

    Low-rank approximations can reduce storage from mn entries to k(m+n) entries and support faster matrix-vector products and data interpretation.

    Source: §1.1 · p. 4

    For CUR, substantial singular-value decay can make C and R ill-conditioned, so evaluating the factorization formula can cause substantial round-off accuracy loss.

    Source: §11.1 · p. 22

Abstract

from arXiv · show

The purpose of this text is to provide an accessible introduction to a set of recently developed algorithms for factorizing matrices. These new algorithms attain high practical speed by reducing the dimensionality of intermediate computations using randomized projections. The algorithms are particularly powerful for computing low-rank approximations to very large matrices, but they can also be used to accelerate algorithms for computing full factorizations of matrices. A key competitive advantage of the algorithms described is that they require less communication than traditional deterministic methods.

1. INTRODUCTION

The text introduces randomized methods for low-rank matrix factorizations, emphasizing practical efficiency, reduced communication, and extensions to full factorizations. It sketches the randomized projection idea and the mathematical tools used throughout.

1.1. Scope and objectives.

  • The text covers randomized low-rank EVD, SVD, interpolative, and CUR factorizations, with later sections also addressing randomized acceleration of full factorizations.The stated focus is target rank k much smaller than the matrix dimensions, with additional full-factorization applications discussed near the end.Source: §1.1 · p. 4

1.2. The key ideas of randomized low-rank approximation.

  • Gaussian random sampling constructs a low-rank basis, after which projection onto that basis yields an approximation close to optimal in many important situations.The text describes computing factors from a random sample and using the pseudoinverse to obtain the second factor.Source: §1.2 · p. 5

1.3. Advantages of randomized methods.

  • Randomized methods can reduce rank-k approximation cost from O(mnk) to O(mn log k + k2(m + n)) [25]*Sec. 6.1, and can reduce communication and enable single-pass processing.The text also notes efficient implementation on GPUs and distributed platforms and [23]*Ch. 4 and.Source: §1.3 · p. 5

1.4. Very brief literature survey.

  • The presentation emphasizes practical usage and communication reduction, while situating the methods among earlier theoretical work on performance guarantees and error bounds.It identifies as early work and describes the presentation as largely following the 2011 survey.Source: §1.4 · p. 5

2.1. Notation.

  • The text fixes conventions for vector and matrix norms, submatrix indexing, transposes, and orthonormal matrices.It uses the Golub and Van Loan notation for submatrices.Source: §2.1 · p. 6

2.2. The singular value decomposition (SVD).

  • The SVD writes A as U D V∗, with orthonormal singular-vector matrices and diagonal singular values ordered from largest to smallest.The truncated SVD gives the best rank-k approximation errors in the spectral and Frobenius norms [13].Source: §2.2 · p. 6

2.3. Orthonormalization.

  • The orth function denotes orthonormalization of a matrix’s columns, typically implemented with an unpivoted QR factorization.The text notes that avoiding pivoting makes efficient implementation easier.Source: §2.3 · p. 6

2.4. The Moore-Penrose pseudoinverse.

  • The Moore-Penrose pseudoinverse generalizes the inverse to rectangular or rank-deficient matrices, and its products with A project onto A’s row and column spaces.For a square nonsingular matrix, A† = A−1.Source: §2.4 · p. 7

3. A TWO-STAGE APPROACH

The factorization is split into randomized range finding and deterministic post-processing. This concentrates the difficult computation in Stage A while making Stage B exact up to floating-point arithmetic, subject to precision limits.

3. A TWO-STAGE APPROACH

  • Stage A constructs an orthonormal basis Q whose columns approximately span the column space of A, so that A ≈ QQ∗A.The target rank k is assumed known at this stage.Source: §3 · p. 7
  • Stage B computes B = Q∗A, takes its SVD, and forms U = QÛ to obtain the requested factors using deterministic techniques.When k is much smaller than both matrix dimensions, the matrices in this post-processing stage have at most k rows or columns.Source: §3 · p. 7
  • Stage B is exact up to floating-point arithmetic, so factorization error is incurred in Stage A unless ε is close to machine precision.If Q satisfies ∥A − QQ∗A∥ ≤ ε, the resulting factorization inherits that error bound.Source: §3 · pp. 7–8
  • The residual bound controls absolute errors in the leading singular values, but does not guarantee relative accuracy or strong approximation of the singular vectors.The stated absolute bound is |σi − D(i, i)| ≤ ε for i = 1, 2, . . . , k.Source: §3 · p. 8

4. A RANDOMIZED ALGORITHM FOR “STAGE A” — THE RANGE FINDING PROBLEM

The range-finding method samples Gaussian directions and uses modest oversampling to form a basis that captures dominant singular directions. The text then extends randomized factorization to single-pass settings, with distinct constructions for Hermitian and general matrices.

4. A RANDOMIZED ALGORITHM FOR “STAGE A” — THE RANGE FINDING PROBLEM

  • With probability almost 1, taking k + 10 samples instead of k can produce a basis comparable to the best possible basis, avoiding pollution from modes beyond the leading k.The approximation error is random; bounds show it is close to optimal when singular values decay reasonably fast.Source: §4 · p. 8

5. SINGLE PASS ALGORITHMS

  • A single-pass modification accesses each matrix entry only once, allowing factorization when A is too large to store, though it tends to be less accurate than the two-pass algorithm.The text says the two-pass method is generally preferable when available, for moderately higher cost.Source: §5 · p. 9

5.1. Hermitian matrices.

  • For a Hermitian matrix, the method samples Y = AG, forms an orthonormal basis Q, and estimates the reduced matrix through an overdetermined system.The system for C is overdetermined because C must be Hermitian; ill-conditioning of Q∗G also reduces accuracy.Source: §5.1 · p. 10

5.2. General matrices.

  • For a general matrix, independent samples of A and A∗ form approximate column- and row-space bases, and a least-squares solution determines the reduced matrix C.The system is overdetermined by about a factor of two, with more aggressive oversampling suggested.Source: §5.2 · p. 11

6. A METHOD WITH COMPLEXITY O(mn log k) FOR GENERAL DENSE MATRICES

For dense matrices stored as arrays, the basic RSVD does not improve asymptotic complexity without fast matrix-vector products. A structured random projection and a modified Stage B reduce the overall cost to O(mn log k), though error theory for the accelerated range finder remains limited.

6. A METHOD WITH COMPLEXITY O(mn log k) FOR GENERAL DENSE MATRICES

  • Replacing Gaussian G with a structured random matrix can make AΩ computable in O(mn log(k)) while retaining enough randomness for its columns to span A’s range accurately.The construction Ω = D F S uses random unit-modulus diagonal phases, a discrete Fourier transform, and random column subsampling; ℓ = 2k is often a good choice.Source: §6 · p. 12
  • The structured projection reduces Stage A from O(mnk) to O(mn log(k)), and Stage B must also avoid directly computing Q∗A to attain that overall cost.The text offers single-pass processing or, typically better, row extraction as Stage B options.Source: §6 · pp. 12–13
  • The accelerated range finder has weaker proven error guarantees than Gaussian sampling, and the text lacks good theory to predict precisely when practice matches Gaussian accuracy.It reports that the accelerated scheme is often as accurate as the Gaussian one in practice.Source: §6 · p. 13

7. THEORETICAL PERFORMANCE BOUNDS

The section bounds the randomized SVD’s error in expectation and in the probability of large deviations. It finds that Frobenius-norm accuracy can be near optimal, while spectral-norm accuracy depends strongly on singular-value decay; a modest-cost modification can improve both.

7. THEORETICAL PERFORMANCE BOUNDS

  • The factors depend on the random draw of G, so the error is a random variable and the analysis bounds its expectation and the likelihood of large deviations.The bounds concern ∥A − QQ* A∥ because the Stage B post-processing is exact up to floating-point arithmetic.Source: §7 · p. 13

7.1. Bounds on the expectation of the error.

  • With p=10, the Frobenius-norm error is within a factor of 1 + k/9 of the theoretically minimal error; with p=k+1, it is within a factor of 2.The theoretically minimal error is a lower bound for any rank-k approximant.Source: §7.1 · p. 13
  • For spectral-norm error, the bound’s second term involves the minimal Frobenius-norm error, which can be much larger, especially for large m or n.The text characterizes this sub-optimality as a true limitation of the basic randomized scheme.Source: §7.1 · pp. 13–14
  • Fast decay of the tail singular values makes spectral- and Frobenius-norm errors similar, while slow decay can leave spectral error poor even when Frobenius error remains fine.With exponential decay, the tail contribution is modest when β is not close to 1; with no decay, the spectral-norm suboptimality can be devastating for large matrices.Source: §7.1 · p. 14
  • A modified RSVD can make errors close to optimal in both spectral and Frobenius norms, at a modest increase in computational cost.The section points to Section 8 and [25]*Sec. 4.5 for this modification.Source: §7.1 · p. 14

7.2. Bounds on the likelihood of large deviations.

  • For p ≥ 4, the likelihood of a large deviation has failure probability at most 3 e^-p.The text notes that this probability decays extraordinarily fast with the over-sampling parameter p.Source: §7.2 · p. 14

8. AN ACCURACY ENHANCED RANDOMIZED SCHEME

Power iteration reshapes the singular-value decay to improve randomized SVD accuracy, while inter-step orthonormalization guards against rounding loss. The section also gives a bound on the accuracy gain and an extended sampling option that uses more computed vectors at higher cost.

8.1. The key idea — power iteration.

  • The basic randomized scheme can be suboptimal for slowly decaying singular values, with its spectral-norm error factor potentially much larger than the minimal σk+1 error.The basic scheme gives accurate results when singular values decay rapidly.Source: §8.1 · p. 14
  • Power iteration preserves A’s left singular vectors while raising its singular values to 2q+1, making even slowly decaying spectra decay fast enough for the approximation.The method samples A(q) rather than A directly.Source: §8.1 · p. 14
  • Orthonormalizing between power iterations prevents rounding errors from erasing information about smaller singular values, but adds orthogonalization calls.Unpivoted Gram–Schmidt can make those extra calls relatively fast.Source: §8.1 · p. 15, Figure 5 · p. 16

8.2. Theoretical results.

  • As q increases, the power scheme drives the factor multiplying σk+1 toward one exponentially fast, reducing the expected degree of sub-optimality.The simplified bound considers a worst case with no singular-value decay beyond the truncation point.Source: §8.2 · p. 16

8.3. Extended sampling matrix.

  • The extended sampling matrix combines AG, A²G, …, A^qG into qℓ columns, using samples the earlier scheme leaves unused.Its QR, matrix-multiplication, and SVD costs grow by O(q²), so it is chiefly useful when applications of A and A* dominate.Source: §8.3 · p. 16

9. THE NYSTR ¨OM METHOD FOR SYMMETRIC POSITIVE DEFINITE MATRICES

For symmetric positive definite matrices, the Nyström method produces a low-rank eigenvalue decomposition. It exploits symmetry to improve accuracy at comparable cost to the basic randomized SVD, using two applications of A to k+p vectors.

9. THE NYSTR ¨OM METHOD FOR SYMMETRIC POSITIVE DEFINITE MATRICES

  • The procedure targets an approximate eigenvalue decomposition of a nonnegative matrix, with an orthonormal eigenvector factor and a nonnegative diagonal factor.The text frames the Nyström method for SPD input.Source: §9 · pp. 16–17
  • The method computes B1 = AQ and B2 = Q*B1, factors the SPD B2 as C*C, and forms F = B1C−1 by triangular solve.F is an approximate Cholesky factor of A.Source: §9 · p. 17
  • For SPD A, Nyström tends to produce substantially more accurate results than the direct formula at comparable cost, exploiting symmetry for one power-iteration step “for free.”Both approaches apply A twice to k+p vectors.Source: §9 · p. 17, Figure 6 · p. 17

10. RANDOMIZED ALGORITHMS FOR COMPUTING THE INTERPOLATORY DECOMPOSITION (ID)

The section presents column-, row-, and double-sided interpolative decompositions, whose selected rows or columns preserve useful structure and interpretation. It relates deterministic ID to pivoted QR and develops randomized sampling methods for computing ID efficiently.

10.1. Structure preserving factorizations.

  • The ID represents A using selected columns and a coefficient matrix; selected columns retain properties such as sparsity or non-negativity and support data interpretation.The section also notes memory savings and preservation of problem physics in some discretization settings.Source: §10.1 · pp. 17–18
  • When A is not exactly rank k, the best rank-k ID can have substantially greater error than the theoretically minimal rank-k approximation.The ID and column-pivoted QR attain the same minimal error as each other.Source: §10.1 · p. 18

10.2. Three flavors of ID: The row, column, and double-sided ID.

  • The row ID represents A using selected rows, while the double-sided ID uses selected rows and columns together with their intersection submatrix.The row factor contains an identity matrix at the selected-row indices.Source: §10.2 · p. 18

10.3. Deterministic techniques for computing the ID.

  • A column ID derived from CPQR has the same error as the truncated QR factorization, and standard CPQR software can compute it with light post-processing.The construction identifies the first k pivot columns as C and forms the interpolation coefficients from the QR factors.Source: §10.3 · pp. 18–19

10.4. Randomized techniques for computing the ID.

  • For randomized ID, q = 0 is fast and accurate with rapid singular-value decay, while larger q can improve accuracy at the cost of longer execution time.The figure suggests q = 1 or 2 for slowly decaying singular values.Source: Figure 8 · p. 21
  • Randomized sampling builds a basis Y for A’s column space, from which a row ID of A can be obtained; structured sampling reduces the stated complexity to O(mn log k).The Gaussian-sampling algorithm has complexity O(mnk).Source: §10.4 · p. 21

11. RANDOMIZED ALGORITHMS FOR COMPUTING THE CUR DECOMPOSITION

The CUR decomposition approximates a matrix using selected columns and rows, preserving useful structure and supporting data interpretation. The section constructs CUR from a double-sided ID and notes that conditioning can compromise accuracy.

11.1. The CUR decomposition.

  • CUR approximates A as C U R, with C and R selected from A’s columns and rows, respectively.Like ID, CUR can preserve sparsity or non-negativity, reduce memory needs, and aid data interpretation.Source: §11.1 · p. 22
  • A common randomized CUR procedure weights rows and columns by leverage scores, samples them, then computes U using pseudoinverses of C and R.The weights may use row and column norms or leading singular vectors.Source: §11.1 · p. 22
  • When the singular values decay substantially, ill-conditioning of C and R can make computing U lose accuracy through accumulated round-off; ID avoids this inversion-related problem.The discussion estimates their condition numbers by roughly σ1(A)/σk(A).Source: §11.1 · p. 22

11.2. Converting a double-sided ID to a CUR decomposition.

  • The section converts a double-sided ID by taking C and R from A’s selected columns and rows, then choosing U = ZR† via a least-squares solution.The formulation targets A ≈ C U R.Source: §11.2 · pp. 22–23

12. ADAPTIVE RANK DETERMINATION WITH UPDATING OF THE MATRIX

The section addresses choosing rank from a tolerance by incrementally building a low-rank approximation. It compares greedy choices of added vectors, describes a blocked randomized variant, and explains residual-norm evaluation.

12.1. Problem formulation.

  • Given tolerance ε, the algorithm seeks rank k and factors Qk and Bk with Qk orthonormal and ∥A − QkBk∥ ≤ ε.The exact minimum rank is the smallest k with σk+1 ≤ ε; the described algorithms typically find a fairly close, not necessarily optimal, rank.Source: §12.1 · p. 23, §12.2 · p. 24

12.2. A greedy updating algorithm.

  • Adding one column to Q and one row to B at each step generalizes Gram-Schmidt and maintains the identity A = QjBj + Aj.The process is designed for dense matrices in RAM and directly updates the matrix while guaranteeing the requested approximation tolerance.Source: §12.1 · p. 23, §12.2 · p. 24

Pick the largest remaining column.

  • Choosing the largest remaining column produces column-pivoted Gram-Schmidt, which can work well for rapidly decaying singular values but can perform poorly, especially with slow decay.The section also notes exceptions even under rapid decay and cites the Kahan counter-example.Source: Pick the largest remaining column · p. 24

Pick the locally optimal vector.

  • The locally optimal choice attains theoretically optimal precision, but is impractical because solving its local minimization problem is computationally hard.The supplied passage states that greediness itself is not the problem.Source: Pick the locally optimal vector · p. 25

A randomized selection strategy.

  • Randomly selecting each update vector makes the greedy method equivalent to the basic randomized SVD, often yielding an efficient, near-optimal-accuracy factorization.Power iterations can improve accuracy at some cost in computational efficiency.Source: A randomized selection strategy · p. 25

12.3. A blocked updating algorithm.

  • Blocking draws b Gaussian vectors per iteration and extends Q and B by b columns and rows, respectively.The algorithm uses a block size b as input.Source: §12.3 · p. 25, Figure 12 · p. 25

12.4. Evaluating the norm of the residual.

  • For Frobenius-norm residuals, orthogonality lets the algorithm compute the remainder norm from the input norm and the factor norm; spectral norms can use power iteration or faster, less accurate randomized estimation.The Frobenius computation uses ∥QjBj∥Fro = ∥Bj∥Fro.Source: §12.4 · p. 26

13. ADAPTIVE RANK DETERMINATION WITHOUT UPDATING THE MATRIX

For matrices that cannot be explicitly updated, randomized residual estimates can guide rank selection without forming the residual. A sequence of small projected sample vectors can signal that the current basis captures most of the range, with a probability guarantee whose usefulness depends on singular-value decay.

13. ADAPTIVE RANK DETERMINATION WITHOUT UPDATING THE MATRIX

  • For sparse, out-of-core, or implicitly defined matrices, randomized norm estimation can bound the residual with probability at least 1 − p without explicitly updating A.This supports the given-tolerance problem when updating A is impractical.Source: §13 · p. 26
  • A run of small projected sample vectors z_j suggests that the current basis captures most of A’s range, enabling adaptive stopping.The sequential range finder is mathematically equivalent to the basic RSVD for the same first p Gaussian samples.Source: §13 · p. 26, Figure 13 · p. 27
  • With α = 1/10, r consecutive small z_j vectors imply a residual bound ∥A − Q_jB_j∥ ≤ ε with probability at least 1 − 10^-r.The randomized range finder also gives an ε bound with probability at least 1 − min{m, n}10^-r.Source: §13 · p. 27, Figure 14 · p. 28
  • The stopping rule works well only when singular values decay reasonably fast.The underlying Gaussian-vector norm estimate is stated for r independent standard Gaussian probes and confidence parameter α.Source: §13 · p. 27

14. `A POSTERIORI ERROR BOUNDS, AND CERTIFICATES OF ACCURACY

A thin independent Gaussian probe can estimate the error of a computed low-rank approximation, providing an accuracy certificate for randomized methods. The certificate can address weak guarantees of fast projections, authenticate single-pass results, and check subsampling methods when prior information is uncertain.

14. A POSTERIORI ERROR BOUNDS, AND CERTIFICATES OF ACCURACY

  • A thin Gaussian test matrix provides a highly reliable bound on the error of a computed low-rank approximation.For α = 0.1 and 10 test columns, the stated bound holds with probability at least 1 − 10^-10.Source: §14 · p. 28
  • For fast Johnson–Lindenstrauss projections, certification restores reliability while retaining O(mn log k) complexity, at only slightly higher practical execution time.If the error is unacceptable, the process can often be repeated with a different projection.Source: §14 · p. 29
  • In a single-pass algorithm, a separately drawn thin Gaussian test can warn that the approximation is unreliable, but the matrix cannot be revisited to redo it.The method compares accumulated samples of AG_cert and A_approxG_cert.Source: §14 · p. 29
  • When sampling columns based on prior knowledge of A, certification can authenticate the result if that prior information is uncertain.Sampling probabilities may use column norms or estimated leverage scores; column norms alone generally lack enough information.Source: §14 · p. 29

15. RANDOMIZED ALGORITHMS FOR COMPUTING A RANK-REVEALING QR DECOMPOSITION

Randomized projections are used to accelerate full and partial column-pivoted QR factorizations, including in regimes where the target rank is not small. The HQRRP method selects groups of pivot columns from a small sampled matrix, while downdating keeps its leading flop count equal to classical Householder QR.

15. RANDOMIZED ALGORITHMS FOR COMPUTING A RANK-REVEALING QR DECOMPOSITION

  • Randomized projections can accelerate full factorizations for any rank, including full CPQR, and are reported faster than traditional algorithms across essentially every computational regime.The section presents a high-level introduction to randomized full-factorization algorithms.Source: §15 · p. 29

15.1. Column pivoted QR decomposition.

  • Classical Householder QR is stable and accurate, but its sequence of rank-1 updates is communication intensive and slow on modern hardware.Blocking groups pivot vectors so most flops become faster matrix-matrix operations.Source: §15.1 · p. 30
  • HQRRP selects groups of b pivot columns by performing CPQR on a small Gaussian sample Y = GA, using oversampling parameter p.The QR of Y fits in fast memory; the algorithm repeats the selection on each remaining lower-right block.Source: §15.1 · p. 31, Figure 16 · p. 31
  • Downdating the sampling matrix from one factorization step to the next obviates drawing a new random matrix at each step and reduces HQRRP’s leading asymptotic flop count to 2mn^2 − (4/3)n^3, identical to classical Householder QR.Source: §15.1 · p. 31

16. A STRONGLY RANK-REVEALING UTV DECOMPOSITION

randUTV is a blocked, incremental factorization designed to provide close-to-optimal low-rank approximations and accurate singular-value estimates. Its flexible right transformations help concentrate mass on the diagonal, while randomized subspace approximations and small-block SVDs underpin each step.

16. A STRONGLY RANK-REVEALING UTV DECOMPOSITION

  • randUTV provides close-to-optimal low-rank approximations and highly accurate singular-value estimates, and numerical-rank k inputs can stop early at cost O(mnk).It is blocked for modern hardware and is not iterative, which simplifies software optimization.Source: §16 · pp. 31–32

16.1. The UTV decomposition.

  • A UTV factorization has unitary U and V with triangular T; SVD is the diagonal-T special case, while CPQR is the permutation-V special case.The flexibility of UTV enables efficient updating procedures.Source: §16.1 · p. 32

16.2. An overview of randUTV.

  • randUTV applies blocked unitary transformations and allows right transformations built mostly from Householder reflectors, concentrating more mass on the diagonal and making off-diagonal entries small.The resulting T is driven toward upper triangular form in successive block steps.Source: §16.2 · p. 32, Figure 17 · p. 32

16.3. A single step block factorization.

  • Each step uses randomized power iteration to approximate dominant right singular vectors, then Householder QR and a small-block SVD to form the transformations.The added application of A gives accuracy similar to q + 1/2 power-iteration steps; the resulting diagonal entries typically approximate the leading singular values.Source: §16.3 · p. 33
  • After processing the first b rows and columns, randUTV repeats the procedure on the remaining block until all blocks are processed.The algorithm builds the factorization incrementally.Source: §16.3 · p. 33, Figure 18 · p. 34

Related papers

Loading 1607.01649v3…