I was trying PCA on a dataset (#samples=24, #dims=42) via eigendecomposition using numpy. I read that for matrices where the number of features exceeds the number of samples, we should use the gram matrix instead of the covariance matrix (Is PCA still done via the eigendecomposition of the covariance matrix when dimensionality is larger than the number of observations?).
My dataset matrix is $X$ (scaled and centered), and I tried the eigendecomposition once on $X^TX$, and once on $XX^T$. Next, I took the dot product of the original data matrix with eigenvectors: $X \cdot U$ to get the transformed dataset. Finally I plotted the scatter plots of the transformed data points with the different principal components (the two 5x5 subplots given below). Clearly, the plots look very different.
I want to know the reason for this difference. From what I have read so far, it is not expected to differ much. I am also attaching my code below the figures.
The figure above is from the $X^TX$ computation, and the figure below is using $XX^T$:
And here is the code...
data_stdx=StandardScaler().fit_transform(data)
cov_data=np.dot(data_stdx.T,data_stdx) # Gives a 42x42 matrix. Top Figure; doesn't appear to be correct.
#cov_data=np.dot(data_stdx,data_stdx.T) # Gives a 24x24 matrix. Bottom Figure; seems correct.
eigenvalues, eigenvectors = np.linalg.eig(cov_data)
P = eigenvectors.real#*np.sqrt(eigenvalues.real) # Not sure whether to multiply or not!
data_new=np.dot(data_stdx,P)
