# 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) >= 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 >= 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.


# Related

• As a side note - in SKLearn, PCA can recieve a number <=1 as n_components, and will keep the top k components that have this number as the ratio of explained variance (what you called "Coverage"). SKLearn docs – Itamar Mushkin Oct 30 '19 at 13:31

There are two ways to perform the PCA:

1. Compute the eigenvalue decomposition of the covariance matrix $$\Sigma$$

2. 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:

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

2. 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 ;)).

• Thanks for the folloup. By changing to "calculated_coverages = (s / (len(s) -1)).cumsum()", k becomes 330. If I remove the StandardScaler, k become 18. So both cases, it still is far away from the scikit learn result... – mon Oct 30 '19 at 6:51
• Ok, I tried it out myself: if you remove StandardScalar AND replace s**2 by s you will get k=153. If you only replace s**2 by s you get k=330 and if you only remove StandardScaler you get k=18. But it works if you do both. – matthiaw91 Oct 30 '19 at 9:35

After looking into the pca.py code of Scikit Learn, understood the original code was incorrect.

I was thinking that SVD was to decompose a co-variance into U, S, V. However apparently it is to decompose input matrix X into U, S, V with no need to create the co-variance matrix.

## Incorrect

covariance_matrx = X_centered.T.dot(X_centered)
U, s, Vt= sp.linalg.svd(covariance_matrx)


## Correct

X_centered = X_train - X_train.mean(axis=0)
U, s, Vt= sp.linalg.svd(X_centered, full_matrices=False)


Without "full_matrices=False", memory error was caused.

• You can do both, it depends on how you interpret the results. – matthiaw91 Oct 30 '19 at 11:43