rsvd (version 1.0.5)

rsvd: Randomized Singular Value Decomposition (rsvd).

Description

The randomized SVD computes the near-optimal low-rank approximation of a rectangular matrix using a fast probablistic algorithm.

Usage

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

Arguments

A

array_like; a real/complex \((m, n)\) input matrix (or data frame) to be decomposed.

k

integer; specifies the target rank of the low-rank decomposition. \(k\) should satisfy \(k << min(m,n)\).

nu

integer, optional; number of left singular vectors to be returned. \(nu\) must be between \(0\) and \(k\).

nv

integer, optional; number of right singular vectors to be returned. \(nv\) must be between \(0\) and \(k\).

p

integer, optional; oversampling parameter (by default \(p=10\)).

q

integer, optional; number of additional power iterations (by default \(q=2\)).

sdist

string \(c( 'unif', 'normal', 'rademacher')\), optional; specifies the sampling distribution of the random test matrix: \('unif'\) : Uniform `[-1,1]`. \('normal\)' (default) : Normal `~N(0,1)`. \('rademacher'\) : Rademacher random variates.

Value

rsvd returns a list containing the following three components:

d

array_like; singular values; vector of length \((k)\).

u

array_like; left singular vectors; \((m, k)\) or \((m, nu)\) dimensional array.

v

array_like; right singular vectors; \((n, k)\) or \((n, nv)\) dimensional array.

Details

The singular value decomposition (SVD) plays an important role in data analysis, and scientific computing. Given a rectangular \((m,n)\) matrix \(A\), and a target rank \(k << min(m,n)\), the SVD factors the input matrix \(A\) as

$$ A = U_{k} diag(d_{k}) V_{k}^\top $$

The \(k\) left singular vectors are the columns of the real or complex unitary matrix \(U\). The \(k\) right singular vectors are the columns of the real or complex unitary matrix \(V\). The \(k\) dominant singular values are the entries of \(d\), and non-negative and real numbers.

\(p\) is an oversampling parameter to improve the approximation. A value of at least 10 is recommended, and \(p=10\) is set by default.

The parameter \(q\) specifies the number of power (subspace) iterations to reduce the approximation error. The power scheme is recommended, if the singular values decay slowly. In practice, 2 or 3 iterations achieve good results, however, computing power iterations increases the computational costs. The power scheme is set to \(q=2\) by default.

If \(k > (min(n,m)/4)\), a deterministic partial or truncated svd algorithm might be faster.

References

• [1] 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. 10.18637/jss.v089.i11.

• [2] N. Halko, P. Martinsson, and J. Tropp. "Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions" (2009). (available at arXiv https://arxiv.org/abs/0909.4061).

See Also

svd, rpca

Examples

# NOT RUN {
library('rsvd')

# Create a n x n Hilbert matrix of order n,
# with entries H[i,j] = 1 / (i + j + 1).
hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") }
H <- hilbert(n=50)

# Low-rank (k=10) matrix approximation using rsvd
k=10
s <- rsvd(H, k=k)
Hre <- s$u %*% diag(s$d) %*% t(s$v) # matrix approximation
print(100 * norm( H - Hre, 'F') / norm( H,'F')) # percentage error
# Compare to truncated base svd
s <- svd(H)
Hre <- s$u[,1:k] %*% diag(s$d[1:k]) %*% t(s$v[,1:k]) # matrix approximation
print(100 * norm( H - Hre, 'F') / norm( H,'F')) # percentage error

# }