Difference between scikit-learn implementations of PCA and TruncatedSVD

Question

I understand the relation between Principal Component Analysis and Singular Value Decomposition at an algebraic/exact level. My question is about the scikit-learn implementation.

The documentation says: "[TruncatedSVD] is very similar to PCA, but operates on sample vectors directly, instead of on a covariance matrix.", which would reflect the algebraic difference between both approaches. However, it later says: "This estimator [TruncatedSVD] supports two algorithm: a fast randomized SVD solver, and a “naive” algorithm that uses ARPACK as an eigensolver on (X * X.T) or (X.T * X), whichever is more efficient.". Regarding PCA, it says: "Linear dimensionality reduction using Singular Value Decomposition of the data to project it ...". And PCA implementation supports the same two algorithms (randomized and ARPACK) solvers plus another one, LAPACK. Looking into the code I can see that both ARPACK and LAPACK in both PCA and TruncatedSVD do svd on sample data X, ARPACK being able to deal with sparse matrices (using svds).

So, aside from different attributes and methods and that PCA can additionally do exact full singular value decomposition using LAPACK, PCA and TruncatedSVD scikit-learn implementations seem to be exactly the same algorithm. First question: Is this correct?

Second question: even though LAPACK and ARPACK use scipy.linalg.svd(X) and scipy.linalg.svds(X), being X the sample matrix, they compute the singular value decomposition or eigen-decomposition of $X^T*X$ or $X*X^T$ internally. While the "randomized" solver doesn't need to compute the product. (This is relevant in connection with numerical stability, see Why PCA of data by means of SVD of the data?). Is this correct?

Relevant code: PCA line 415. TruncatedSVD line 137.

Answer 1

PCA and TruncatedSVD scikit-learn implementations seem to be exactly the same algorithm.

No: PCA is (truncated) SVD on centered data (by per-feature mean substraction). If the data is already centered, those two classes will do the same.

In practice TruncatedSVD is useful on large sparse datasets which cannot be centered without making the memory usage explode.

• numpy.linalg.svd and scipy.linalg.svd both rely on LAPACK _GESDD described here: http://www.netlib.org/lapack/lug/node32.html (divide and conquer driver)

• scipy.sparse.linalg.svds relies on ARPACK to do a eigen value decomposition of XT . X or X . X.T (depending on the shape of the data) via the Arnoldi iteration method. The HTML user guide of ARPACK has a broken formatting which hides the computational details but the Arnoldi iteration is well described on wikipedia: https://en.wikipedia.org/wiki/Arnoldi_iteration

Here is the code for the ARPACK-based SVD in scipy:

https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/eigen/arpack/arpack.py#L1642 (search for the string for "def svds" in case of line change in the source code).

Answer 2

There is also a difference in how attribute explained_variance_ is calculated.

Let the data matrix $\mathbf X$ be of $n \times p$ size, where $n$ is the number of samples and $p$ is the number of variables. And $\mathbf{X}_c$ is the centered data matrix, i.e. column means have been subtracted and are now equal to zero in this matrix. Assume that we are reducing the dimensionality of the data from $p$ to $k \lt p$.

Then for sklearn.decomposition.PCA we have the following expressions: $$\mathbf{X}_c \approx \mathbf{U}_k \mathbf{S}_k \mathbf{V}_k^T \qquad (\text{truncated SVD of } \mathbf{X}_c);$$ $$\mathbf{L}_k = \frac{1}{n-1} \mathbf{S}^2_k \quad \Longleftrightarrow \quad \lambda_j = \frac{s_j^2}{n-1}, \quad \forall j =1,\ldots,k; \qquad(*)$$

where $\mathbf{L}_k = \mathrm{diag}(\lambda_1, \ldots, \lambda_k)$ is the matrix of $k$ largest eigenvalues of the covariance matrix $\mathbf{C} = \frac{1}{n-1} \mathbf{X}_c^T\mathbf{X}_c$, and $\mathbf{S}_k = \mathrm{diag}(s_1, \ldots, s_k)$ is the matrix of $k$ largest singular values of $\mathbf{X}_c$. You can simply prove $(*)$ if you substitute truncated SVD of $\mathbf{X}_c$ in the expression for the covariance matrix $\mathbf{C}$ and compare the result with the truncated eigendecomposition $\mathbf{C} \approx \mathbf{V}_k \mathbf{L} \mathbf{V}_k^T$ (here it was done for full decompositions). Matrix $\mathbf{L}_k$ is called explained_variance_ attribute in sklearn.decomposition.PCA.

But for sklearn.decomposition.TruncatedSVD only the following holds: $$\mathbf{X} \approx \tilde{\mathbf{U}}_k \tilde{\mathbf{S}}_k \tilde{\mathbf{V}}_k^T \qquad (\text{truncated SVD of } \mathbf{X}).$$ In this case we can't get simple equiation like $(*)$ that will link $\mathbf{L}_k$ and $\tilde{\mathbf{S}}_k$, because substitution of truncated SVD of $\mathbf{X}$ in the expression $\mathbf{C} = \frac{1}{n-1} \mathbf{X}_c^T\mathbf{X}_c = \frac{1}{n-1}\mathbf{X}^T\mathbf{X} - \frac{n}{n-1}\bar{x}\bar{x}^T$ will not be very useful. So explained_variance_ in sklearn.decomposition.TruncatedSVD is calculated instead by np.var(X_transformed, axis=0), where X_transformed = $\mathbf{X} \tilde{\mathbf{V}}_k$ – matrix of scores (new features).