Singular Value Decomposition (SVD)

Prelimenary

Let bold-face lower-case letters (like $\mathbf{a}$) refer to vectors, bold-face capital letters $\mathbf{A}$ refer to matrices, and italic lower-case letters (like $\mathcal{a}$) refer to scalars. Think $\mathbf{A}$ is a square or rectangular. Its rank is $p$. We will diagonalize this $A$, but not by $X^{-1} AX$. The eigenvectors in $X$ have three big problems:

They are usually not orthogonal
There are not always enough eigenvectors, and
$Ax = \lambda x$ requires $A$ to be a square matrix

The singular vectors of $A$ solve all those three problems in a perfect way.

Defenition

The SVD of a matrix is a sort of change of coordinates that makes the matrix simple, a generalization of diagonalization. It has some interesting algebraic properties and conveys important geometrical and theoritical insights about linear transformations.
The singular value decomposition (SVD) of a matrix is a factorization of that matrix into three matrices, where the factorization has the form $\mathbf{U\Sigma V}$.
A factorization is called the eigendecomposition of A, also called the spectral decomposition of $\mathbf{A}$.
$\mathbf{U}$ is an $m \times p$ matrix, $\mathbf{\Sigma}$ is a $p\times p$ diagonal matrix, and $\mathbf{V}$ is an $n \times p$ matrix, with $\mathbf{V}^T$ being the transpose of $V$, a $p \times n$ matrix. The value $p$ is called the rank.
The diagonal entries of $\mathbf{\Sigma}$ are referred to as the singular values of $\mathbf{A}$. The columns of $\mathbf{U}$ are typically called the left-singular vectors of $\mathbf{A}$, and the columns of $V$ are called the right-singular vectors of $\mathbf{A}$.
Unfortunately not all matrices can be diagonalized. SVD is a way to do something like diagonalization for any matrix, even non-square matrices.

Application

In many applications, the data matrix $A$ is close to a matrix of low rank and it is useful to find a low rank matrix which is a good approximation to the data matrix . We will show that from the singular value decomposition of $A$, we can get the matrix $B$ of rank $p$ which best approximates $A$.
Principal component analysis (PCA)
Clustering a mixture of spherical Gaussians
To solve discrete optimization problem
For spectral decomposition
It can be used to reduce the dimensionality , i.e., the number of columns, of a dataset.
It is often used in digital signal processing for noise reduction, image compression, and other areas.

Implementation

In this code, I will try to calculate the SVD using Numpy and Scipy. I will calculate SVD, and also perform pseudo-inverse. In the end, I will apply SVD for compressing the image.

# Imports

from skimage.color import rgb2gray
from skimage import data
import matplotlib.pyplot as plt
import numpy as np
from numpy import array
from scipy.linalg import svd

"""
Singular Value Decomposition (SVD)
"""

# define a matrix
A = array([[1, 2], [3, 4], [5, 6], [7, 8]])
print(A)

Perform SVD

U, s, VT = svd(A)

Print different components

print("U: ", U)
print("Singular  array ", s)
print("V^{T}: ", VT)

Calculate pseudo-inverse

from numpy import zeros
from numpy import diag
from numpy import dot

# create the reciprocal of s
d = 1.0/s
# create m x n matrix of zeroes 
D = zeros(A.shape)
# populate D with n x n diagonal matrix 
D[:A.shape[1], :A.shape[1]]=diag(d)

print("Sigma is ", D)

# calculate pseudo-inverse
B = VT.T.dot(D.T).dot(U.T)
print("The pseudo-inverse of matrix A is ", B)

SVD on Image Compression

cat = data.chelsea()
plt.imshow(cat)

Convert image to grayscale

gray_cat = rgb2gray(cat)
plt.imshow(gray_cat)

Calculate the SVD and plot the image

U, s, VT = svd(gray_cat, full_matrices=False)
s = np.diag(s)
print("U: ", U)
print("Singular  array ", s)
print("V^{T}: ", VT)
fig, ax  = plt.subplots(5, 2, figsize=(8,20))

curr_fig = 0
for p in [5, 10, 70, 100, 200]:
    cat_approx = U[:, :p] @ s[0:p, :p] @ VT[:p, :]
    ax[curr_fig][0].imshow(cat_approx, cmap='gray')
    ax[curr_fig][0].set_title("k = "+str(p))
    ax[curr_fig, 0].axis('off')
    ax[curr_fig][1].set_title("Original Image")
    ax[curr_fig][1].imshow(gray_cat, cmap='gray')
    ax[curr_fig, 1].axis('off')
    curr_fig += 1
plt.show()