Principal Component Analysis in Python: Analytical Mistake

Question

I'm implementing a Principal Component Analysis for face recognition in Python without making use of the already defined PCA methods in numpy or OpenCV. But my results are literally rubbish.

I read the documentation and algorithmic description at OpenCV's doc. But some things aren't clear.

• If X = {x1, x2, ..., xn} then I believe it isn't the X matrix used when calculating the Covariance matrix? But that you've to subtract the mean etc. as explained in the first two steps.
• The eigenvectors must be sorted corresponding to the eigenvalues in descending order. Do you have to sort while looking at the absolute value of the corresponding eigenvalue or not? Anyway the order is not the main issue as I plot all the eigenvectors, so I can rearrange them accordingly. I think I made an analytical mistake.

I implemented the following:

MU = X.mean(axis=0)
for i in range(n):
    X[i,:] -= MU

S = (np.dot(X, X.T) / float(n))
#Hermitian (or symmetric) matrix.
eigenvalues, eigenvectors = numpy.linalg.eigh(S)
eigenvectors = numpy.dot(X.T, eigenvectors)
for i in range(n):
    eigenvectors[:,i] = normalize_vector(eigenvectors[:,i])

Note that the sample values are stored in the rows and not the columns. So X is of shape nxd with n the number of samples and d the dimension of the samples.

The above image is the reference. The first is the mean, the following three are the largest eigenvectors. The below image is my result. The first is the mean, the following are all the eigenvectors in some order. But they seem not to match the result.

The cv2.PCACompute(X, 6) still produces better results.

Answer 1

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:

1. 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 with

   eigenvectors = eigenvectors[:,indexes][0:nb_components]

2. 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 to nb_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)