fast.svd: Fast Singular Value Decomposition

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

matrix

tol

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

Value

A list with the following components:
da vector containing the positive singular values
ua matrix with the corresponding left singular vectors
va matrix with the corresponding left singular vectors

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}$$

Examples

Run this code

# load corpcor library
library("corpcor")


# 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)