rsvd (version 0.6)

reigen: Randomized Spectral Decomposition of a matrix (reigen).

Description

Computes the approximate low-rank eigendecomposition of a symmetric matrix.

Usage

reigen(A, k = NULL, p = 10, q = 1, sdist = "unif")

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)$.

p

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

q

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

sdist

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

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

.

Value

reigen returns a list containing the following two components: returns a list containing the following two components:

Details

The eigenvalue decomposition plays a central role in data analysis and scientific computing. Randomized eigen (reigen) is a fast algorithm to compute the the approximate low-rank eigenvalue decomposition of $A'A$ given the rectangular $(m,n)$ matrix $A$ using a probablistic algorithm. Given a target rank $k << n$, reigen factors the input matrix $A$ as $A'A = V * diag(d) * V'$. The eigenvectors are the columns of the real or complex unitary matrix $V$. The eigenvalues $d$ are non-negative and real numbers.

The parameter $p$ is a oversampling parameter to improve the approximation. A value between 2 and 10 is recommended and $p=10$ is set as 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 archive 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 eigen 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

rsvd, rpca, eigen

Examples

library(rsvd)
set.seed(123)

#Create real random test matrix with dimension (m, n) and rank k
m = 10
n = 5
k = 3
A <- matrix(runif(m*k), m, k)
A <- A %*% t(A)
A <- A[,1:n]

AtA = t(A) %*% A

# Randomized low-rank eigenvalue decomposition k=3
reigen.out <- reigen(A, k=3)
AtA.re = reigen.out$vectors %*% diag(reigen.out$values) %*% t(reigen.out$vectors)
100 * norm( AtA - AtA.re, 'F') / norm( AtA,'F') #Percentage reconstruction error << 1e-8
print(reigen.out$values) # print eigenvalues

# Compare with the deterministic eigenvalue decomposition
eigen.out <- eigen(AtA)
AtA.re2 = eigen.out$vectors %*% diag(eigen.out$values) %*% t(eigen.out$vectors)
100 * norm( AtA - AtA.re2, 'F') / norm( AtA,'F') #Percentage reconstruction error << 1e-8
print(eigen.out$values) # print eigenvalues
all.equal(eigen.out$values[1:k], reigen.out$values)