Image Reconstruction using Singular Value Decomposition (SVD) in Python

In this post, we will explore the use of SVD on Image analysis. We will mainly use SVD on images to get main components/singular vectors capturing the image and use part of them to reconstruct the image.

Singular Value Decomposition (SVD) is one of the commonly used dimensionality reduction techniques. SVD/PCA is the mainstay of common unsupervised learning methodologies in Machine Learning/Data Science.

One of the interesting applications of SVD you may not have heard is image compression and reconstruction. If you think about it a bit, it is not surprising that SVD can capture all the redundancies in an image nicely. Obviously the idea is not new. It has been around for a while and commonly taught at Computer Science courses as an application of SVD before Data Science was a thing.

In this post, we will see step-by-step example of performing SVD on an image and use top singular vectors or principal components to reconstruct it. If you are new to SVD in Python, check out the post on Singular Value Decomposition in Python

Let us load the packages needed to perform SVD on images. In addition to the standard Pandas and NumPy, we need PIL for image manipulation.

import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
import seaborn as sns
%matplotlib inline

We use the NY City image from Wikipedia. For the sake of simplicity, we have downloaded the image and work with it locally.

#ny_file= "https://en.wikipedia.org/wiki/Boroughs_of_New_York_City#/media/File:5_Boroughs_Labels_New_York_City_Map.svg"
ny_file = "5_Boroughs_Labels_New_York_City_Map.svg.png"

We can load the image using Image module’s open() function.

img = Image.open(ny_file)

Let us check the image by opening the image object using show() function from Image module. This will open the image in a separate window.

# view the image from Python
img.show()

We can also check the image size with size(). It will give you the width and height of the image.

>img.size
(1280, 1266)

An image is stored as a matrix. An image in RGB color model stores an image in three matrices one each for Red, Green, and Blue color. In this post on using SVD to reconstruct an image, we will only deal with one of the colo matrices, not all three. However, the same principle holds for dealing with other two matrices as well.

Let us first extract image corresponding to Red color matrix. With PIL we can get a color matrix using getdata() function. We specify band=0 to get red color image.

red_band =img.getdata(band=0)

Let us convert the red color image as a numpy array containing the numerical values corresponding to each pixel i.e. matrix element.

# convert to numpy array 
>img_mat = np.array(list(red_band), float) 
>img_mat.size
1620480

Now we will convert the one dimensional numpy array and 2d-numpy matrix using the image size.

# get image shape
img_mat.shape = (img.size[1], img.size[0])
# conver to 1d-array to matrix
img_mat = np.matrix(img_mat)
img_mat

matrix([[161., 165., 165., ..., 158., 158., 158.],
        [177., 247., 247., ..., 158., 158., 158.],
        [177., 247., 247., ..., 158., 158., 158.],
        ...,
        [158., 158., 158., ..., 158., 158., 158.],
        [158., 158., 158., ..., 158., 158., 158.],
        [158., 158., 158., ..., 158., 158., 157.]])

Let us check how the numpy matrix looks like as image using PIL’s imshow() function.

plt.imshow(img_mat)

Let us also compare how the original RGB image looks in comparison to the image using the single layer of the image by plotting side-by-side.

fig, axs = plt.subplots(1, 2,figsize=(10,10))
axs[0].imshow(img)
axs[0].set_title('Original Image', size=16)
axs[1].imshow(img_mat)
axs[1].set_title(' "R" band image', size=16)
plt.tight_layout()
plt.savefig('Original_image_and_R_band_image_for_SVD.jpg',dpi=150)
Original Image and R Band Image for SVD
Original Image and R Band Image for SVD


Let us center and scale the data before applying SVD. This will help us put each variable in the same scale.
[sourcecode language="Python"]
# scale the image matrix befor SVD
img_mat_scaled= (img_mat-img_mat.mean())/img_mat.std()

We can use NumPy’s linalg module’s svd function to perform singular value decomposition (SVD) on the scaled image matrix as before.

# Perform SVD using np.linalg.svd
U, s, V = np.linalg.svd(img_mat_scaled) 

Performing singular value decomposition (SVD) on matrix will factorize or decompose the matrix in three matrices, U, s, and V. The columns of both U and V matrices are orthonormal and called right and left singular vectors. And the matrix s is a diagonal matrix with only positive numbers and it correspond to eigen values.

Singular Value Decomposition (SVD) in Python

Let us check the dimension of U and V matrices. We can see that both U and V are square matrices and their dimensions matches the image size.

U.shape
(1266, 1266)
V.shape
(1280, 1280)

We can also see that the eigen values are simply a vector here ordered in descending order.

s
array([7.28401576e+02, 5.61698279e+02, 4.94065979e+02, ...,
       5.45091892e-14, 5.45091892e-14, 5.45091892e-14])

We can use eigen values from SVD to compute the amount of variances explained by each singular vector.

# Compute Variance explained by each singular vector
var_explained = np.round(s**2/np.sum(s**2), decimals=3)

The first singular vector or principal component explains most of variation in the image. In this example it explains 32% of the total variation and the second one explains close to 20% of the variation.

# Variance explained top Singular vectors
var_explained[0:20]
array([0.327, 0.195, 0.151, 0.068, 0.041, 0.028, 0.021, 0.015, 0.012,
       0.009, 0.008, 0.008, 0.007, 0.006, 0.006, 0.005, 0.004, 0.004,
       0.004, 0.003])

A bar plot made using variance explained nicely captures how each vector is contributing to the variation in the image.

Here we use Seaborn in Python to make a barplot using variance explained by top 20 singular vectors.

sns.barplot(x=list(range(1,21)),
            y=var_explained[0:20], color="dodgerblue")
plt.xlabel('Singular Vector', fontsize=16)
plt.ylabel('Variance Explained', fontsize=16)
plt.tight_layout()
plt.savefig('svd_scree_plot.png',dpi=150, figsize=(8,6))
#plt.savefig("Line_Plot_with_Pandas_Python.jpg")

Such a plot is called as Scree plot and widely used to guess the number of vectors one needs to capture most of the variations.

Scree Plot: Proportion of variance Explained

Reconstructing Image with top-K Singular vectors

The top K singular vectors captures most of the variation. Therefore instead of using all the singular vectors and multiplying them as shown in SVD decomposition, we can reconstruct the image with top K singular vectors.

Reconstructing Image with top K PCs

Let us use the top 5 singular vectors and reconstruct the matrix using matrix multiplication as shown above. Let us also visualize the reconstructed image.

num_components = 5
reconst_img_5 = np.matrix(U[:, :num_components]) * np.diag(s[:num_components]) * 
                np.matrix(V[:num_components, :])
plt.imshow(reconst_img_5)
plt.savefig('reconstructed_image_with_5_SVs.png',dpi=150, figsize=(8,6))

We can see that the top 5 components is not enough to reconstruct the image,

Reconstructing an Image with 5 SVs

Let us use top 50 singular vectors and see how the reconstructed image look like.

num_components = 50
reconst_img_50 = np.matrix(U[:, :num_components]) * np.diag(s[:num_components]) * np.matrix(V[:num_components, :])
plt.imshow(reconst_img_50)
plt.title('Reconstructed Image: 50 SVs', size=16)
plt.savefig('reconstructed_image_with_50_SVs.png',dpi=150, figsize=(8,6))

With top 50 singular vectors, we can see that we captured the essence of original image nicely.

Reconstructed Image with 50 SVs

The quality of reconstructed image would improve as we use more top singular vectors. Here is a compraison of the reconstructed image using different number of top components.

fig, axs = plt.subplots(2, 2,figsize=(10,10))
axs[0, 0].imshow(reconst_img_5)
axs[0, 0].set_title('Reconstructed Image: 5 SVs', size=16)
axs[0, 1].imshow(reconst_img_50)
axs[0, 1].set_title('Reconstructed Image: 50 SVs', size=16)
axs[1, 0].imshow(reconst_img_100)
axs[1, 0].set_title('Reconstructed Image: 100 SVs', size=16)
axs[1, 1].imshow(reconst_img_1000)
axs[1, 1].set_title('Reconstructed Image: 1000 SVs', size=16)
plt.tight_layout()
plt.savefig('reconstructed_images_using_different_SVs.jpg',dpi=150)

We can see the improvement in the quality of image as we add more singular vectors initally and then it kind of saturates suggesting that we don’t gain much adding more components as the variance explained is small after the top components

Reconstructed Images from SVD with different SVs

To summarize, in this post we saw how we can use SVD to decompose an image and reconstruct using top singular vectors from SVD using Python.