The goal of this notebook is to understand the Singular Value Decomposition algorithm (SVD) a little bit better. It is based on a very nice discussion by Andrew Gibiansky: http://andrew.gibiansky.com/blog/mathematics/cool-linear-algebra-singular-value-decomposition/.
Instead of going into the details that Andrew goes in, we will concentrate on the two applications of SVD he demonstrates:
First, let's open the image of the tiger featured in the blog post by Andrew.
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
import urllib
from scipy.linalg import svd
from skimage.color import rgb2gray
url = 'http://andrew.gibiansky.com/blog/mathematics/cool-linear-algebra-singular-value-decomposition/images/tiger.jpg'
tiger = plt.imread(urllib.request.urlopen(url), format='jpeg')
plt.imshow(tiger, cmap='Greys_r')
<matplotlib.image.AxesImage at 0x82c7198>
Let's check the shape of the tiger.
tiger.shape
(1000, 1600, 3)
As in the blog post, let's convert the tiger to black and white.
tiger = rgb2gray(tiger)
from skimage.color import rgb2gray
tiger.shape
(1000, 1600)
plt.imshow(tiger, cmap='gray')
<matplotlib.image.AxesImage at 0x88128d0>
Now, let's perform the SVD.
U, s, Vh = svd(tiger)
U.shape
(1000, 1000)
Vh.shape
(1600, 1600)
s.shape
(1000,)
Let's look at the magnitude of the eigenvalues.
plt.semilogy(s, label='eigenvalue magnitude')
plt.legend()
<matplotlib.legend.Legend at 0x8821b38>
Let's look at the cumulated sum of the singular values, which can be interpreted as the (unnormed) sum of information contained in these eigenvalues.
plt.plot(np.cumsum(s), label='cumulative sum of eigenvalues')
plt.legend()
<matplotlib.legend.Legend at 0x88c6cf8>
Finally, let's look at approximations of the original image obtained by keeping a certain number of eigenvalues (and setting zeros for the rest). We can use diagsvd
as a helper for this.
from scipy.linalg import diagsvd
help(diagsvd)
Help on function diagsvd in module scipy.linalg.decomp_svd: diagsvd(s, M, N) Construct the sigma matrix in SVD from singular values and size M, N. Parameters ---------- s : (M,) or (N,) array_like Singular values M : int Size of the matrix whose singular values are `s`. N : int Size of the matrix whose singular values are `s`. Returns ------- S : (M, N) ndarray The S-matrix in the singular value decomposition
sigma = diagsvd(s, *tiger.shape)
sigma.shape
(1000, 1600)
reconstructed_tiger = U.dot(sigma.dot(Vh))
plt.figure(figsize=(10, 4))
plt.subplot(121)
plt.imshow(tiger, cmap='gray')
plt.subplot(122)
plt.imshow(reconstructed_tiger, cmap='gray')
<matplotlib.image.AxesImage at 0x94ff588>
np.allclose(tiger, reconstructed_tiger)
True
Let's now introduce some degree of approximation by reducing the number of singular values in the reconstruction.
from ipywidgets import interact
@interact
def approximate_tiger(n=(0, s.size - 1)):
"""Uses n singular values for approximating the tiger image."""
sigma = diagsvd(s * (np.arange(s.size) < n).astype(np.int), *tiger.shape)
reconstructed_tiger = U.dot(sigma.dot(Vh))
plt.figure(figsize=(10, 4))
plt.subplot(121)
plt.imshow(tiger, cmap='gray')
plt.subplot(122)
plt.imshow(reconstructed_tiger, cmap='gray')
It's quite interesting to see how much information is contained only in a few singular values.
approximate_tiger(1)
approximate_tiger(2)
approximate_tiger(4)
approximate_tiger(10)
approximate_tiger(25)
approximate_tiger(50)
Let's now move to the second application listed by Andrew, linear regression. Let's generate the test data.
intercept = -10
slope = 3
npts = 50
noise = 80
xs = 10 + np.random.random((npts, 1)) * 90;
ys = slope * xs + intercept + np.random.random((npts, 1)) * noise;
plt.plot(xs, ys, 'o', label='points to fit')
plt.legend()
<matplotlib.legend.Legend at 0xda939b0>
Let's build the $A$ matrix referred to in the post.
A = np.c_[xs, ys, - np.ones_like(xs)]
A.shape
(50, 3)
Now, let's perform the SVD:
U, s, Vh = svd(A)
s
array([ 1471.92586099, 49.29409757, 2.76407261])
We have three singular values, but only the smallest one contains our linear regression information. We can thus read the coefficients directly from the last value of Vh
:
Vh.shape
(3, 3)
fit = Vh[:, 1]
fit
array([-0.95945345, 0.28177047, 0.00738112])
plt.plot(xs, -fit[0] / fit[1] * xs + fit[2] / fit[1], 'o')
plt.plot(xs, ys, 'o')
[<matplotlib.lines.Line2D at 0xda38198>]
Interestingly, the right proper value for this linear regression is not the one described in Andrew's article. I wonder there this comes from.