rsvd (version 0.6)

rsvd: Randomized Singular Value Decomposition (rsvd).

Description

Compute the near-optimal low-rank singular value decomposition (SVD) of a rectangular matrix.

Usage

rsvd(A, k = NULL, nu = NULL, nv = NULL, p = 10, q = 1, sdist = "unif", vt = FALSE)

Arguments

A

array_like a real/complex input matrix (or data frame), with dimensions $(m, n)$.

k

int, optional determines the target rank of the low-rank decomposition and should satisfy $k << min(m,n)$.

nu

int, optional the number of left singular vectors to be computed. This must be between $0$ and $k$.

nv

int, optional the number of right singular vectors to be computed. This must be between $0$ and $k$.

p

int, optional oversampling parameter for (default $p=10$).

q

int, optional number of power iterations (default $q=1$).

sdist

str $c('normal', 'unif', 'col')$, optional Specifies the sampling distribution. $'unif'$ : (default) Uniform `[-1,1]`. $'normal$' : Normal `~N(0,1)`. $'col'$ : Random column sampling.

vt

bool ($TRUE$, $FALSE$), optional $TRUE$ : returns the transposed right singular vectors $vt$. $FALSE$ : (default) returns the right singular vectors as $v$, this is the format as svd returns $v$.

.............

.

Value

rsvd returns a list containing the following three components: returns a list containing the following three components:

Details

The singular value decomposition (SVD) plays a central role in data analysis and scientific computing. Randomized SVD (rSVD) is a fast algorithm to compute the approximate low-rank SVD of a rectangular $(m,n)$ matrix $A$ using a probablistic algorithm. Given a target rank $k << n$, rsvd factors the input matrix $A$ as $A = U * diag(d) * V'$. The right singluar vectors are the columns of the real or complex unitary matrix $U$ . The left singular vectors are the columns of the real or complex unitary matrix $V$. The singular values $d$ are non-negative and real numbers.

The parameter $p$ is a oversampling parameter to improve the approximation. A value between 5 and 10 is recommended and $p=10$ is set by default.

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

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

References

• [1] N. Halko, P. Martinsson, and J. Tropp. "Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions" (2009). (available at arXiv http://arxiv.org/abs/0909.4061).
• [2] S. Voronin and P.Martinsson. "RSVDPACK: Subroutines for computing partial singular value decompositions via randomized sampling on single core, multi core, and GPU architectures" (2015). (available at `arXiv http://arxiv.org/abs/1502.05366).

See Also

svd, rpca

Examples

# Create a n by 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