fast.svd: Fast Singular Value Decomposition

Description

fast.svd returns the singular value decomposition of a rectangular real matrix $$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). The main difference to the native version svd is that fast.svd is substantially faster for "fat" (small n, large p) and "thin" (large n, small p) matrices. In this case the decomposition of $M$ can be greatly sped up by first computing the SVD of either $M M'$ (fat matrices) or $M' M$ (thin matrices), rather than that of $M$. A second difference to svd is that fast.svd only returns the positive singular values (thus the dimension of $D$ always equals the rank of $M$). Note that the singular vectors computed by fast.svd 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)