rsvd: Randomized singular value decomposition

In rsvd: Randomized Singular Value Decomposition

Randomized singular value decomposition

Description

Approximate the leading singular values and singular vectors of a matrix. Random matrix products identify a subspace in which the singular value decomposition (SVD) is computed. This can reduce the factorization work when the requested rank is small relative to the matrix dimensions.

Usage

rsvd(A, k = NULL, nu = NULL, nv = NULL, p = 10, q = 2, sdist = "normal")

Arguments

A	A real or complex numeric matrix or a numeric data frame. Sparse and dense matrices from the Matrix package are supported; see Details. The input must have positive dimensions and finite entries.
k	A positive integer giving the target rank. The default, NULL, uses the smaller matrix dimension; larger requests are capped at that value. Use a rank much smaller than both dimensions for the main computational savings.
nu	An integer giving the number of left singular vectors to return. The default, NULL, returns k vectors. Integer requests are clamped to the range from zero to the effective target rank; zero omits u.
nv	An integer giving the number of right singular vectors to return. The default, NULL, returns k vectors. Integer requests are clamped to the range from zero to the effective target rank; zero omits v.
p	A non-negative integer oversampling parameter. The default is 10. The sketch uses up to k + p columns, bounded by the matrix dimensions.
q	A non-negative integer giving the number of power (subspace) iterations. The default is 2; zero skips the power iterations.
sdist	A character string selecting the sampling distribution. "normal" (the default) uses standard normal entries, "unif" uses uniform entries on [0,1], and "rademacher" uses equally likely values of -1 and 1. For complex input, the real and imaginary parts are sampled independently.

Details

A low-rank approximation represents the columns of a matrix using a small set of shared directions. When these directions capture the structure needed for an analysis, storing and working with the factors can save both memory and computation. For an m \times n input matrix A, the rank-k reconstruction has the form

A \approx U_k D_k V_k^{H}.

The columns of U_k and V_k are orthonormal and the diagonal of D_k contains the returned singular values. Each pair of left and right vectors contributes one matrix component, weighted by its singular value. The superscript H denotes the conjugate transpose, which is the ordinary transpose for real data. Reconstruction requires both sets of vectors and uses Conj(t(v)) for complex input.

To find these directions, the algorithm samples a subspace of width \ell = \min(k+p,m,n). It refines the sample with q power iterations and intermediate orthogonalization, then computes an SVD of the matrix projected onto that subspace. The large input is accessed through matrix products. When the subspace is small, the final SVD acts on a smaller dense matrix.

Oversampling gives the algorithm extra directions from which to identify the dominant rank-k subspace; only k singular values are returned. Increasing p or q can improve accuracy, at the cost of more computation and, for p, more intermediate storage. The function uses the supplied parameters throughout. It does not adapt them to an error tolerance or estimate a numerical rank.

The accuracy achievable at a given rank depends on the singular values of the input. A full-rank matrix can have a useful low-rank approximation when these values decay rapidly. With a flat spectrum and a small target rank, even an optimal approximation leaves a large residual. If \sigma_j are the singular values of the original matrix in decreasing order, its optimal rank-k Frobenius error is

\left(\sum_{j>k}\sigma_j^2\right)^{1/2}.

This expression measures the contribution of the discarded components. It is zero when k reaches the full available rank. Use this truncation error as a reference for the randomized reconstruction residual, with the same preprocessing and target rank for both approximations. An arbitrary full-rank matrix cannot in general be recovered exactly at a smaller rank.

The decomposition acts directly on the supplied matrix, without centering or scaling. See rpca for principal component analysis with those preprocessing options.

Real-valued sparse Matrix inputs are multiplied by thin dense factors without converting the entire input to a dense matrix. Logical and pattern matrices are treated as zero/one data. Complex input is supported for ordinary R matrices. Dense storage is still needed for the random samples, intermediate bases, projected matrices and returned singular vectors. Explicitly reconstructing the approximation also creates a dense matrix.

Omitting k requests all available singular values, even if nu or nv requests only a few vectors. The resulting subspace no longer provides the storage savings of a small-rank approximation. When the target rank is a substantial fraction of the matrix dimensions, a deterministic method may be faster. Use svd as a reference for a full SVD.

Value

A list of class "rsvd" with the following components. Here k denotes the target rank after clamping to the matrix dimensions.

d

A numeric vector of length k containing the approximate singular values in decreasing order. The values are non-negative.

u

An ordinary R matrix with nrow(A) rows and nu orthonormal columns, after applying the default and clamping rules above. Omitted when no left singular vectors are requested.

v

An ordinary R matrix with ncol(A) rows and nv orthonormal columns, after applying the default and clamping rules above. Omitted when no right singular vectors are requested.

Note

Singular vectors are defined only up to sign, or a constant of modulus one for complex data; multiplying a corresponding left and right vector by the same factor leaves the reconstruction unchanged. Repeated singular values also allow rotations within their singular subspace, and closely spaced values can make individual directions sensitive. Compare reconstruction residuals or subspaces, since different valid singular vectors can represent the same approximation.

Use set.seed for reproducibility within a fixed numerical environment; identical vectors are not guaranteed across BLAS or package versions.

Author(s)

N. Benjamin Erichson

References

• N. B. Erichson, S. Voronin, S. L. Brunton and J. N. Kutz. 2019. Randomized Matrix Decompositions Using R. Journal of Statistical Software, 89(11), 1-48. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.18637/jss.v089.i11")}.

• N. Halko, P. Martinsson and J. Tropp. 2009. Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions. https://arxiv.org/abs/0909.4061.

See Also

svd, rqb, rpca

Examples

# A low-rank signal with small dense noise.
set.seed(42)
A <- matrix(rnorm(60 * 5), 60, 5) %*% matrix(rnorm(5 * 40), 5, 40)
A <- A + 0.01 * matrix(rnorm(60 * 40), 60, 40)
k <- 5L
fit <- rsvd(A, k = k, p = 5, q = 2)
A_approx <- sweep(fit$u, 2, fit$d, "*") %*% Conj(t(fit$v))

# Compare residuals at the same target rank on this small matrix.
d <- svd(A, nu = 0, nv = 0)$d
best <- sqrt(sum(d[-seq_len(k)]^2))
c(randomized = norm(A - A_approx, "F"), optimal = best)

# Complex data require a conjugate transpose for reconstruction.
Z <- A + 1i * A[, rev(seq_len(ncol(A)))]
fit_complex <- rsvd(Z, k = k, p = 5, q = 2)
Z_approx <- sweep(fit_complex$u, 2, fit_complex$d, "*") %*%
  Conj(t(fit_complex$v))
sqrt(sum(Mod(Z - Z_approx)^2) / sum(Mod(Z)^2))

# Sparse input stays sparse; the returned factors are dense.
S <- Matrix::Diagonal(x = c(9, 4, 1, rep(0, 17)))
fit_sparse <- rsvd(S, k = 3, p = 2, q = 1)
fit_sparse$d

rsvd documentation built on Nov. 2, 2026, 5:06 p.m.