24
$\begingroup$

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.

$\endgroup$
7
  • 1
    $\begingroup$ could you add a link to the code $\endgroup$
    – seanv507
    Oct 10, 2016 at 21:49
  • 1
    $\begingroup$ drake - I think I agree with you on first Q. don't understand the second. what do you mean 'they compute the singular value decomposition or eigen-decomposition of XT∗XXT∗X or X∗XTX∗XT internally' - you have just shown the code where it is all done using SVD on X? - numerical issues refer to first computing covariance matrix (call it C) then finding eigenvectors of C $\endgroup$
    – seanv507
    Oct 11, 2016 at 0:24
  • $\begingroup$ @seanv507 Regarding 2nd question -- I guess that scipy.linalg.svd(X) computes the svd by doing the eigen-decomposition of $X^T * X$ or/and $X * X^T$. Same thing for linalg.svds(X). Quoting: "a fast randomized SVD solver, and a “naive” algorithm that uses ARPACK as an eigensolver on (X * X.T) or (X.T * X)". See also last line in docs.scipy.org/doc/scipy/reference/generated/… . The only way I can understand the first quote is that the randomized algorithm is the only one that doesn't compute the covariance/gram matrix $\endgroup$
    – drake
    Oct 11, 2016 at 0:49
  • 1
    $\begingroup$ I would guess the ARPACK approach has to do with something like Arnoldi iteration, so it only has to do matrix-vector products. (In principle these sort of iterative methods do not even an explicit $X$, just a pair of routines Xtimes() and Xt_times(). This is common for large sparse matrices in PDE solvers, for example.) $\endgroup$
    – GeoMatt22
    Oct 11, 2016 at 1:01
  • $\begingroup$ @GeoMatt22 Can you elaborate on your comment? Do you mean that ARPACK or LAPACK approaches don't suffer from numerical instabilities because they don't need to compute the covariance matrix? $\endgroup$
    – drake
    Oct 11, 2016 at 1:10

2 Answers 2

25
$\begingroup$

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

$\endgroup$
13
  • 4
    $\begingroup$ One can support sparse data efficiently (TruncatedSVD), the other cannot (PCA). This is why we have 2 classes. $\endgroup$
    – ogrisel
    Oct 12, 2016 at 14:43
  • 1
    $\begingroup$ If that's the reason, then I'd call them SVD and SparseSVD (or similar) to avoid confusion. $\endgroup$
    – drake
    Oct 12, 2016 at 22:37
  • 2
    $\begingroup$ But people want PCA and they might not know that PCA is just SVD on centered data. $\endgroup$
    – ogrisel
    Oct 13, 2016 at 7:51
  • 6
    $\begingroup$ @drake I disagree that "the procedures are different (PCA uses covariance matrix and SVD uses data matrix)". PCA is a name for the type of analysis. One can use different algorithms and implementations to carry it out. EIG of cov matrix is one method, SVD of centered data matrix is another method, and then EIG and SVD can be performed by various methods as well. Does not matter - all of that is PCA. $\endgroup$
    – amoeba
    Oct 18, 2016 at 12:10
  • 1
    $\begingroup$ @amoeba Thanks for the clarification/correction on the terminology. What you say makes more sense to me, given that SVD and EIG are algebraic theorems/methods with a wider scope than PCA $\endgroup$
    – drake
    Oct 18, 2016 at 15:56
6
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.