Version info: Code for this page was tested in R Under development (unstable) (2012-07-05 r59734)
On: 2012-08-08
With: knitr 0.6.3
Singular value decomposition (SVD) is a type of matrix factorization.
For more details on SVD, the Wikipedia page
is a good starting point. On this page, we provide four examples of data
analysis using SVD in R.
Example 1: SVD to find a generalized inverse of a non-full-rank matrix
For a square matrix A with a non-zero determinant, there exists an inverse matrix B such that AB = I and
BA = I. For a matrix that is not square, generalized inverse matrices have
some (but not all) of the properties of an inverse matrix. SVD can be used to
find a generalized inverse matrix. In the example below, we use SVD to
find a generalized inverse B to the matrix A such that ABA = A. We compare our
generalized inverse with the one generated by the ginv command.
library(MASS) a <- matrix(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1), 9, 4) a.svd <- svd(a) a.svd$d
## [1] 3.46e+00 1.73e+00 1.73e+00 1.92e-16
ds <- diag(1/a.svd$d[1:3]) u <- a.svd$u v <- a.svd$v us <- as.matrix(u[, 1:3]) vs <- as.matrix(v[, 1:3]) (a.ginv <- vs %*% ds %*% t(us))
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] ## [1,] 0.0833 0.0833 0.0833 0.0833 0.0833 0.0833 0.0833 0.0833 ## [2,] 0.2500 0.2500 0.2500 -0.0833 -0.0833 -0.0833 -0.0833 -0.0833 ## [3,] -0.0833 -0.0833 -0.0833 0.2500 0.2500 0.2500 -0.0833 -0.0833 ## [4,] -0.0833 -0.0833 -0.0833 -0.0833 -0.0833 -0.0833 0.2500 0.2500 ## [,9] ## [1,] 0.0833 ## [2,] -0.0833 ## [3,] -0.0833 ## [4,] 0.2500
# using the function ginv defined in MASS ginv(a)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] ## [1,] 0.0833 0.0833 0.0833 0.0833 0.0833 0.0833 0.0833 0.0833 ## [2,] 0.2500 0.2500 0.2500 -0.0833 -0.0833 -0.0833 -0.0833 -0.0833 ## [3,] -0.0833 -0.0833 -0.0833 0.2500 0.2500 0.2500 -0.0833 -0.0833 ## [4,] -0.0833 -0.0833 -0.0833 -0.0833 -0.0833 -0.0833 0.2500 0.2500 ## [,9] ## [1,] 0.0833 ## [2,] -0.0833 ## [3,] -0.0833 ## [4,] 0.2500
Example 2: Image processing
The code below requires the ReadImages package. It reads in a jpeg (pansy.jpg)
and plots it in R, first in color (when the image is stored as three
matrices–one red, one green, one blue) and then in grayscale (when the image is
stored as one matrix). Then, using SVD, we can essentially compress the image.
Note that we can recover the image to varying degrees of detail as we recreate
the image from different numbers of dimensions from our SVD matrices. You can
see how many dimensions are needed before you have an image that cannot be
differentiated from the original.
library(ReadImages) x <- read.jpeg("pansy.jpg") dim(x)
## [1] 600 465 3
plot(x, useRaster = TRUE)
r <- imagematrix(x, type = "grey") plot(r, useRaster = TRUE)
r.svd <- svd(r) d <- diag(r.svd$d) dim(d)
## [1] 465 465
u <- r.svd$u v <- r.svd$v plot(1:length(r.svd$d), r.svd$d)
# first approximation u1 <- as.matrix(u[-1, 1]) v1 <- as.matrix(v[-1, 1]) d1 <- as.matrix(d[1, 1]) l1 <- u1 %*% d1 %*% t(v1) l1g <- imagematrix(l1, type = "grey") plot(l1g, useRaster = TRUE)
# more approximation depth <- 5 us <- as.matrix(u[, 1:depth]) vs <- as.matrix(v[, 1:depth]) ds <- as.matrix(d[1:depth, 1:depth]) ls <- us %*% ds %*% t(vs) lsg <- imagematrix(ls, type = "grey")
## Warning: Pixel values were automatically clipped because of range over.
plot(lsg, useRaster = TRUE)
Example 3: Principal components analysis using SVD
This example uses the Stata auto dataset. PCA can be achieved using SVD.
Below, we first use the prcomp command in R and then plot the variances
of the principal components (i.e. the square roots of the eigenvalues). These
values can also be found through spectral decomposition on the correlation
matrix or by SVD on the
variable matrix after standardizing each variable.
library(foreign) auto <- read.dta("http://statistics.ats.ucla.edu/stat/data/auto.dta") pca.m1 <- prcomp(~trunk + weight + length + headroom, data = auto, scale = TRUE) screeplot(pca.m1)
# spectral decomposition: eigen values and eigen vectors xvars <- with(auto, cbind(trunk, weight, length, headroom)) corr <- cor(xvars) a <- eigen(corr) (std <- sqrt(a$values))
## [1] 1.738 0.807 0.526 0.225
(rotation <- a$vectors)
## [,1] [,2] [,3] [,4] ## [1,] -0.507 -0.233 0.825 0.0921 ## [2,] -0.522 0.454 -0.268 0.6708 ## [3,] -0.536 0.390 -0.137 -0.7358 ## [4,] -0.428 -0.767 -0.479 -0.0057
# svd approach df <- nrow(xvars) - 1 zvars <- scale(xvars) z.svd <- svd(zvars) z.svd$d/sqrt(df)
## [1] 1.738 0.807 0.526 0.225
z.svd$v
## [,1] [,2] [,3] [,4] ## [1,] 0.507 -0.233 0.825 -0.0921 ## [2,] 0.522 0.454 -0.268 -0.6708 ## [3,] 0.536 0.390 -0.137 0.7358 ## [4,] 0.428 -0.767 -0.479 0.0057
Example 4: Metric multi-dimensional scaling with SVD
This example uses the Stata cerealnut dataset.
Multi-dimensional scaling can also be achieved using SVD. The plots generated
using cmdscale and the coordinates generated from the SVD steps are
mirrored about the x/x1 = 0 axis, but are otherwise identical.
cnut <- read.dta("http://statistics.ats.ucla.edu/stat/data/cerealnut.dta") # centering the variables mds.data <- as.matrix(sweep(cnut[, -1], 2, colMeans(cnut[, -1]))) dismat <- dist(mds.data) mds.m1 <- cmdscale(dismat, k = 8, eig = TRUE) mds.m1$eig
## [1] 1.58e+05 1.09e+05 1.06e+04 3.83e+02 6.98e+01 1.25e+01 5.76e+00 ## [8] 2.22e+00 4.51e-12 4.51e-12 4.12e-12 3.19e-12 3.15e-12 2.30e-12 ## [15] 2.09e-12 1.25e-12 1.13e-12 8.79e-13 2.97e-13 -1.38e-12 -1.45e-12 ## [22] -1.57e-12 -1.88e-12 -5.92e-12 -2.52e-11
mds.m1 <- cmdscale(dismat, k = 2, eig = TRUE) x <- mds.m1$points[, 1] y <- mds.m1$points[, 2] plot(x, y) text(x + 20, y, label = cnut$brand)
# eigenvalues xx <- svd(mds.data %*% t(mds.data)) xx$d
## [1] 1.58e+05 1.09e+05 1.06e+04 3.83e+02 6.98e+01 1.25e+01 5.76e+00 ## [8] 2.22e+00 1.58e-11 9.90e-12 7.19e-12 4.71e-12 4.15e-12 3.03e-12 ## [15] 2.77e-12 2.08e-12 1.97e-12 1.50e-12 1.26e-12 1.05e-12 7.93e-13 ## [22] 7.35e-13 2.09e-13 1.65e-13 8.38e-14
# coordinates xxd <- xx$v %*% sqrt(diag(xx$d)) x1 <- xxd[, 1] y1 <- xxd[, 2] plot(x1, y1) text(x1 + 20, y1, label = cnut$brand)
