In [30]:
import numpy as np
import pandas as pd
from scipy.ndimage import imread
from matplotlib import pylab as plt
SVD is one of the matrix factorization tecniques. It factors a matrix into three parts with which we can reconstruct the initial matrix. However, reconstructing original matrix is not mostly the primary aim. Rather, we factorize matrices in order to achive following goals:
In a simple terms, factorization can be defined as breaking something into its building blocks, in other terms, its factors. Using SVD, we can decompose a matrix into three separate matrices as follows:
$$ A_{m x n} = U_{m x r} * \Sigma_{r x r} * (V_{n x r})^{T} $$where
In [39]:
A = np.mat([
[4, 5, 4, 1, 1],
[5, 3, 5, 0, 0],
[0, 1, 0, 1, 1],
[0, 0, 0, 0, 1],
[1, 0, 0, 4, 5],
[0, 1, 0, 5, 4],
])
In [40]:
U, S, V = np.linalg.svd(A)
U.shape, S.shape, V.shape
Out[40]:
In [44]:
U
Out[44]:
In [48]:
S
Out[48]:
In [49]:
np.diag(S)
Out[49]:
As you can see, the singular values are sorted descendingly.
In [55]:
V
Out[55]:
In [71]:
def reconstruct(U, S, V, rank):
return U[:,0:rank] * np.diag(S[:rank]) * V[:rank]
r = len(S)
reconstruct(U, S, V, r)
Out[71]:
We use all the dimentions to get back to the original matrix. As a result, we obtain the matrix which is almost identical. Let's calculate the difference between the two matrices.
In [72]:
def calcError(A, B):
return np.sum(np.power(A - B, 2))
calcError(A, reconstruct(U, S, V, r))
Out[72]:
Expectedly, the error is infinitesimal.
However, most of the time this is not our intention. Instead of using all the dimentions(rank), we only use some of them, which have more variance, in other words, which provides more information. Let's see what we will get when using only the first three most significant dimentions.
In [74]:
reconstruct(U, S, V, 3)
Out[74]:
In [76]:
calcError(A, reconstruct(U, S, V, 3))
Out[76]:
Again, the reconstructed matrix is very similar to the original one. And the total error is still small.
Now we can ask the question that which rank should we pick? There is trade-off that when you use more rank, you get closer to the original matrix and have less error, however you need to keep more data. On the other hand, if you use less rank, you will have much error but save space and remove the redundant dimentions and noise.
In [79]:
reconstruct(U, S, V, 2)
Out[79]:
In [80]:
calcError(A, reconstruct(U, S, V, 2))
Out[80]:
In [33]:
A = np.mat([
[4, 5, 4, 0, 4, 0, 0, 1, 0, 1, 2, 1],
[5, 3, 5, 5, 0, 1, 0, 0, 2, 0, 0, 2],
[0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0],
[0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1],
[1, 0, 1, 0, 1, 5, 0, 0, 4, 5, 4, 0],
[0, 1, 1, 0, 0, 4, 3, 5, 5, 3, 4, 0],
])
In [11]:
def reconstruct(U, S, V, rank):
return U[:,0:rank] * np.diag(S[:rank]) * V[:rank]
In [13]:
for rank in range(1, len(S)):
rA = reconstruct(U, S, V, rank)
error = calcError(A, rA)
coverage = S[:rank].sum() / S.sum()
print("with rank {}, coverage: {:.4f}, error: {:.4f}".format(rank, coverage, error))
As it can be seen above, more rank is used, less error occur. From another perspective, we get closer to the original data by increasing rank number.
On the other hand, after a certain rank, using more rank will not contribute as much as
Let's compare a reconstructed column to the original one by just naked eyes. Even it is reconstructed using only 4 dimention, we almost, with some error, get the original data.
In [29]:
print("Original:\n", A[:,10])
print("Reconstructed:\n", reconstruct(U, S, V, 4)[:,10])
In [185]:
imread("data/pacman.png", flatten=True).shape
Out[185]:
In [186]:
A = np.mat(imread("data/pacman.png", flatten=True))
In [187]:
U, S, V = np.linalg.svd(A)
In [188]:
A.shape, U.shape, S.shape, V.shape
Out[188]:
In [189]:
for rank in range(1, len(S)):
rA = reconstruct(U, S, V, rank)
error = calcError(A, rA)
coverage = S[:rank].sum() / S.sum()
print("with rank {}, coverage: {:.4f}, error: {:.4f}".format(rank, coverage, error))
In [192]:
for i in range(1, 50, 5):
rA = reconstruct(U, S, V, i)
print(rA.shape)
plt.imshow(rA, cmap='gray')
plt.show()
In [ ]:
plt.imshow(data, interpolation='nearest')
In [164]:
128 * 128
Out[164]:
In [165]:
- (10*128*2)
Out[165]:
In [167]:
from PIL import Image
A = np.mat(imread("data/noise.png", flatten=True))
img = Image.open('data/noise.png')
imggray = img.convert('LA')
In [171]:
imgmat = np.array(list(imggray.getdata(band=0)), float)
In [180]:
imgmat = np.array(list(imggray.getdata(band=0)), float)
imgmat.shape = (imggray.size[1], imggray.size[0])
imgmat = np.matrix(imgmat)
plt.figure(figsize=(9,6))
plt.imshow(imgmat, cmap='gray');
plt.show()
In [181]:
U, S, V = np.linalg.svd(imgmat)
In [183]:
for i in range(1, 10, 1):
rA = reconstruct(U, S, V, i)
print(rA.shape)
plt.imshow(rA, cmap='gray');
plt.show()
In [ ]: