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    0

The 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    0

Now 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

We approximate \(M\) now with the first terms of the singular value decomposition. But instead of using the numerical values of the decomposition for image display we round the elements of the matrix elements to whole numbers \(0\) and \(1\) to reconstruct the pattern.

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