Repository navigation
2 PCA implementations that give same results but different from Python scikit-learn implementation ... #81
Description
We have 2 implementations of PCA :
- one based on Jacobi Transformation of the covariance matrix
- another one based on SVD.
They give the same results but the result are different from the one you can find with sci-kit learn in Python:
import numpy as np
from sklearn.decomposition import PCA
X = np.array([[-1, -1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])
pca = PCA(n_components=2)
pca.fit(X)
pca.components_
pca.transform(X)pca.components_ returns :
array([[-0.83849224, -0.54491354],
[ 0.54491354, -0.83849224]])pca.transform(X) returns:
array([[ 1.38340578, 0.2935787 ],
[ 2.22189802, -0.25133484],
[ 3.6053038 , 0.04224385],
[-1.38340578, -0.2935787 ],
[-2.22189802, 0.25133484],
[-3.6053038 , -0.04224385]])I try to implement a flipsvd method like this one : https://github.com/scikit-learn/scikit-learn/blob/4c65d8e615c9331d37cbb6225c5b67c445a5c959/sklearn/utils/extmath.py#L609
but fails until now.
Please have a look to tests of PMPrincipalComponentAnalyserTest.
Activity
hemalvarambhia commentedon Nov 23, 2018
Maybe this is an odd suggestion, but to debug this can we write a simpler test, say, with the identity matrix as
X? I wrote this test in Pharo and the output is:( -1/sqrt(2), 1/sqrt(2) )
( 1/sqrt(2), 1/sqrt(2) )
which, interestingly, with
numpythe value ofpca.components:[[-0.70710678 0.70710678]
[ 0.70710678 0.70710678]]
hemalvarambhia commentedon Jan 16, 2019
Hi, @SergeStinckwich I looked into the failing test for the Jacobi Transformation form of PCA a little deeper and noted that transposing and negating the matrix of eigenvectors results in the test passing:
The Scikit Learn library, from what I can see, uses the SVD implementation which involves flipping signs. For this test, then, does it make sense to directly compare the Jacobi matrix of eigenvectors with their SVD computed one? Might we instead assert that the negated transpose of Jacobi equals the SVD one?
Please forgive my ignorance- I'm not too familiar with this domain, so I'm sure I'm missing something!
hemalvarambhia commentedon Jan 16, 2019
@SergeStinckwich I noted a similar thing with PolyMath's PCA implementation using SVD, where transposing and negating the returned matrix from
transformMatrixalso ets the test passing:This phenomenon is interesting, I think, perhaps worth investigating.
Also I note that
flipEigenvectorsSignhas been commented out in thefit:message because it doesn't work. Could you elaborate on this a little more, please?SergeStinckwich commentedon Jan 17, 2019
Thank you for spending time to investigate more on this. I forget a little bit about the details of implementation. I will have to spend some time to remember what I have done here.
hemalvarambhia commentedon Jan 17, 2019
Thank you for your help on this, @SergeStinckwich. The code is really well-written compared to the Python library. I read this paper on the subject and your work was very easy to follow from that. So far the difference between your work and the Python one is just
transpose negate.hemalvarambhia commentedon Mar 11, 2019
Hi, @SergeStinckwich,
I made some more progress on this, where I extracted out the algorithm from SciKit-Learn and tried to determine the outputs of each step. Here is something you can run on the command line to see this at a more granular level: scikit-learn flip algorithm
Given that in SciKit-Learn,
Uis(the first two columns from PolyMath match SciKit-Learn's)
the first thing to note is the
max_abs_columnvariable, which in SciKitLearn is a 2-element array:[2 0], whereas in PolyMath it is#(3 4 5...). For the former, this leads to a sub matrix:[-0.57221483 0.53396977], which then leads to the signs being[-1 1]. This is the key difference because in PolyMath, we get#(-1-1 -1 1 1 -1). Our ourargMaxOnColumnsmessage yields-0.53396...from row 4, column 2 rather than 0.53396977 at row 1.hemalvarambhia commentedon Mar 12, 2019
On further investigation, when we compute the
maxAbsColumnin PolyMath, we get the element at row 4 as the highest value because it is 0.5339697737377911, whereas at row 1 it is 0.5339697737377903. With SciKitLearn the columns aren't as precise as you know, so perhaps we need to do the same and round off.SergeStinckwich commentedon Mar 27, 2019
In Mathematica if we use SVD:
SergeStinckwich commentedon Mar 27, 2019
Action points now:
hemalvarambhia commentedon Mar 27, 2019
I tried Wolfram Alpha and our V matrix matches theirs.
SergeStinckwich commentedon Mar 27, 2019
Why Mathematica SVD is different from another math engine
Mathematica SVD documentation: https://reference.wolfram.com/language/ref/SingularValueDecomposition.html
hemalvarambhia commentedon Mar 27, 2019
Wolfram Alpha:
hemalvarambhia commentedon Mar 27, 2019
Polymath calculation: here we see it matches Wolfram Alpha exactly:
hemalvarambhia commentedon Mar 27, 2019
The next example we could try is in section 3.1 of this article
Here is the answer according to Wolfram Alpha. You can see that the V matrix is similar in that the elements have the same sign and magnitude.
hemalvarambhia commentedon Apr 14, 2019
@SergeStinckwich In relation to the first of our action points, to find another example as an acceptance test, I tried the example in section 3.1 of the PCA tutorial. The PCA of scikit learn matches the output published in the above paper almost exactly:.
This is not the case for PolyMath, sadly, for both implementations. However, as we've discovered, the problem is a bit further up, in how we're computing the eigenvectors.
hemalvarambhia commentedon Apr 19, 2019
@SergeStinckwich Out of curiosity, I used the mean centred data from the PCA tutorial:
and for the SVD-based implementation of PCA, PolyMath's output is correct up to negation:
Similarly for the Jacobi implementation, again our output is correct up to the negation:
hemalvarambhia commentedon Apr 20, 2019
On point 3, Scikit-Learn uses scipy for the SVD part of the PCA. In turn scipy delegates to LAPACK driver routines and according to the documentation, by default, it uses
gessd(the available options aregessdandgesvd). I could not see anything in the scipy source code (decomp_svd.py) that indicated the decision taken in flipping signs.