First question is quite well answer here.
For the second one the answer is you need to order the eigenvectors according to the value of their corresponding eigenvalue. For that, you can use python np.argsort()
and then reverse this order (argsort order from smaller to bigger):
indexes = np.argsort(eigenvalues)[::-1]
eigval = eigenvalues[indexes]
eigvec = eigenvectors[:,indexes]
Checking your code, I've just found a couple of problems:
Now, you are getting all the eigenvectors, you are ignoring the
nb_components
parameters, you should take just the vectors you are asked for. This is made witheigenvectors = eigenvectors[:,indexes][0:nb_components]
For normalizing vectors (inside pca function), you are using a for loop from 0 to
n
, but if you are asked only for, lets say, 3 eigenvectors, you only will have 3 columns. To solve it, iterate from 0 tonb_components
.
Apart from that, your code works perfect. I tried it only taking 3 Principals components and I got 6/6 as final result. In my opinion, the difference when showing the eigenvectors is just a representation problem when casting from float to uint8 to use imshow.
And about the negative eigenvalues, it is just a matter of eigh
. As eigenvalues shows the variance in a direction, we care about absolute value but if we change a sign, we also have to change the "direcction" (eigenvector). You can make this multiplying negative eigenvalues and their corresponding eigenvectors with -1.0
(see this):
s = np.where(eigenvalues < 0)
eigenvalues[s] = eigenvalues[s] * -1.0
eigenvectors[:,s] = eigenvectors[:,s] * -1.0
You can also solve this using numpy.linalg.svd
(docs), but is should be slower than numpy.linalg.eigh
.
To sum up, this is the code I came up from yours (I deleted all comments when coping here to make it shorter):
def pca(X, nb_components=0):
[n,d] = X.shape
if (nb_components <= 0) or (nb_components>6):
nb_components = n
MU = X.mean(axis=0)
for i in range(n):
X[i,:] -= MU
S = (np.dot(X, X.T) / float(n))
eigenvalues, eigenvectors = np.linalg.eigh(S)
s = np.where(eigenvalues < 0)
eigenvalues[s] = eigenvalues[s] * -1.0
eigenvectors[:,s] = eigenvectors[:,s] * -1.0
indexes = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[indexes]
eigenvectors = eigenvectors[:,indexes][:,0:nb_components]
eigenvectors = np.dot(X.T, eigenvectors)
for i in range(nb_components):
eigenvectors[:,i] = normalize_vector(eigenvectors[:,i])
return (eigenvalues, eigenvectors, MU)