# PCA in numpy and sklearn produces different results [closed]

Am i misunderstanding something. This is my code

using sklearn

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sklearn import decomposition
from sklearn import datasets
from sklearn.preprocessing import StandardScaler

pca = decomposition.PCA(n_components=3)

x = np.array([
[0.387,4878, 5.42],
[0.723,12104,5.25],
[1,12756,5.52],
[1.524,6787,3.94],
])
pca.fit_transform(x)


Output:

array([[ -4.25324997e+03,  -8.41288672e-01,  -8.37858943e-03],
[  2.97275001e+03,  -1.25977271e-01,   1.82476780e-01],
[  3.62475003e+03,  -1.56843494e-01,  -1.65224286e-01],
[ -2.34425007e+03,   1.12410944e+00,  -8.87390454e-03]])


Using numpy methods

x_std = StandardScaler().fit_transform(x)
cov = np.cov(x_std.T)
ev , eig = np.linalg.eig(cov)
a = eig.dot(x_std.T)


Output

array([[ 0.06406894,  0.94063993, -1.62373172],
[-0.35357757,  0.7509653 ,  0.63365168],
[ 0.29312477,  0.6710958 ,  1.11766206],
[-0.00361615, -2.36270102, -0.12758202]])


I have kept all 3 components but it doesnt seem to allow me to retain my original data.

May I know why is it so?

If I want to obtain back my original matrix what should I do?

• Your numpy code is wrong IMHO (also, it uses X which is not defined). Recheck your math. Sep 20, 2016 at 6:12
• I am using ipython notebook so i could only copy by cells. My maths is wrong? Which part @Anony-Mousse Sep 20, 2016 at 8:41
• @Anony-Mousse Yes i realised my error but it still doesnt match up Sep 20, 2016 at 8:50
• @aceminer I am curious why you calculate covariance matrix of x_std.T, not x_std? Mar 3, 2019 at 18:35
• @EvgeniNabokov its been too long. Sry i cannot remember already Mar 6, 2019 at 6:58

The difference is because decomposition.PCA does not standardize your variables before doing PCA, whereas in your manual computation you call StandardScaler to do the standardization. Hence, you are observing this difference: PCA on correlation or covariance?

If you replace

pca.fit_transform(x)


with

x_std = StandardScaler().fit_transform(x)
pca.fit_transform(x_std)


you will get the same result as with manual computation...

...but only up to the order of the PCs. That is because when you run

ev , eig = np.linalg.eig(cov)


you get eigenvalues not necessarily in the decreasing order. I get

array([ 0.07168571,  2.49382602,  1.43448827])


So you will want to order them manually. Sklearn does that for you.

Regarding reconstructing original variables, please see How to reverse PCA and reconstruct original variables from several principal components?

• Just would like to check. Is it really necessary to standardize the matrix by its standard deviation? I seen examples where they do not do it Sep 20, 2016 at 13:50
• It is not necessary, it is just one way of doing it. See the link that I put in the first paragraph: stats.stackexchange.com/questions/53 -- it really is all about this question. If you standardize, you do PCA on correlations. If you do not, you do PCA on covariances. Sep 20, 2016 at 14:07

Here is a nice implementation with discussion and explanation of PCA in python. This implementation leads to the same result as the scikit PCA. This is another indicator that your PCA is wrong.

import numpy as np
from scipy import linalg as LA

x = np.array([
[0.387,4878, 5.42],
[0.723,12104,5.25],
[1,12756,5.52],
[1.524,6787,3.94],
])

#centering the data
x -= np.mean(x, axis = 0)

cov = np.cov(x, rowvar = False)

evals , evecs = LA.eigh(cov)


you need to sort the eigenvalues (and eigenvectors accordingly) descending

idx = np.argsort(evals)[::-1]
evecs = evecs[:,idx]
evals = evals[idx]

a = np.dot(x, evecs)


Generally, I recommend you to you check your code by implementation a simple example (as simple as possible) and calculating by hand the correct resuls (and intermediate results). This helps you to identify the problem.

• Love this answer. It solved my problem! Jul 12, 2018 at 15:28
• Note that linalg.eigh assumes the input is symmetric (that is, it is optimized to read only those values from one "triangle" half of the matrix). The OP uses linalg.eig instead, which makes no such assumption. Feb 28 at 22:12