Why does eig give an inverted PC1 value compared to eigh, but PC2 is the same?

Question

I am comparing eig and eigh and I would like clarification of why the results are partly inverted!

f1 and f2 are my 2 features.

pc1 and pc2 are the principal components using eig.

pc1h and pc2h are using eigh.

As can be seen, pc2 and pc2h appear very similar.

But pc1 appears inverted to pc1h. Also pc1h is similar to f2 in shape, which I would expect.

Why is pc1 upside-down?

import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
import matplotlib.pyplot as plt

#Generate a dummy dataset.
np.random.seed(300)
X = np.random.randint(10,50,100).reshape(50,2)
X[:,1] =X[:,1]*10
def PCA(X, num_components):
    # Step-1
    X_meaned = X - np.mean(X, axis=0)

    # Step-2
    cov_mat = np.cov(X_meaned, rowvar=False)

    # Step-3
    eigen_values, eigen_vectors = np.linalg.eig(cov_mat)

    # Step-4
    sorted_index = np.argsort(eigen_values)[::-1]
    sorted_eigenvalue = eigen_values[sorted_index]
    sorted_eigenvectors = eigen_vectors[:, sorted_index]

    # Step-5
    eigenvector_subset = sorted_eigenvectors[:, 0:num_components]

    # Step-6
    X_reduced = np.dot(eigenvector_subset.transpose(), X_meaned.transpose()).transpose()

    return X_reduced

def PCAh(X, num_components):
    # Step-1
    X_meaned = X - np.mean(X, axis=0)

    # Step-2
    cov_mat = np.cov(X_meaned, rowvar=False)

    # Step-3
    eigen_values, eigen_vectors = np.linalg.eigh(cov_mat)

    # Step-4
    sorted_index = np.argsort(eigen_values)[::-1]
    sorted_eigenvalue = eigen_values[sorted_index]
    sorted_eigenvectors = eigen_vectors[:, sorted_index]

    # Step-5
    eigenvector_subset = sorted_eigenvectors[:, 0:num_components]

    # Step-6
    X_reduced = np.dot(eigenvector_subset.transpose(), X_meaned.transpose()).transpose()

    return X_reduced

X_reduced = PCA(X,2)
X_reducedh = PCAh(X,2)

df_y = pd.DataFrame()
df_y['f1'] = X[:,0]
df_y['f2'] = X[:,1]
df_y['pc1'] = X_reduced[:,0]
df_y['pc2'] = X_reduced[:,1]
df_y['pc1h'] = X_reducedh[:,0]
df_y['pc2h'] = X_reducedh[:,1]
df_y.plot()
plt.show(block=True)

Answer 1

I don't see any other answer to this "why?" than "why not"

Understand that there is always an infinite number of eigenvectors for a given eigenvalues. So the choice of which eigenvector is arbitrary.

For matrix M

[[-1/2,  3/2],
 [ 3/2, -1/2]]

Eigenvalues are λ₁=-2 and λ₂=-2. And example of eigenvectors are u₁=[1,-1] and u₂=[1,1].

You can check that M@u₂ = M@[1,1] = [-1/2+3/2, 3/2-1/2] = [1,1] = λ₂u₂.
And M@u₁ = M@[1,-1] = [-1/2-3/2, 3/2+1/2 ] = [-2, 2] = -2[1,-1] = λ₁u₁.

So, sure, since you have 2 different eigenvalues in a 2D space, you know that each eigenspace is dimension 1. So, all eigenvectors associated to the same eigenvalue are colinear.

But still, that let room for an infinite number of them. For example, v₁=[-2,2] and v₂=[-3,-3] are also eigenvectors. M@[-3,-3] = [3/2-9/2, -9/2+3/2] = λ₂[-3,-3]
and M@[-2,2] = [1+3, -3-1] = [4,-4] = λ₁[-2,2]

So if you have a data X=[1,2] in the canonical space, the representation of that X in the system with basis u₁,u₂ = [1,-1], [1,1] is

X_reduced_v1 = [-0.5, 1.5]
# meaning that X is -0.5[1,-1]+1.5[1,1], which is indeed [1,2]

If you represent the same data X=[1,2] in another system with basis v₁,v₂=[-2,2],[-3,-3], then you will have

X_reduced_v2 = [1/4, -1/2]
# Meaning that X is 1/4[-2,2]-1/2[-3,-3], which is indeed [1, 2]

So, both reduction are correct. And both are done in a basis made of eigenvectors of my matrix. So both represent equally well X in terms of eigenvectors.

Now, np.eig and np.eigh both returns normalized vector. So since in my example eigenspace are dimension 1, that reduce a lot the infinite number of possibilities. Which is convenient in your case, because you are obviously doing this for PCA. But even then, there is surely no longer an infinite number of possible eigenvectors, but there are still 2 possibility for each eigenvectors.

In my example eigenvectors with norm 1 could be U₁=[1/√2, -1/√2] for the one associated with λ₁=-2, and U₂=[1/√2, 1/√2] for the one associated with λ₂=1. And in such system, my vector [1,2] would have coordinates

X_reduced_v3 = [-0.70710678,  2.12132034])
# You can check that indeed, -0.70710678*U₁+2.21132034*U₂ = [1,2]

But V₁=[-1/√2, 1/√2] and V₂=[-1/√2, -1/√2] are also eigenvectors of norm 1 associated with eigenvalues λ₁=-2 and λ₂=1.

In system V₁,V₂, X coordinates become

X_reduced_v4 = [ 0.70710678, -2.12132034]

There is absolutely no criterion left to choose between those 2 systems (nor the 2 systems I haven't mentioned, made of vectors [1/√2, -1/√2], [-1/√2, -1/√2], and vectors [-1/√2, -1/√2], [1/√2, 1/√2]. So 4 different possible basis). Those 4 basis are made of eigenvectors of my matrix. Those 4 basis are made of vector of norm 1. Those 4 basis are sorted by their eigenvalues. There is no other criterion, not even an arbitrary one, to choose among those 4.

So, it really doesn't matter which one you, or the diagonalization algorithm you use, choose. The only thing that matters is that you are consistent. If you choose basis U₁,U₂, then you must say that X_reduced is [-0.70710678, 2.12132034] (that is what you call [pca₁, pca₂]).

If you choose basis V₁, V₂, then you must say that X_reduced is [0.70710678, -2.12132034]

So that, what ever basis (eigenvectors) W₁,W₂ you choose, X=pca₁W₁+pca₂W₂ remains true. There are multiple, equally valid, W₁,W₂,pca₁,pca₂ that do that. Which one is chosen is totally arbitray.

In your case eig happens to choose one, while eigh choose another. That's all.