SVD, pseudoinverse, and PCA: SVD, pseudoinverse and PCA in R
SVD and digital image compression
R has the function svd on board to numerically determine the singular value decomposition. First we calculate the example from the lesson text:
> A <- matrix(c(2, 0, 0, -3), nrow=2, byrow=TRUE); A
[,1] [,2]svdB <- svd(B)
[1,] 2 0
[2,] 0 -3
> svdA <- svd(A)
> U <- svdA$u; U # matrix U in singular value decomposition (SVD)
[,1] [,2]
[1,] 0 1
[2,] 1 0
> D <- svdA$d; D
[1] 3 2
> S <- diag(D) # matrix S in singular value decomposition
> V <- svdA$v; V # matrix V in singular value decomposition
[,1] [,2]
[1,] 0 1
[2,] -1 0
> U %*% S %*% t(V) - A # check of svd(A)
[,1] [,2]
[1,] 0 0
[2,] 0 0The following example of a singular matrix illustrates the difference between exact and numerical arithmetic. \[\begin{aligned}\matrix{1 & -1\\ -1 & 1} &= \matrix{-\frac{1}{2}\sqrt{2} & \frac{1}{2}\sqrt{2} \\ \frac{1}{2}\sqrt{2} & \frac{1}{2}\sqrt{2}}\matrix{2 & 0\\ 0 & 0}\matrix{-\frac{1}{2}\sqrt{2} & \frac{1}{2}\sqrt{2} \\ -\frac{1}{2}\sqrt{2} & -\frac{1}{2}\sqrt{2}}\\ \\ &=\frac{1}{2}\sqrt{2} \matrix{-1 & 1\\ 1 & 1}\matrix{2 & 0\\ 0 & 0}\matrix{-1 & 1\\ -1 & -1}\frac{1}{2}\sqrt{2} \end{aligned} \]
> B <- matrix(c(1, -1, -1, 1), nrow=2, byrow=TRUE); B
[,1] [,2]
[1,] 1 -1
[2,] -1 1
> svdB <- svd(B) # singular value decomposition
> U <- svdB$u; U # matrix U in SVD
[,1] [,2]
[1,] -0.7071068 0.7071068
[2,] 0.7071068 0.7071068
> S <- diag(svdB$d); S # matrix S in SVD
[,1] [,2]
[1,] 2 0
[2,] 0 0
> V <- svdB$v; V # matrix V in SVD
[,1] [,2]
[1,] -0.7071068 0.7071068
[2,] 0.7071068 0.7071068
> U %*% S %*% t(V) - B
[,1] [,2]
[1,] -3.330669e-16 4.440892e-16
[2,] 0.000000e+00 -2.220446e-16
> round(.Last.value, 3) # rounding of numbers
[,1] [,2]
[1,] 0 0
[2,] 0 0Now an example of image compression using a pattern consisting of white and black boxes on a \(9\times 9\) grid. We specify the figure with a matrix consisting of zeros (for white boxes) and ones (for black boxes): \[M = \matrix{ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 1 & 0 & 1 & 1 &1 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1\\ 0 & 1 & 0 & 1 & 1 &1 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0}\] In R we can visualise this pattern.
> M <- matrix(c(0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 1, 0, 1, 0, 0, 0,
0, 0, 1, 0, 0, 0, 1, 0, 0,
0, 1, 0, 1, 1, 1, 0, 1, 0,
1, 0, 0, 1, 0, 1, 0, 0, 1,
0, 1, 0, 1, 1, 1, 0, 1, 0,
0, 0, 1, 0, 0, 0, 1, 0, 0,
0, 0, 0, 1, 0, 1, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0), nrow=9, byrow=TRUE) # patroon
> dev.new(width=240, height=240) % square graphical output
> # pattern as gray picture (dark gray=255, light gray=0)
> image(255*M, col=gray.colors(2, rev=TRUE), axes=FALSE)
Use col=gray.colors(2, rev=TRUE, start=0, end=1) when you prefer a black-and-white picture
Suppose we want to send someone this pattern. Then we can send the matrix \(M\), that is. 81 numbers one and zero. But you can send fewer numbers when you use the singular value decomposition, while the image can still be reconstructed.
The eigenvalues of the matrix \(M\) are as follows in increasing order:
> ev <- eigen(M)
> round(sort(ev$values), 3)
[1] -2.000 -0.903 0.000 0.000 0.000 0.000 1.194
[8] 2.000 3.709
Thus, the singular values of \(M\) are in descending order, rounded to 1 decimal place, equal to \(3.7\), \(2.0\), \(1.2\), \(0.9\) and \(0\).
> svdM <- svd(M)
> round(svdM$d, 4)
[1] 3.7093 2.0000 2.0000 1.1939 0.9032 0.0000 0.0000
[8] 0.0000 0.0000
Suppose we only use the highest singular value \(3.7\). Then the picture looks already a bit like the original figure:
> S1 <- diag(c(S[1,1],0,0,0,0,0,0,0,0))
> M1 <- round(U %*% S1 %*% t(V))
> dev.new(width=240, height=240)
NULL
> image(255*M1, col=gray.colors(2, rev=TRUE), axes=FALSE)To do so we only need to transmit, in short SVD notation, \(1+2\times 9 = 19\) numeric values, namely the singular value and the two vectors \(\vec{u}_1\) and \(\vec{v}_1\).
If we also include the singular value \(2\), we get even back the original figure.
> M3 <- round(U[,1:3] %*% S[1:3,1:3] %*% t(V[,1:3]))
> dev.new(width=240, height=240)
> image(255*M3, col=gray.colors(2, rev=TRUE), axes=FALSE)
Now it suffices to transmit \(3+6\times 9 = 57\) numerical values, namely \(3\) eigenvalues and \(2\times 3= 6\) vectors. This does not seem so much less than \(81\), it is already a reduction of approximately \(30\%\)
\(\phantom{xx}\)
The gain is more spectacular when transmitting large digital images, for example in the picture commonly used in computer vision from \(512\times 512\) pixels.
We read the image into the software environment, calculate the singular value decomposition of the corresponding \(512\times 512\) matrix with grayscale values, and approximate the original image with 8, 16, 32, 64, and 128 terms from the singular value decomposition.
> library(imager) # load R package for image processing
> Lena<- load.image("G:/lena512.bmp") # load bitmap picture
> dev.new(height=512, width=512)
NULL
> plot(Lena, axes=FALSE) # display the picture
> svdLena <- svd(scale(Lena)) # SVD of scaled matrix
> # definition of function to use only highest singular values
> compressLena <- function(dim) {
U <- svdLena $u [,1:dim];
S <-diag(svdLena $d )[1:dim,1:dim];
V <- svdLena $v [,1:dim];
plot(as.cimg(U %*% S %*% t(V)), axes=FALSE)
}
> # display various compressed images
> compressLena(8)
> compressLena(16)
> compressLena(32)
> compressLena(64)
> compressLena(128)Below we present the images for different values of the number of terms \(k\) extracted from the SVD. Clearly less information you need while still having a picture that differs on first sight little form the original (at \(k=64\) the compression factor is about \(\frac{1}{4}\)). With colour pictures, the gain is even greater.
| \(k=8\) | \(k=16\) |
 |
 |
| \(k=32\) | \(k=64\) |
 |
 |
| \(k=128\) | original: \(k=512\) |
 |
 |
