Randomized Singular Value Decomposition (rsvd).
The randomized SVD computes the near-optimal low-rank approximation of a rectangular matrix using a fast probablistic algorithm.
rsvd(A, k = NULL, nu = NULL, nv = NULL, p = 10, q = 2, sdist = "normal")
A |
array_like; |
k |
integer; |
nu |
integer, optional; |
nv |
integer, optional; |
p |
integer, optional; |
q |
integer, optional; |
sdist |
string c( 'unif', 'normal', 'rademacher'), optional; |
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 diag(d) t(V)
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.
rsvd returns a list containing the following three components:
array_like;
singular values; vector of length (k).
array_like;
left singular vectors; (m, k) or (m, nu) dimensional array.
array_like;
right singular vectors; (n, k) or (n, nv) dimensional array.
The singular vectors are not unique and only defined up to sign (a constant of modulus one in the complex case). If a left singular vector has its sign changed, changing the sign of the corresponding right vector gives an equivalent decomposition.
N. Benjamin Erichson, erichson@berkeley.edu
[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. doi: 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).
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 errorPlease choose more modern alternatives, such as Google Chrome or Mozilla Firefox.