0

I recently learned about different methods of PCA. I decided to manually implement PCA in Python with Eigendecomposition of cov(X) and the Singular Value Decomposition of X and compare the results. I heard that the main difference is that SVD should give a more accurate result while taking longer to execute. Here is my implementation:

import numpy as np
from scipy.linalg import svd,eig
SIZE=(3,3) #the data shape

def svd_(data):
    mean=np.mean(data,axis=0)
    centered=data-mean #center data using column means
    print(f'{centered=}')
    U, s, Vh=svd(centered)
    print(f'{U=}\n{s=}\n,{Vh=}')

def eig_(data):
    mean=np.mean(data,axis=0)
    centered=data-mean #center data using column means
    print(f'{centered=}')
    cov=np.cov(centered.T)
    eigen_values, eigen_vectors=eig(cov)
    print(f'{eigen_values=}\n{eigen_vectors=}')


def main():
    np.random.seed(42) #fixed seed
    data=np.random.random(size=SIZE)
    print(f'data is\n{data}')
    svd_(data)
    print()
    eig_(data)

if __name__=='__main__':
    main()

I get the following output:

data is
[[0.37454012 0.95071431 0.73199394]
 [0.59865848 0.15601864 0.15599452]
 [0.05808361 0.86617615 0.60111501]]

#svd_ below
centered=array([[ 0.03077938,  0.29307794,  0.23562612],
       [ 0.25489775, -0.50161772, -0.3403733 ],
       [-0.28567713,  0.20853978,  0.10474719]])
U=array([[ 0.41479721,  0.7032851 ,  0.57735027],
       [-0.81646137,  0.00758237,  0.57735027],
       [ 0.40166416, -0.71086748,  0.57735027]])
s=array([8.05432409e-01, 2.49318619e-01, 2.61896384e-17])
,Vh=array([[-0.38500217,  0.76341895,  0.5186182 ],
       [ 0.90910975,  0.21687009,  0.35564986],
       [ 0.15903707,  0.60840683, -0.77752707]])  

#eig_ below
centered=array([[ 0.03077938,  0.29307794,  0.23562612],
       [ 0.25489775, -0.50161772, -0.3403733 ],
       [-0.28567713,  0.20853978,  0.10474719]])
eigen_values=array([ 3.24360683e-01+0.j,  3.10798868e-02+0.j, -4.80648600e-18+0.j])
eigen_vectors=array([[ 0.38500217, -0.90910975, -0.15903707],
       [-0.76341895, -0.21687009, -0.60840683],
       [-0.5186182 , -0.35564986,  0.77752707]])

The main problem I have is that the values for Principal Components are too different between each method. I would expect U and eigen_vectors to be close to each other. Am I missing something?

TLDR. Assuming that U, and eigen_vectors represent principal components for the original matrix X, why are they not equal?

Alex.Kh
  • 155
  • 4
  • 1
    It's not clear what you're specifically asking about. The eigenvalues are numerical. The linked thread discusses why the sign of the eigenvectors is arbitrary (if $v$ is an eigenvector to a matrix, then $-v$ is also). If you have a different question, please [edit] to explain what you know and what you would like to know. – Sycorax Jan 14 '22 at 20:36
  • @Sycorax Thanks, for the review. I've modified the question to make it easier to answer. Hope this helps now – Alex.Kh Jan 14 '22 at 20:40
  • @Sycorax That's not quite true. Notice that the squares of the singular values of the column-centered matrix exactly equal (to numerical precision) twice ($=3-1$) the eigenvalues of the covariance matrix. This relationship holds generally for rectangular matrices. The issue in the question is that *of course* the eigenvalues of the square of a matrix (which is what the covariance is proportional to) will not equal the singular values of that matrix, as is evident for almost all $1\times 1$ matrices. The svd method in Python also returns the *transpose* of the eigenvectors. – whuber Jan 14 '22 at 20:51
  • @whuber I'm not sure where we are in disagreement. – Sycorax Jan 14 '22 at 20:54
  • @Sycorax In a now deleted comment, to which I was reacting, you wrote "the svd (or eigendata) of X⊤X is not the same as the svd (or eigendata) of XX⊤ in general." The eigenvalues are *very closely* related! – whuber Jan 14 '22 at 20:56
  • 1
    Yes, I agree. I was reacting to an edit: OP wrote in a (removed) edit that they found $-V^T$ to match `eigenvectors`, which prompted my remark. If we take $X = USV^T$, one covariance eigen factorization is $\Sigma \propto X^TX=VS^2 V^T$. But if we expect the eigenvectors to be $U$, then something has gone wrong; perhaps we've transposed $X$. That's what I was trying to point out. (And likewise there's a similar discussion to be had about scaling factors...) This is why I suggest using a non-square $X$: you'll have the wrong number of vectors when you've transposed by mistake! – Sycorax Jan 14 '22 at 21:02
  • 1
    @Sycorax Thank you for the explanation. (I never saw that edit you were responding to...) – whuber Jan 14 '22 at 21:06
  • @whuber not sure why the question got closed again even after getting an answer that was deleted but I think you got the gist. If `numpy.linalg.svd` returns the transpose of the vectors, it explains why in my implementation it is so. As @sycorax mentioned, I indeed found that in my case `eigenvectors`=$-V^T$. I don't do any transpose of the original matrix $X$ myself, so I am unsure where the minus sign is coming from. Should I assume it is happening purely due to the SVD matrix not being unique? – Alex.Kh Jan 14 '22 at 21:08
  • I think you'll be able to answer your own question if you do the following: (1) write down a non-square $X$ with shape $(n,p)$; (2) label whether rows (n) or columns (p) are the features (and likewise the observations); (3) write down the algebra. There's enough ambiguity in this presentation that it is unclear what, exactly, is the sticking point for you. If you get stuck, the duplicate threads work through the answer. – Sycorax Jan 14 '22 at 21:10
  • I have just reviewed the three duplicates and believe you can find thorough explanations of your situation in all three. As @Sycorax remarked earlier, eigenvectors are merely unit *basis* vectors of eigenspaces. At *best* they are defined only up to sign. See https://stats.stackexchange.com/a/34401/919 for an example. – whuber Jan 14 '22 at 21:11
  • @whuber fair enough! Thanks for pointing in the right direction – Alex.Kh Jan 14 '22 at 21:12

0 Answers0