fast.svd: Fast Singular Value Decomposition

Description

Any rectangular real matrix M can be decomposed into

$$M = U D V^{'},$$

where $U$ and $V$ are orthogonal matrices with $U' U = I$ and $V' V = I$, and $D$ is a diagonal matrix containing the singular values (see svd). If the matrix $M$ is either "fat" (small n, large p) or "thin" (large n, small p) then the 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$. Using this simple trick 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. Note that unlike the native svd function fast.svd returns only the *positive* singular values (i.e. the dimension of $D$ equals the rank of $M$) and their associated singular vectors. Also note that the latter may differ in sign from those computed by 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)