2
$\begingroup$

Input matrix M x N, M rows (number of samples), N cols (dimesionality of data).

But I don't understand why we must transpose matrix when put it in np.cov?

def pca():
    #M x N
    xs= np.loadtxt("data_3d.txt",delimiter=" ", skiprows=1, usecols=(0,1,2))
    print xs.shape
    # print xs

    #get mean
    mean= np.mean(xs,axis=0)
    # print mean.shape
    # print mean

    #N x M
    data= (xs-mean).T # why need transpose?
    print data.shape
    # print data

    #N x N
    covData=np.cov(data)#calculate covariance matrix
    print covData.shape

    eigenvalues, eigenvectors = np.linalg.eig(covData)
    print eigenvalues.shape # N long
    print eigenvectors.shape # N x N
    print eigenvalues
    print eigenvectors

And how then project and reconstruct data? By projection we need to obtain matrix M x k where k< N. By reconstruction we must obtain data before we put it in np.cov , so it must be N x M.

#sort and get k largest eigenvalues
k=2
idx = eigenvalues.argsort()[-k:][::-1]
print idx

eigenvalues = eigenvalues[idx] # k long 
eigenvectors = eigenvectors[:,idx] # N x k
print eigenvalues.shape
print eigenvectors.shape
print eigenvalues
print eigenvectors

#projection and reconstruction
pr= np.dot(data.T,eigenvectors) # (M N) * (N k) => (M k)
rec= np.dot(eigenvectors, pr.T) #(N k) * (k M) => (N M)
print (data-rec) # test reconstruction error
$\endgroup$

1 Answer 1

3
$\begingroup$

If you take a look at the np.cov documentation ( http://docs.scipy.org/doc/numpy/reference/generated/numpy.cov.html ) you can see that your input matrix shall have each row corresponding to a variable and each column corresponding to a different sample, and that's exactly the transpose of your matrix.

For the reconstruction error, can you more precise about what type of error you are seeing?

But to reconstruct you should multiply the (M k) matrix with the inverse of your eigenvectors matrix.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.