Singular Value Decomposition: SVD

$\begin{matrix} \text{The matrix} \\ \text{A} \end{matrix}$ $\begin{matrix} \text{Orthogonal basis for} \\ \text{column space of A} \end{matrix}$ $\begin{matrix} \text{Singular} \\ \text{values of A} \end{matrix}$ $\begin{matrix} \text{Orthogonal basis for} \\ \text{row space of A} \end{matrix}$

$\: \: \: \left[ \begin{matrix} \\ \\ A \\ m\times n \\ \\ \end{matrix} \right] \: \: = \: \:$ $\left[ \begin{matrix} \\ \\ \quad \quad U \quad \quad \\ m\times m \\ \\ \end{matrix} \right] \quad$ $\left[ \begin{matrix} \\ \\ \sum \\ m\times n \\ \\ \end{matrix} \right] \quad \quad$ $\left[ \begin{matrix} \\ \quad V^T \quad \\ n\times n \\\end{matrix} \right] $

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
In [2]:
ein = Image.open('einstein.jpg')
print(np.shape(ein))

ein = np.mean(ein,2) # 2: 3rd dimention: 3 colors (rgb)
print(np.shape(ein))

plt.imshow(ein)
plt.show()

(1312, 1024, 3)
(1312, 1024)
In [8]:
U,s,V = np.linalg.svd(ein) # U: matrix, s: singular value, V: matrix

print(np.shape(U)) # square matrix
print(np.shape(s)) # vector
print(np.shape(V)) # square matrix

fig,ax = plt.subplots(1,3)

ax[0].imshow(U)
ax[0].set_title('U')

ax[1].imshow(np.diag(np.log(s)))
ax[1].set_title('$\\Sigma$')
print(s)
print(np.log(s))

ax[2].imshow(V)
ax[2].set_title('V')

plt.show()

(1312, 1312)
(1024,)
(1024, 1024)
[1.23753227e+05 2.54921703e+04 1.99859941e+04 ... 1.93076577e+00
 1.90566344e+00 1.86651116e+00]
[11.72604476 10.14612664  9.90278701 ...  0.6579167   0.64483021
  0.624071  ]
In [4]:
plt.plot(s,'ks-',markerfacecolor='w')
plt.xlim([-1,50])
plt.title('Eigenspectrum of Einstein')
plt.show()

Exercise 1

Reconstruct the image from the SVD.

In [5]:
S = np.zeros(np.shape(ein))

for i in range(0,len(s)):
    S[i,i] = s[i] # diagonalize
    
rein = U@S@V

plt.subplot(1,2,1)
plt.imshow(ein)
plt.axis('off')
plt.title('original')

plt.subplot(1,2,2)
plt.imshow(rein)
plt.axis('off')
plt.title('SVD version')

plt.show()

Exercise 2

Reconstruct the image but shuffle the singular values.

In [6]:
# shuffle the order of the singular values
randorder = np.random.permutation(len(s))
print(np.random.permutation(5))
print(randorder)

S = np.zeros(np.shape(ein))
for i in range(0,len(s)):
    S[i,i] = s[randorder[i]] # diagonalize
    
rein = U@S@V

plt.subplot(1,2,1)
plt.imshow(ein)
plt.axis('off')
plt.title('original')

plt.subplot(1,2,2)
plt.imshow(rein)
plt.axis('off')
plt.title('SVD version')

plt.show()

[0 3 4 1 2]
[393 519 728 ... 682 644 194]

Exercise 3

Partially reconstruct the image using the first N components.

In [7]:
S = np.zeros(np.shape(ein))
for i in range(10,len(s)):
    S[i,i] = s[i]
    
rein = U@S@V

plt.subplot(1,2,1)
plt.imshow(ein)
plt.axis('off')
plt.title('original')

plt.subplot(1,2,2)
plt.imshow(rein)
plt.axis('off')
plt.title('SVD version')

plt.show()