Repository: scikit-learn/scikit-learn
Pull request: [MRG+1] Made PCA expose the singular values
<!--
Thanks for contributing a pull request! Please ensure you have taken a look at
the contribution guidelines: https://github.com/scikit-learn/scikit-learn/blob/master/CONTRIBUTING.md#Contributing-Pull-Requests
-->
#### Reference Issue

<!-- Example: Fixes #1234 -->

Fixes #6955.
#### What does this implement/fix? Explain your changes.

I've made it so that the singular values from the underlying SVD are stored in the PCA classes under the attribute "singular_values_".

This was a trivial fix, and all I did was to save the singular values in an instance variable. I also added the attribute to the list of attributes in the documentations and added doc tests for them. I also added unit tests (`test_singular_values`) that cover the estimators PCA, IncrementalPCA and TruncatedSVD to the corresponding tests in sklearn/decomposition/tests.
#### Any other comments?

It may appear like there are loads of commits in this fix, but the commits listed below are some old stuff I was working on years ago, that are no longer in my fork. This fix only add changes to the modules: "decomposition/incremental_pca.py", "decomposition/pca.py" and "decomposition/truncated_pca.py" and to their corresponding unit tests, and the changes are made on the latest commit to the main scikit-learn repository.

<!--
Please be aware that we are a loose team of volunteers so patience is
necessary; assistance handling other issues is very welcome. We value
all user contributions, no matter how minor they are. If we are slow to
review, either the pull request needs some benchmarking, tinkering,
convincing, etc. or more likely the reviewers are simply busy. In either
case, we ask for your understanding during the review process.
For more information, see our FAQ on this topic:
http://scikit-learn.org/dev/faq.html#why-is-my-pull-request-not-getting-any-attention.

Thanks for contributing!
-->


Problem statement:
Store singular values in PCA instances
Is there a particular reason why classes such as PCA, TruncatedPCA do not store the singular values as instance variables after they have been computed? It would seem useful to store them since they are already computed so that they can be later accessed if need be. One use case would be to scale the right singular vectors by their singular values. I am happy to make the changes and push if there are no objections.

Hints from issue discussion:
From a quick look at it, it looks like you can compute them from the `explained_variance_` attribute.

``` py
n_samples = X.shape[0]
self.explained_variance_ = (S ** 2) / n_samples
```

`S` being the diagonal matrix with the singular values.

True, but a few comments:
- It should be possible to access the singular values directly as a class attribute (especially since it is computed anyways) rather than indirectly
- By PCA instances i meant to include PCA, TruncatedPCA and others. While indirectly computing the singular values from PCA is straightforward using `explained_variance_`, it is not possible (to my reading) from TruncatedPCA whose calculation for the `explained_variance_` is more convoluted. 
- It seems like IncrementalPCA already exposes a class attribute `self.singular_values_`, so it would make sense to expose a similar attribute in other PCA related classes as well. 

I guess there is no harm in storing S.

Commit message: ENH: Added the singular values to PCA by a singular_values_ instance variable as per Issue #6955.

Patch to review:
diff --git a/sklearn/decomposition/incremental_pca.py b/sklearn/decomposition/incremental_pca.py
index d95c465d2fc4a..d5771a283f9e6 100644
--- a/sklearn/decomposition/incremental_pca.py
+++ b/sklearn/decomposition/incremental_pca.py
@@ -71,7 +71,12 @@ class IncrementalPCA(_BasePCA):
     explained_variance_ratio_ : array, shape (n_components,)
         Percentage of variance explained by each of the selected components.
         If all components are stored, the sum of explained variances is equal
-        to 1.0
+        to 1.0.
+
+    singular_values_ : array, [n_components]
+        The singular values corresponsing to each of the selected components.
+        The singular values corresponds to the 2-norms of the ``n_components``
+        variables in the lower-dimensional space.
 
     mean_ : array, shape (n_features,)
         Per-feature empirical mean, aggregate over calls to ``partial_fit``.
@@ -166,6 +171,7 @@ def fit(self, X, y=None):
         self.singular_values_ = None
         self.explained_variance_ = None
         self.explained_variance_ratio_ = None
+        self.singular_values_ = None
         self.noise_variance_ = None
 
         X = check_array(X, copy=self.copy, dtype=[np.float64, np.float32])
diff --git a/sklearn/decomposition/pca.py b/sklearn/decomposition/pca.py
index b9b171301faeb..7bee0e87aee67 100644
--- a/sklearn/decomposition/pca.py
+++ b/sklearn/decomposition/pca.py
@@ -202,6 +202,11 @@ class PCA(_BasePCA):
         If ``n_components`` is not set then all components are stored and the
         sum of explained variances is equal to 1.0.
 
+    singular_values_ : array, [n_components]
+        The singular values corresponsing to each of the selected components.
+        The singular values corresponds to the 2-norms of the ``n_components``
+        variables in the lower-dimensional space.
+
     mean_ : array, [n_features]
         Per-feature empirical mean, estimated from the training set.
 
@@ -250,22 +255,28 @@ class PCA(_BasePCA):
     >>> pca.fit(X)
     PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,
       svd_solver='auto', tol=0.0, whiten=False)
-    >>> print(pca.explained_variance_ratio_) # doctest: +ELLIPSIS
+    >>> print(pca.explained_variance_ratio_)  # doctest: +ELLIPSIS
     [ 0.99244...  0.00755...]
+    >>> print(pca.singular_values_)  # doctest: +ELLIPSIS
+    [ 6.30061...  0.54980...]
 
     >>> pca = PCA(n_components=2, svd_solver='full')
     >>> pca.fit(X)                 # doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
     PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,
       svd_solver='full', tol=0.0, whiten=False)
-    >>> print(pca.explained_variance_ratio_) # doctest: +ELLIPSIS
+    >>> print(pca.explained_variance_ratio_)  # doctest: +ELLIPSIS
     [ 0.99244...  0.00755...]
+    >>> print(pca.singular_values_)  # doctest: +ELLIPSIS
+    [ 6.30061...  0.54980...]
 
     >>> pca = PCA(n_components=1, svd_solver='arpack')
     >>> pca.fit(X)
     PCA(copy=True, iterated_power='auto', n_components=1, random_state=None,
       svd_solver='arpack', tol=0.0, whiten=False)
-    >>> print(pca.explained_variance_ratio_) # doctest: +ELLIPSIS
+    >>> print(pca.explained_variance_ratio_)  # doctest: +ELLIPSIS
     [ 0.99244...]
+    >>> print(pca.singular_values_)  # doctest: +ELLIPSIS
+    [ 6.30061...]
 
     See also
     --------
@@ -385,6 +396,7 @@ def _fit_full(self, X, n_components):
         explained_variance_ = (S ** 2) / n_samples
         total_var = explained_variance_.sum()
         explained_variance_ratio_ = explained_variance_ / total_var
+        singular_values_ = S.copy()  # Store the singular values.
 
         # Postprocess the number of components required
         if n_components == 'mle':
@@ -409,6 +421,7 @@ def _fit_full(self, X, n_components):
         self.explained_variance_ = explained_variance_[:n_components]
         self.explained_variance_ratio_ = \
             explained_variance_ratio_[:n_components]
+        self.singular_values_ = singular_values_[:n_components]
 
         return U, S, V
 
@@ -463,6 +476,7 @@ def _fit_truncated(self, X, n_components, svd_solver):
         total_var = np.var(X, axis=0)
         self.explained_variance_ratio_ = \
             self.explained_variance_ / total_var.sum()
+        self.singular_values_ = S.copy()  # Store the singular values.
         if self.n_components_ < n_features:
             self.noise_variance_ = (total_var.sum() -
                                     self.explained_variance_.sum())
@@ -520,9 +534,11 @@ def score(self, X, y=None):
         return np.mean(self.score_samples(X))
 
 
-@deprecated("RandomizedPCA was deprecated in 0.18 and will be removed in 0.20. "
+@deprecated("RandomizedPCA was deprecated in 0.18 and will be removed in "
+            "0.20. "
             "Use PCA(svd_solver='randomized') instead. The new implementation "
-            "DOES NOT store whiten ``components_``. Apply transform to get them.")
+            "DOES NOT store whiten ``components_``. Apply transform to get "
+            "them.")
 class RandomizedPCA(BaseEstimator, TransformerMixin):
     """Principal component analysis (PCA) using randomized SVD
 
@@ -549,8 +565,8 @@ class RandomizedPCA(BaseEstimator, TransformerMixin):
         .. versionchanged:: 0.18
 
     whiten : bool, optional
-        When True (False by default) the `components_` vectors are multiplied by
-        the square root of (n_samples) and divided by the singular values to
+        When True (False by default) the `components_` vectors are multiplied
+        by the square root of (n_samples) and divided by the singular values to
         ensure uncorrelated outputs with unit component-wise variances.
 
         Whitening will remove some information from the transformed signal
@@ -569,8 +585,13 @@ class RandomizedPCA(BaseEstimator, TransformerMixin):
 
     explained_variance_ratio_ : array, [n_components]
         Percentage of variance explained by each of the selected components.
-        k is not set then all components are stored and the sum of explained
-        variances is equal to 1.0
+        If k is not set then all components are stored and the sum of explained
+        variances is equal to 1.0.
+
+    singular_values_ : array, [n_components]
+        The singular values corresponsing to each of the selected components.
+        The singular values corresponds to the 2-norms of the ``n_components``
+        variables in the lower-dimensional space.
 
     mean_ : array, [n_features]
         Per-feature empirical mean, estimated from the training set.
@@ -584,8 +605,10 @@ class RandomizedPCA(BaseEstimator, TransformerMixin):
     >>> pca.fit(X)                 # doctest: +ELLIPSIS +NORMALIZE_WHITESPACE
     RandomizedPCA(copy=True, iterated_power=2, n_components=2,
            random_state=None, whiten=False)
-    >>> print(pca.explained_variance_ratio_) # doctest: +ELLIPSIS
+    >>> print(pca.explained_variance_ratio_)  # doctest: +ELLIPSIS
     [ 0.99244...  0.00755...]
+    >>> print(pca.singular_values_)  # doctest: +ELLIPSIS
+    [ 6.30061...  0.54980...]
 
     See also
     --------
@@ -663,6 +686,7 @@ def _fit(self, X):
         self.explained_variance_ = exp_var = (S ** 2) / n_samples
         full_var = np.var(X, axis=0).sum()
         self.explained_variance_ratio_ = exp_var / full_var
+        self.singular_values_ = S.copy()  # Store the singular values.
 
         if self.whiten:
             self.components_ = V / S[:, np.newaxis] * sqrt(n_samples)
diff --git a/sklearn/decomposition/truncated_svd.py b/sklearn/decomposition/truncated_svd.py
index 3624a6153cbd4..c2c98b00a9664 100644
--- a/sklearn/decomposition/truncated_svd.py
+++ b/sklearn/decomposition/truncated_svd.py
@@ -71,26 +71,33 @@ class TruncatedSVD(BaseEstimator, TransformerMixin):
     ----------
     components_ : array, shape (n_components, n_features)
 
-    explained_variance_ratio_ : array, [n_components]
-        Percentage of variance explained by each of the selected components.
-
     explained_variance_ : array, [n_components]
         The variance of the training samples transformed by a projection to
         each component.
 
+    explained_variance_ratio_ : array, [n_components]
+        Percentage of variance explained by each of the selected components.
+
+    singular_values_ : array, [n_components]
+        The singular values corresponsing to each of the selected components.
+        The singular values corresponds to the 2-norms of the ``n_components``
+        variables in the lower-dimensional space.
+
     Examples
     --------
     >>> from sklearn.decomposition import TruncatedSVD
     >>> from sklearn.random_projection import sparse_random_matrix
     >>> X = sparse_random_matrix(100, 100, density=0.01, random_state=42)
     >>> svd = TruncatedSVD(n_components=5, n_iter=7, random_state=42)
-    >>> svd.fit(X) # doctest: +NORMALIZE_WHITESPACE
+    >>> svd.fit(X)  # doctest: +NORMALIZE_WHITESPACE
     TruncatedSVD(algorithm='randomized', n_components=5, n_iter=7,
             random_state=42, tol=0.0)
-    >>> print(svd.explained_variance_ratio_) # doctest: +ELLIPSIS
+    >>> print(svd.explained_variance_ratio_)  # doctest: +ELLIPSIS
     [ 0.0782... 0.0552... 0.0544... 0.0499... 0.0413...]
-    >>> print(svd.explained_variance_ratio_.sum()) # doctest: +ELLIPSIS
+    >>> print(svd.explained_variance_ratio_.sum())  # doctest: +ELLIPSIS
     0.279...
+    >>> print(svd.singular_values_)  # doctest: +ELLIPSIS
+    [ 2.6318... 2.2215... 2.1939... 2.1010... 1.9317...]
 
     See also
     --------
@@ -185,6 +192,8 @@ def fit_transform(self, X, y=None):
         else:
             full_var = np.var(X, axis=0).sum()
         self.explained_variance_ratio_ = exp_var / full_var
+        self.singular_values_ = Sigma  # Store the singular values.
+
         return X_transformed
 
     def transform(self, X):


Code review comments:
diff_hunk: @@ -663,6 +686,7 @@ def _fit(self, X):
         self.explained_variance_ = exp_var = (S ** 2) / n_samples
         full_var = np.var(X, axis=0).sum()
         self.explained_variance_ratio_ = exp_var / full_var
+        self.singular_values_ = S.copy()  # Store the singular values.
path: sklearn/decomposition/pca.py
text: @user1:
why do you need to copy?


@author:
No need to copy here. I've removed it.

diff_hunk: @@ -71,7 +71,12 @@ class IncrementalPCA(_BasePCA):
     explained_variance_ratio_ : array, shape (n_components,)
         Percentage of variance explained by each of the selected components.
         If all components are stored, the sum of explained variances is equal
-        to 1.0
+        to 1.0.
+
+    singular_values_ : array, [n_components]
path: sklearn/decomposition/incremental_pca.py
text: @user1:
please say `array, shape (n_components,)`


@author:
Done.

diff_hunk: @@ -71,26 +71,33 @@ class TruncatedSVD(BaseEstimator, TransformerMixin):
     ----------
     components_ : array, shape (n_components, n_features)
 
-    explained_variance_ratio_ : array, [n_components]
-        Percentage of variance explained by each of the selected components.
-
     explained_variance_ : array, [n_components]
         The variance of the training samples transformed by a projection to
         each component.
 
+    explained_variance_ratio_ : array, [n_components]
+        Percentage of variance explained by each of the selected components.
+
+    singular_values_ : array, [n_components]
+        The singular values corresponsing to each of the selected components.
+        The singular values corresponds to the 2-norms of the ``n_components``
path: sklearn/decomposition/truncated_svd.py
text: @user1:
Maybe say `are the 2-norms` or `are equal to`. You just used "corresponding". Also, there's a typo in "corresponding".


@author:
Done.