Question

SVD is a decomposition of the data matrix \(X = U S V^T\) where \(U\) and \(V\) are orthogonal matrices and \(S\) is a diagnonal matrix.

Recall that the transpose of an orthogonal matrix is also its inverse, so if we multiply on the right by \(X^T\), we get the follwoing simplification

Compare with the eigendecomposition of a matrix \(A = W \Lambda W^{-1}\), we see that SVD gives us the eigendecomposition of the matrix \(XX^T\), which as we have just seen, is basically a scaled version of the covariance for a data matrix with zero mean, with the eigenvectors given by \(U\) and eigenvealuse by \(S^2\) (scaled by \(n-1\))..

u, s, v = np.linalg.svd(x)

e2 = s**2/(n-1)
v2 = u
plt.scatter(x[0,:], x[1,:], alpha=0.2)
for e_, v_ in zip(e2, v2):
    plt.plot([0, 3*e_*v_[0]], [0, 3*e_*v_[1]], 'r-', lw=2)
plt.axis([-3,3,-3,3]);

v1 # from eigenvectors of covariance matrix

array([[ 0.9205, -0.3909],
       [ 0.3909,  0.9205]])

v2 # from SVD

array([[-0.9205, -0.3909],
       [-0.3909,  0.9205]])

e1 # from eigenvalues of covariance matrix

array([ 0.7204,  0.1161])

e2 # from SVD

array([ 0.7204,  0.1161])

Exercises

The exercise is meant to help your understanding of what is going on when PCA is performed on a data set and how it can remove linear redundancy. You would normally use a library function to perform PCA to reduce the dimensionality of a real data set.

1 . Create a data set of 100 3-vectors such that the first componnent \(a_0 \sim N(0,1)\), \(a_1 \sim a_0 + N(0,3)\) and \(a_2 = 2a_0 + a_1\), where \(N(\mu, \sigma)\) is the normal distribution with mean \(\mu\) and standard deviaiton \(\sigma\). The data set should be a \(3 \times 100\) matrix. Normalzie so that the row means are 0.

a0 = np.random.normal(0, 1, 100)
a1 = a0 + np.random.normal(0, 3, 100)
a2 = 2*a0 + a1
A = np.row_stack([a0, a1, a2])
A = A - A.mean(0)
A.shape

(3, 100)

2 . Find the eigenvecors ane eigenvalues of the coveriance matrix of the data set using spectral decomposition.

import scipy.linalg as la

M = np.cov(A)
e, v = la.eig(M)
idx = np.argsort(e)[::-1]
e = e[idx]
e = np.real_if_close(e)
v = v[:, idx]

print M
print e
print v

[[ 4.1765 -1.4026 -2.7739]
 [-1.4026  1.3427  0.0598]
 [-2.7739  0.0598  2.7141]]
[  6.5711e+00   1.6621e+00   8.2823e-16]
[[-0.7906  0.204   0.5774]
 [ 0.2186 -0.7867  0.5774]
 [ 0.572   0.5827  0.5774]]

3 . Find the eigenvecors ane eigenvalues of the coveriance matrix of the data set using SVD. Cehck that they are equivalent to those found using spectral decomposition.

u, s, v = la.svd(A)
print s**2/(A.shape[1] - 1)
print u

[  6.5995e+00   1.6711e+00   3.4685e-31]
[[-0.7899  0.2066  0.5774]
 [ 0.2161 -0.7874  0.5774]
 [ 0.5739  0.5808  0.5774]]

4 . What percent of the total variability is explained by the principal compoennts? Given how the dataset was constructed, do these make sense? Reduce the dimenisionality of the system so that over 99% of the total variability is retained.

print "Explained variance", np.cumsum(e)/e.sum()
# note: a2 is a linear combination of a1
# and a0 explains part of the variance of a1 by construction

u.dot(np.diag(e).dot(u.T))

B1 = u.T.dot(A) # in PCA coordinates

e[2:] = 0
A1 = u.dot(np.diag(e).dot(B1)) # in original coorindate with dimension reduction

Explained variance [ 0.7981  1.      1.    ]

5 . Plot the data points in the origianla and PCA coordiantes as a set of scatter plots. Your final figure should have 2 rows of 3 plots each, where the columns show the (0,1), (0,2) and (1,2) proejctions.

plt.figure(figsize=(12, 8))
dims = [(0,1), (0,2), (1,2)]
for k, dim in enumerate(dims):
    plt.subplot(2, 3, k+1)
    plt.scatter(A[dim[0], :], A[dim[1], :])
    plt.subplot(2, 3, k+4)
    plt.scatter(B1[dim[0], :], B1[dim[1], :])

6 . Use the decomposition.PCA() function from the sklearn package to perfrom the decomposiiton.

from sklearn.decomposition import PCA

pca = PCA(copy=True)
pca.fit(A.T)
B2 = pca.transform(A.T)
B2 = B2.T

plt.figure(figsize=(12, 8))
dims = [(0,1), (0,2), (1,2)]
for k, dim in enumerate(dims):
    plt.subplot(2, 3, k+1)
    plt.scatter(A[dim[0], :], A[dim[1], :])
    plt.subplot(2, 3, k+4)
    plt.scatter(B2[dim[0], :], B2[dim[1], :])