Singular Value Decomposition of a Matrix

Dr. Hari V S
Department of Electronics and Communication
College of Engineering Chengannur

What You will Learn

• You will learn to the singular value decomposition and implement it using Python.
• You will understand SVD as a tool for trend analysis and dimensionality reduction.
• You will use SVD as a primitive image compression scheme.

Singular Value Decomposition of a Matrix

A matrix $\mathbf{A}$ can be decomposed as \begin{equation*} \tilde{\mathbf{A}}=\sum_{i=0}^{r}\lambda_{i}U_{i}V_{i}^{T} \end{equation*} where $U_{i}$ and $V_{i}$ are the singular vectors and $\lambda_{i}$ are the eigen values with $\lambda_{i}< \lambda_{j}$ for $i>j$. The parameter $r$ is the number of eigen values of the matrix. Now let us consider a matrix of random integers and its SV decomposition using the code below.
import numpy as np
from pylab import *
A = np.mat(np.random.randint(1000, size=(1000, 1000)))
[a,b,c]=svd(A)
bnorm=b*((max(b))**-1)
figure(1) plot(bnorm[0:10],'*')
plot(bnorm[0:10])
show()
A $1000\times 1000$ matrix of random integers in the range $[0,1000]$ is created using the mat function from the numerical Python (numpy) module and its SVD is computed. The array b contains the $1000$ eigen values in the decreasing order. The arrays a and b contain the $\mathbf{U}$ and $\mathbf{V}$ with each vector corresponding to the eigen value in b. It can be seen that the eigen values are in the decreasing order. Figure below shows the first 10 normalized eigen values. It is observed that the first component is very strong and the remaining values are relatively weak. It means that the data matrix can be represented with the first few eigen values, resulting in good data compression. But it should not be forgotten that the matrix we used is purely random and the correlation between neighbouring elements is almost zero.

Let us look at the compression of the above matrix if only a few (say 10) principal eigen values and compute the Frobenius norm between the matrix and its approximation. The Frobenius norm, expressed as \begin{equation}\zeta=\sum_{i=1}^{N}\sum_{j=1}^{N}|a_{i,j}-\tilde{a_{i,j}}|^{2}\end{equation} is computed and plotted for different rank values using the code below. You can give diffrent values for rho and observe the compression.
import numpy as np
from pylab import *
A = np.mat(np.random.randint(50, size=(1000, 1000)))
a,b,c= np.linalg.svd(A)
bnorm=b*((max(b))**-1)
figure(1)
plot(bnorm[0:10],'*')
plot(bnorm[0:10])
title('Normalized eigen values')
grid('on')
show()
rho=10
frobnormarr=[]
for k in range(0,rho):

approxmat=np.dot(a[:,:k],np.dot(np.diag(b[:k]),c[:k,:]))

frobnorm=linalg.norm(A-approxmat,2)

frobnormarr.append(frobnorm)
figure(2)
plot(frobnormarr)
plot(frobnormarr,'*')
show()

Compression Ratio

You require $10^{6}$ locations for storing a $1000\times 1000$ matrix. Let us look at the compression based on SVD with say 10 principal eigen values. You require 10 locations for the eigen values, $10\times 1000$ for storing $\mathbf{U}$ (Each $\mathbf{U_{i}}$/$\mathbf{V_{i}}$ has $1000$ values) and $10\times 1000$ for storing $\mathbf{V}$, making a total of $20010$ values to store. So
\begin{equation} \text{Compression ratio}=\frac{\text{Size of the uncompressed data}}{\text{Size of the compressed data}}\\ =\frac{10^{6}}{20010}=49.98 \end{equation}

SVD of a Gray Scale Image

An image is a two dimensional array or matrix of numbers. If this image comes out of a camera or scanner, it is called a natural image. If it is a single array of integers in the range $[0,255]$, it is a gray scale image. Color images contain three such arrays each corresponding to the red, green and blue spectral components. Some hyper-spectral images such as satellite images contain 16 such channels.
The code below shows the SVD of a ($512\times512$) gray scale image. The image is read using the scipy.misc module and converted into a matrix. The rest of the steps in compression is similar to the previous code.

import scipy.misc as mis
import numpy as np
from pylab import *
img = mis.ascent()
A=mat(img)
a,b,c= np.linalg.svd(A)
bnorm=b*((max(b))**-1)
rho=40
frobnormarr=[]
for k in range(0,rho):

approxmat=np.dot(a[:,:k],np.dot(np.diag(b[:k]),c[:k,:]))

frobnorm=linalg.norm(A-approxmat,2)

frobnormarr.append(frobnorm)
figure()
subplot(221)
imshow(img)
title("Original image")
gray()
subplot(222)
imshow(approxmat)
title("Compressed image with 40 eigen values")
gray()
subplot(223)
plot(bnorm[0:10],'*')
plot(bnorm[0:10])
title('Normalized eigen values')
grid('on')
subplot(224)
plot(frobnormarr)
plot(frobnormarr,'*')
grid('on')
title("Frobenius norm between the original and compressed images")
show()

The top left figure shows the original image and the top right is the compressed image with the first $40$ eigen values. You may compare and recognize the same picture with loss of some features. The bottom left figure shows the variation of eigen values and the bottom right figure shows the variation in Frobenius norm with rank.

What You Learned

• You learned to implement the singular value decomposition of a matrix using Python.
• You used SVD as tool for trend analysis and dimensionality reduction.
• You learned to use SVD as a primitive image compression scheme.