How to decide the number of primary components for PCA Background

Trying to identify the number of primary components to use (k) for PCA for MNIST aiming at 95%.

from sklearn.datasets import fetch_openml
mnist = fetch_openml('mnist_784', version=1)

# Split data into training and test
X, y = mnist["data"], mnist["target"]
X_train, y_train = X[:60000], y[:60000]

COVERAGE=0.95


If I follow Coursera Machine Learning - Principal Component Analysis Algorithm it is 67.



from sklearn.preprocessing import StandardScaler
X_centered = StandardScaler().fit_transform(X_train - X_train.mean(axis=0))
covariance_matrx = X_centered.T.dot(X_centered) 
U, s, Vt= sp.linalg.svd(covariance_matrx)

calculated_coverages = ((s ** 2) / (len(s) -1)).cumsum()
calculated_coverages = calculated_coverages / calculated_coverages[-1]
k = np.argmax(np.array(calculated_coverages) &gt;= COVERAGE)    
print("k-th component to cover {0} is {1}".format(calculated_coverages[k], k))



  k-th component to cover 0.9507022719172283 is 66


However, if I use explained_variance_ratio_ from scikit learn, it is 154.

from sklearn.decomposition import PCA
pca = PCA()
pca.fit(X_train)

contributions = pca.explained_variance_ratio_
coverages = pca.explained_variance_ratio_.cumsum()
k = np.argmax(coverages &gt;= COVERAGE)

print("k-th primary compoent for 95% coverage is {}".format(k + 1))



  k-th primary compoent for 95% coverage is 154


When I look at scikit-learn/sklearn/decomposition/_pca.py, it looks the logic is the same.

    U, S, V = linalg.svd(X, full_matrices=False)
    # flip eigenvectors' sign to enforce deterministic output
    U, V = svd_flip(U, V)

    components_ = V

    # Get variance explained by singular values
    explained_variance_ = (S ** 2) / (n_samples - 1)
    total_var = explained_variance_.sum()
    explained_variance_ratio_ = explained_variance_ / total_var
    singular_values_ = S.copy()  # Store the singular values.


Question

Please help understand why they are different.

Related


Deciding on the number of components in PCA
Python scikit learn pca.explained_variance_ratio_ cutoff
There are two ways to perform the PCA:


Compute the eigenvalue decomposition of the covariance matrix $\Sigma$
Compute the singular value decomposition of the data matrix $X$


Numerically, you can do both by calling svd() on either of them, as for positive semi-definite matrices (like $\Sigma$) svd() gives you the eigenvalue decomposition.

There is a difference, though, when it comes to interpreting the result:


The singular values in s are the eigenvalues of $\Sigma$, i.e. the variances along the PCs
The singular values in s are the singular values of $X$, i.e. the square roots of the variances along the PCs


In sklearn they go with method 2. Hence they need to square the singular values to compute the coverage. In coursera they go with method 1 so no need to square s. In the slide you show and in the video you linked they just sum the values up.

Without having run your code, my guess would be that if you change the line

calculated_coverages = ((s ** 2) / (len(s) -1)).cumsum()


to 

calculated_coverages = (s / (len(s) -1)).cumsum()


you will get better results.

Addendum: On second thought I'm not sure how StandardScaler() impacts the results of the PCA either. When comparing, make sure that it is applied in both your implementation of the PCA and the implementation provided by sklearn (and maybe leave a comment if this mattered or not, pretty please ;)).