fast.svd {GeneTS}R Documentation

Efficient Computation of the Singular Value Decomposition for Fat and Thin Matrices

Description

Any rectangular real matrix M can be decomposed as

M = U D V',

where U and V are orthogonal, V' means V transposed, and D is a diagonal matrix with the singular values (see svd).

If the matrix M is either "fat" (small n, large p) or "thin" (large n, small p) then then decomposition of M into the U, D, V matrices can be greatly sped up by computing the SVD of either M M' (fat matrices) or M' M (thin matrices), rather than that of M.

This is the trick that allows fast.svd to be substantially faster than the native svd for fat and thin matrices. For near-square matrices fast.svd simply uses the standard svd function.

In contrast to the native svd function note that fast.svd returns only *positive* singular values (i.e. rank D equals rank M) and their associated singular vectors. Also note that the latter may differ in sign from svd.

Usage

fast.svd(m, tol)

Arguments

m matrix
tol tolerance - singular values larger than tol are considered non-zero (default value: tol = max(dim(m))*max(D)*.Machine$double.eps)

Details

For "fat" M (small n, large p) the SVD decomposition of M M' yields

M M' = U D^2 U'

As the matrix M M' has dimension n x n only, this is faster to compute than SVD of M. The V matrix is subsequently obtained by

V = M' U D^(-1)

Similarly, for "thin" M (large n, small p), the decomposition of M' M yields

M' M = V D^2 V'

which is also quick to compute as M' M has only dimension p x p. The U matrix is then computed via

U = M V D^(-1)

Value

A list with the following components:

d a vector containing the *positive* singular values
u a matrix with the corresponding left singular vectors
v a matrix with the corresponding left singular vectors

Author(s)

Korbinian Strimmer (http://www.stat.uni-muenchen.de/~strimmer/).

See Also

svd, solve.

Examples

# load GeneTS library
library(GeneTS)

# generate a "fat" data matrix
n <- 50
p <- 5000
X <- matrix(rnorm(n*p), n, p)

# compute SVD
system.time( s1 <- svd(X) ) 
system.time( s2 <- fast.svd(X) )

eps <- 1e-10
sum(abs(s1$d-s2$d) > eps)
sum(abs(abs(s1$u)-abs(s2$u)) > eps)
sum(abs(abs(s1$v)-abs(s2$v)) > eps)

[Package GeneTS version 2.3 Index]