forked from 170010011/fr
617 lines
23 KiB
Python
617 lines
23 KiB
Python
|
""" Principal Component Analysis.
|
||
|
"""
|
||
|
|
||
|
# Author: Alexandre Gramfort <alexandre.gramfort@inria.fr>
|
||
|
# Olivier Grisel <olivier.grisel@ensta.org>
|
||
|
# Mathieu Blondel <mathieu@mblondel.org>
|
||
|
# Denis A. Engemann <denis-alexander.engemann@inria.fr>
|
||
|
# Michael Eickenberg <michael.eickenberg@inria.fr>
|
||
|
# Giorgio Patrini <giorgio.patrini@anu.edu.au>
|
||
|
#
|
||
|
# License: BSD 3 clause
|
||
|
|
||
|
from math import log, sqrt
|
||
|
import numbers
|
||
|
|
||
|
import numpy as np
|
||
|
from scipy import linalg
|
||
|
from scipy.special import gammaln
|
||
|
from scipy.sparse import issparse
|
||
|
from scipy.sparse.linalg import svds
|
||
|
|
||
|
from ._base import _BasePCA
|
||
|
from ..utils import check_random_state
|
||
|
from ..utils._arpack import _init_arpack_v0
|
||
|
from ..utils.extmath import fast_logdet, randomized_svd, svd_flip
|
||
|
from ..utils.extmath import stable_cumsum
|
||
|
from ..utils.validation import check_is_fitted
|
||
|
from ..utils.validation import _deprecate_positional_args
|
||
|
|
||
|
|
||
|
def _assess_dimension(spectrum, rank, n_samples):
|
||
|
"""Compute the log-likelihood of a rank ``rank`` dataset.
|
||
|
|
||
|
The dataset is assumed to be embedded in gaussian noise of shape(n,
|
||
|
dimf) having spectrum ``spectrum``.
|
||
|
|
||
|
Parameters
|
||
|
----------
|
||
|
spectrum : ndarray of shape (n_features,)
|
||
|
Data spectrum.
|
||
|
rank : int
|
||
|
Tested rank value. It should be strictly lower than n_features,
|
||
|
otherwise the method isn't specified (division by zero in equation
|
||
|
(31) from the paper).
|
||
|
n_samples : int
|
||
|
Number of samples.
|
||
|
|
||
|
Returns
|
||
|
-------
|
||
|
ll : float
|
||
|
The log-likelihood.
|
||
|
|
||
|
Notes
|
||
|
-----
|
||
|
This implements the method of `Thomas P. Minka:
|
||
|
Automatic Choice of Dimensionality for PCA. NIPS 2000: 598-604`
|
||
|
"""
|
||
|
|
||
|
n_features = spectrum.shape[0]
|
||
|
if not 1 <= rank < n_features:
|
||
|
raise ValueError("the tested rank should be in [1, n_features - 1]")
|
||
|
|
||
|
eps = 1e-15
|
||
|
|
||
|
if spectrum[rank - 1] < eps:
|
||
|
# When the tested rank is associated with a small eigenvalue, there's
|
||
|
# no point in computing the log-likelihood: it's going to be very
|
||
|
# small and won't be the max anyway. Also, it can lead to numerical
|
||
|
# issues below when computing pa, in particular in log((spectrum[i] -
|
||
|
# spectrum[j]) because this will take the log of something very small.
|
||
|
return -np.inf
|
||
|
|
||
|
pu = -rank * log(2.)
|
||
|
for i in range(1, rank + 1):
|
||
|
pu += (gammaln((n_features - i + 1) / 2.) -
|
||
|
log(np.pi) * (n_features - i + 1) / 2.)
|
||
|
|
||
|
pl = np.sum(np.log(spectrum[:rank]))
|
||
|
pl = -pl * n_samples / 2.
|
||
|
|
||
|
v = max(eps, np.sum(spectrum[rank:]) / (n_features - rank))
|
||
|
pv = -np.log(v) * n_samples * (n_features - rank) / 2.
|
||
|
|
||
|
m = n_features * rank - rank * (rank + 1.) / 2.
|
||
|
pp = log(2. * np.pi) * (m + rank) / 2.
|
||
|
|
||
|
pa = 0.
|
||
|
spectrum_ = spectrum.copy()
|
||
|
spectrum_[rank:n_features] = v
|
||
|
for i in range(rank):
|
||
|
for j in range(i + 1, len(spectrum)):
|
||
|
pa += log((spectrum[i] - spectrum[j]) *
|
||
|
(1. / spectrum_[j] - 1. / spectrum_[i])) + log(n_samples)
|
||
|
|
||
|
ll = pu + pl + pv + pp - pa / 2. - rank * log(n_samples) / 2.
|
||
|
|
||
|
return ll
|
||
|
|
||
|
|
||
|
def _infer_dimension(spectrum, n_samples):
|
||
|
"""Infers the dimension of a dataset with a given spectrum.
|
||
|
|
||
|
The returned value will be in [1, n_features - 1].
|
||
|
"""
|
||
|
ll = np.empty_like(spectrum)
|
||
|
ll[0] = -np.inf # we don't want to return n_components = 0
|
||
|
for rank in range(1, spectrum.shape[0]):
|
||
|
ll[rank] = _assess_dimension(spectrum, rank, n_samples)
|
||
|
return ll.argmax()
|
||
|
|
||
|
|
||
|
class PCA(_BasePCA):
|
||
|
"""Principal component analysis (PCA).
|
||
|
|
||
|
Linear dimensionality reduction using Singular Value Decomposition of the
|
||
|
data to project it to a lower dimensional space. The input data is centered
|
||
|
but not scaled for each feature before applying the SVD.
|
||
|
|
||
|
It uses the LAPACK implementation of the full SVD or a randomized truncated
|
||
|
SVD by the method of Halko et al. 2009, depending on the shape of the input
|
||
|
data and the number of components to extract.
|
||
|
|
||
|
It can also use the scipy.sparse.linalg ARPACK implementation of the
|
||
|
truncated SVD.
|
||
|
|
||
|
Notice that this class does not support sparse input. See
|
||
|
:class:`TruncatedSVD` for an alternative with sparse data.
|
||
|
|
||
|
Read more in the :ref:`User Guide <PCA>`.
|
||
|
|
||
|
Parameters
|
||
|
----------
|
||
|
n_components : int, float or 'mle', default=None
|
||
|
Number of components to keep.
|
||
|
if n_components is not set all components are kept::
|
||
|
|
||
|
n_components == min(n_samples, n_features)
|
||
|
|
||
|
If ``n_components == 'mle'`` and ``svd_solver == 'full'``, Minka's
|
||
|
MLE is used to guess the dimension. Use of ``n_components == 'mle'``
|
||
|
will interpret ``svd_solver == 'auto'`` as ``svd_solver == 'full'``.
|
||
|
|
||
|
If ``0 < n_components < 1`` and ``svd_solver == 'full'``, select the
|
||
|
number of components such that the amount of variance that needs to be
|
||
|
explained is greater than the percentage specified by n_components.
|
||
|
|
||
|
If ``svd_solver == 'arpack'``, the number of components must be
|
||
|
strictly less than the minimum of n_features and n_samples.
|
||
|
|
||
|
Hence, the None case results in::
|
||
|
|
||
|
n_components == min(n_samples, n_features) - 1
|
||
|
|
||
|
copy : bool, default=True
|
||
|
If False, data passed to fit are overwritten and running
|
||
|
fit(X).transform(X) will not yield the expected results,
|
||
|
use fit_transform(X) instead.
|
||
|
|
||
|
whiten : bool, default=False
|
||
|
When True (False by default) the `components_` vectors are multiplied
|
||
|
by the square root of n_samples and then divided by the singular values
|
||
|
to ensure uncorrelated outputs with unit component-wise variances.
|
||
|
|
||
|
Whitening will remove some information from the transformed signal
|
||
|
(the relative variance scales of the components) but can sometime
|
||
|
improve the predictive accuracy of the downstream estimators by
|
||
|
making their data respect some hard-wired assumptions.
|
||
|
|
||
|
svd_solver : {'auto', 'full', 'arpack', 'randomized'}, default='auto'
|
||
|
If auto :
|
||
|
The solver is selected by a default policy based on `X.shape` and
|
||
|
`n_components`: if the input data is larger than 500x500 and the
|
||
|
number of components to extract is lower than 80% of the smallest
|
||
|
dimension of the data, then the more efficient 'randomized'
|
||
|
method is enabled. Otherwise the exact full SVD is computed and
|
||
|
optionally truncated afterwards.
|
||
|
If full :
|
||
|
run exact full SVD calling the standard LAPACK solver via
|
||
|
`scipy.linalg.svd` and select the components by postprocessing
|
||
|
If arpack :
|
||
|
run SVD truncated to n_components calling ARPACK solver via
|
||
|
`scipy.sparse.linalg.svds`. It requires strictly
|
||
|
0 < n_components < min(X.shape)
|
||
|
If randomized :
|
||
|
run randomized SVD by the method of Halko et al.
|
||
|
|
||
|
.. versionadded:: 0.18.0
|
||
|
|
||
|
tol : float, default=0.0
|
||
|
Tolerance for singular values computed by svd_solver == 'arpack'.
|
||
|
Must be of range [0.0, infinity).
|
||
|
|
||
|
.. versionadded:: 0.18.0
|
||
|
|
||
|
iterated_power : int or 'auto', default='auto'
|
||
|
Number of iterations for the power method computed by
|
||
|
svd_solver == 'randomized'.
|
||
|
Must be of range [0, infinity).
|
||
|
|
||
|
.. versionadded:: 0.18.0
|
||
|
|
||
|
random_state : int, RandomState instance or None, default=None
|
||
|
Used when the 'arpack' or 'randomized' solvers are used. Pass an int
|
||
|
for reproducible results across multiple function calls.
|
||
|
See :term:`Glossary <random_state>`.
|
||
|
|
||
|
.. versionadded:: 0.18.0
|
||
|
|
||
|
Attributes
|
||
|
----------
|
||
|
components_ : ndarray of shape (n_components, n_features)
|
||
|
Principal axes in feature space, representing the directions of
|
||
|
maximum variance in the data. The components are sorted by
|
||
|
``explained_variance_``.
|
||
|
|
||
|
explained_variance_ : ndarray of shape (n_components,)
|
||
|
The amount of variance explained by each of the selected components.
|
||
|
|
||
|
Equal to n_components largest eigenvalues
|
||
|
of the covariance matrix of X.
|
||
|
|
||
|
.. versionadded:: 0.18
|
||
|
|
||
|
explained_variance_ratio_ : ndarray of shape (n_components,)
|
||
|
Percentage of variance explained by each of the selected components.
|
||
|
|
||
|
If ``n_components`` is not set then all components are stored and the
|
||
|
sum of the ratios is equal to 1.0.
|
||
|
|
||
|
singular_values_ : ndarray of shape (n_components,)
|
||
|
The singular values corresponding to each of the selected components.
|
||
|
The singular values are equal to the 2-norms of the ``n_components``
|
||
|
variables in the lower-dimensional space.
|
||
|
|
||
|
.. versionadded:: 0.19
|
||
|
|
||
|
mean_ : ndarray of shape (n_features,)
|
||
|
Per-feature empirical mean, estimated from the training set.
|
||
|
|
||
|
Equal to `X.mean(axis=0)`.
|
||
|
|
||
|
n_components_ : int
|
||
|
The estimated number of components. When n_components is set
|
||
|
to 'mle' or a number between 0 and 1 (with svd_solver == 'full') this
|
||
|
number is estimated from input data. Otherwise it equals the parameter
|
||
|
n_components, or the lesser value of n_features and n_samples
|
||
|
if n_components is None.
|
||
|
|
||
|
n_features_ : int
|
||
|
Number of features in the training data.
|
||
|
|
||
|
n_samples_ : int
|
||
|
Number of samples in the training data.
|
||
|
|
||
|
noise_variance_ : float
|
||
|
The estimated noise covariance following the Probabilistic PCA model
|
||
|
from Tipping and Bishop 1999. See "Pattern Recognition and
|
||
|
Machine Learning" by C. Bishop, 12.2.1 p. 574 or
|
||
|
http://www.miketipping.com/papers/met-mppca.pdf. It is required to
|
||
|
compute the estimated data covariance and score samples.
|
||
|
|
||
|
Equal to the average of (min(n_features, n_samples) - n_components)
|
||
|
smallest eigenvalues of the covariance matrix of X.
|
||
|
|
||
|
See Also
|
||
|
--------
|
||
|
KernelPCA : Kernel Principal Component Analysis.
|
||
|
SparsePCA : Sparse Principal Component Analysis.
|
||
|
TruncatedSVD : Dimensionality reduction using truncated SVD.
|
||
|
IncrementalPCA : Incremental Principal Component Analysis.
|
||
|
|
||
|
References
|
||
|
----------
|
||
|
For n_components == 'mle', this class uses the method of *Minka, T. P.
|
||
|
"Automatic choice of dimensionality for PCA". In NIPS, pp. 598-604*
|
||
|
|
||
|
Implements the probabilistic PCA model from:
|
||
|
Tipping, M. E., and Bishop, C. M. (1999). "Probabilistic principal
|
||
|
component analysis". Journal of the Royal Statistical Society:
|
||
|
Series B (Statistical Methodology), 61(3), 611-622.
|
||
|
via the score and score_samples methods.
|
||
|
See http://www.miketipping.com/papers/met-mppca.pdf
|
||
|
|
||
|
For svd_solver == 'arpack', refer to `scipy.sparse.linalg.svds`.
|
||
|
|
||
|
For svd_solver == 'randomized', see:
|
||
|
*Halko, N., Martinsson, P. G., and Tropp, J. A. (2011).
|
||
|
"Finding structure with randomness: Probabilistic algorithms for
|
||
|
constructing approximate matrix decompositions".
|
||
|
SIAM review, 53(2), 217-288.* and also
|
||
|
*Martinsson, P. G., Rokhlin, V., and Tygert, M. (2011).
|
||
|
"A randomized algorithm for the decomposition of matrices".
|
||
|
Applied and Computational Harmonic Analysis, 30(1), 47-68.*
|
||
|
|
||
|
Examples
|
||
|
--------
|
||
|
>>> import numpy as np
|
||
|
>>> from sklearn.decomposition import PCA
|
||
|
>>> X = np.array([[-1, -1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])
|
||
|
>>> pca = PCA(n_components=2)
|
||
|
>>> pca.fit(X)
|
||
|
PCA(n_components=2)
|
||
|
>>> print(pca.explained_variance_ratio_)
|
||
|
[0.9924... 0.0075...]
|
||
|
>>> print(pca.singular_values_)
|
||
|
[6.30061... 0.54980...]
|
||
|
|
||
|
>>> pca = PCA(n_components=2, svd_solver='full')
|
||
|
>>> pca.fit(X)
|
||
|
PCA(n_components=2, svd_solver='full')
|
||
|
>>> print(pca.explained_variance_ratio_)
|
||
|
[0.9924... 0.00755...]
|
||
|
>>> print(pca.singular_values_)
|
||
|
[6.30061... 0.54980...]
|
||
|
|
||
|
>>> pca = PCA(n_components=1, svd_solver='arpack')
|
||
|
>>> pca.fit(X)
|
||
|
PCA(n_components=1, svd_solver='arpack')
|
||
|
>>> print(pca.explained_variance_ratio_)
|
||
|
[0.99244...]
|
||
|
>>> print(pca.singular_values_)
|
||
|
[6.30061...]
|
||
|
"""
|
||
|
@_deprecate_positional_args
|
||
|
def __init__(self, n_components=None, *, copy=True, whiten=False,
|
||
|
svd_solver='auto', tol=0.0, iterated_power='auto',
|
||
|
random_state=None):
|
||
|
self.n_components = n_components
|
||
|
self.copy = copy
|
||
|
self.whiten = whiten
|
||
|
self.svd_solver = svd_solver
|
||
|
self.tol = tol
|
||
|
self.iterated_power = iterated_power
|
||
|
self.random_state = random_state
|
||
|
|
||
|
def fit(self, X, y=None):
|
||
|
"""Fit the model with X.
|
||
|
|
||
|
Parameters
|
||
|
----------
|
||
|
X : array-like of shape (n_samples, n_features)
|
||
|
Training data, where n_samples is the number of samples
|
||
|
and n_features is the number of features.
|
||
|
|
||
|
y : Ignored
|
||
|
|
||
|
Returns
|
||
|
-------
|
||
|
self : object
|
||
|
Returns the instance itself.
|
||
|
"""
|
||
|
self._fit(X)
|
||
|
return self
|
||
|
|
||
|
def fit_transform(self, X, y=None):
|
||
|
"""Fit the model with X and apply the dimensionality reduction on X.
|
||
|
|
||
|
Parameters
|
||
|
----------
|
||
|
X : array-like of shape (n_samples, n_features)
|
||
|
Training data, where n_samples is the number of samples
|
||
|
and n_features is the number of features.
|
||
|
|
||
|
y : Ignored
|
||
|
|
||
|
Returns
|
||
|
-------
|
||
|
X_new : ndarray of shape (n_samples, n_components)
|
||
|
Transformed values.
|
||
|
|
||
|
Notes
|
||
|
-----
|
||
|
This method returns a Fortran-ordered array. To convert it to a
|
||
|
C-ordered array, use 'np.ascontiguousarray'.
|
||
|
"""
|
||
|
U, S, Vt = self._fit(X)
|
||
|
U = U[:, :self.n_components_]
|
||
|
|
||
|
if self.whiten:
|
||
|
# X_new = X * V / S * sqrt(n_samples) = U * sqrt(n_samples)
|
||
|
U *= sqrt(X.shape[0] - 1)
|
||
|
else:
|
||
|
# X_new = X * V = U * S * Vt * V = U * S
|
||
|
U *= S[:self.n_components_]
|
||
|
|
||
|
return U
|
||
|
|
||
|
def _fit(self, X):
|
||
|
"""Dispatch to the right submethod depending on the chosen solver."""
|
||
|
|
||
|
# Raise an error for sparse input.
|
||
|
# This is more informative than the generic one raised by check_array.
|
||
|
if issparse(X):
|
||
|
raise TypeError('PCA does not support sparse input. See '
|
||
|
'TruncatedSVD for a possible alternative.')
|
||
|
|
||
|
X = self._validate_data(X, dtype=[np.float64, np.float32],
|
||
|
ensure_2d=True, copy=self.copy)
|
||
|
|
||
|
# Handle n_components==None
|
||
|
if self.n_components is None:
|
||
|
if self.svd_solver != 'arpack':
|
||
|
n_components = min(X.shape)
|
||
|
else:
|
||
|
n_components = min(X.shape) - 1
|
||
|
else:
|
||
|
n_components = self.n_components
|
||
|
|
||
|
# Handle svd_solver
|
||
|
self._fit_svd_solver = self.svd_solver
|
||
|
if self._fit_svd_solver == 'auto':
|
||
|
# Small problem or n_components == 'mle', just call full PCA
|
||
|
if max(X.shape) <= 500 or n_components == 'mle':
|
||
|
self._fit_svd_solver = 'full'
|
||
|
elif n_components >= 1 and n_components < .8 * min(X.shape):
|
||
|
self._fit_svd_solver = 'randomized'
|
||
|
# This is also the case of n_components in (0,1)
|
||
|
else:
|
||
|
self._fit_svd_solver = 'full'
|
||
|
|
||
|
# Call different fits for either full or truncated SVD
|
||
|
if self._fit_svd_solver == 'full':
|
||
|
return self._fit_full(X, n_components)
|
||
|
elif self._fit_svd_solver in ['arpack', 'randomized']:
|
||
|
return self._fit_truncated(X, n_components, self._fit_svd_solver)
|
||
|
else:
|
||
|
raise ValueError("Unrecognized svd_solver='{0}'"
|
||
|
"".format(self._fit_svd_solver))
|
||
|
|
||
|
def _fit_full(self, X, n_components):
|
||
|
"""Fit the model by computing full SVD on X."""
|
||
|
n_samples, n_features = X.shape
|
||
|
|
||
|
if n_components == 'mle':
|
||
|
if n_samples < n_features:
|
||
|
raise ValueError("n_components='mle' is only supported "
|
||
|
"if n_samples >= n_features")
|
||
|
elif not 0 <= n_components <= min(n_samples, n_features):
|
||
|
raise ValueError("n_components=%r must be between 0 and "
|
||
|
"min(n_samples, n_features)=%r with "
|
||
|
"svd_solver='full'"
|
||
|
% (n_components, min(n_samples, n_features)))
|
||
|
elif n_components >= 1:
|
||
|
if not isinstance(n_components, numbers.Integral):
|
||
|
raise ValueError("n_components=%r must be of type int "
|
||
|
"when greater than or equal to 1, "
|
||
|
"was of type=%r"
|
||
|
% (n_components, type(n_components)))
|
||
|
|
||
|
# Center data
|
||
|
self.mean_ = np.mean(X, axis=0)
|
||
|
X -= self.mean_
|
||
|
|
||
|
U, S, Vt = linalg.svd(X, full_matrices=False)
|
||
|
# flip eigenvectors' sign to enforce deterministic output
|
||
|
U, Vt = svd_flip(U, Vt)
|
||
|
|
||
|
components_ = Vt
|
||
|
|
||
|
# 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.
|
||
|
|
||
|
# Postprocess the number of components required
|
||
|
if n_components == 'mle':
|
||
|
n_components = \
|
||
|
_infer_dimension(explained_variance_, n_samples)
|
||
|
elif 0 < n_components < 1.0:
|
||
|
# number of components for which the cumulated explained
|
||
|
# variance percentage is superior to the desired threshold
|
||
|
# side='right' ensures that number of features selected
|
||
|
# their variance is always greater than n_components float
|
||
|
# passed. More discussion in issue: #15669
|
||
|
ratio_cumsum = stable_cumsum(explained_variance_ratio_)
|
||
|
n_components = np.searchsorted(ratio_cumsum, n_components,
|
||
|
side='right') + 1
|
||
|
# Compute noise covariance using Probabilistic PCA model
|
||
|
# The sigma2 maximum likelihood (cf. eq. 12.46)
|
||
|
if n_components < min(n_features, n_samples):
|
||
|
self.noise_variance_ = explained_variance_[n_components:].mean()
|
||
|
else:
|
||
|
self.noise_variance_ = 0.
|
||
|
|
||
|
self.n_samples_, self.n_features_ = n_samples, n_features
|
||
|
self.components_ = components_[:n_components]
|
||
|
self.n_components_ = 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, Vt
|
||
|
|
||
|
def _fit_truncated(self, X, n_components, svd_solver):
|
||
|
"""Fit the model by computing truncated SVD (by ARPACK or randomized)
|
||
|
on X.
|
||
|
"""
|
||
|
n_samples, n_features = X.shape
|
||
|
|
||
|
if isinstance(n_components, str):
|
||
|
raise ValueError("n_components=%r cannot be a string "
|
||
|
"with svd_solver='%s'"
|
||
|
% (n_components, svd_solver))
|
||
|
elif not 1 <= n_components <= min(n_samples, n_features):
|
||
|
raise ValueError("n_components=%r must be between 1 and "
|
||
|
"min(n_samples, n_features)=%r with "
|
||
|
"svd_solver='%s'"
|
||
|
% (n_components, min(n_samples, n_features),
|
||
|
svd_solver))
|
||
|
elif not isinstance(n_components, numbers.Integral):
|
||
|
raise ValueError("n_components=%r must be of type int "
|
||
|
"when greater than or equal to 1, was of type=%r"
|
||
|
% (n_components, type(n_components)))
|
||
|
elif svd_solver == 'arpack' and n_components == min(n_samples,
|
||
|
n_features):
|
||
|
raise ValueError("n_components=%r must be strictly less than "
|
||
|
"min(n_samples, n_features)=%r with "
|
||
|
"svd_solver='%s'"
|
||
|
% (n_components, min(n_samples, n_features),
|
||
|
svd_solver))
|
||
|
|
||
|
random_state = check_random_state(self.random_state)
|
||
|
|
||
|
# Center data
|
||
|
self.mean_ = np.mean(X, axis=0)
|
||
|
X -= self.mean_
|
||
|
|
||
|
if svd_solver == 'arpack':
|
||
|
v0 = _init_arpack_v0(min(X.shape), random_state)
|
||
|
U, S, Vt = svds(X, k=n_components, tol=self.tol, v0=v0)
|
||
|
# svds doesn't abide by scipy.linalg.svd/randomized_svd
|
||
|
# conventions, so reverse its outputs.
|
||
|
S = S[::-1]
|
||
|
# flip eigenvectors' sign to enforce deterministic output
|
||
|
U, Vt = svd_flip(U[:, ::-1], Vt[::-1])
|
||
|
|
||
|
elif svd_solver == 'randomized':
|
||
|
# sign flipping is done inside
|
||
|
U, S, Vt = randomized_svd(X, n_components=n_components,
|
||
|
n_iter=self.iterated_power,
|
||
|
flip_sign=True,
|
||
|
random_state=random_state)
|
||
|
|
||
|
self.n_samples_, self.n_features_ = n_samples, n_features
|
||
|
self.components_ = Vt
|
||
|
self.n_components_ = n_components
|
||
|
|
||
|
# Get variance explained by singular values
|
||
|
self.explained_variance_ = (S ** 2) / (n_samples - 1)
|
||
|
total_var = np.var(X, ddof=1, 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_ < min(n_features, n_samples):
|
||
|
self.noise_variance_ = (total_var.sum() -
|
||
|
self.explained_variance_.sum())
|
||
|
self.noise_variance_ /= min(n_features, n_samples) - n_components
|
||
|
else:
|
||
|
self.noise_variance_ = 0.
|
||
|
|
||
|
return U, S, Vt
|
||
|
|
||
|
def score_samples(self, X):
|
||
|
"""Return the log-likelihood of each sample.
|
||
|
|
||
|
See. "Pattern Recognition and Machine Learning"
|
||
|
by C. Bishop, 12.2.1 p. 574
|
||
|
or http://www.miketipping.com/papers/met-mppca.pdf
|
||
|
|
||
|
Parameters
|
||
|
----------
|
||
|
X : array-like of shape (n_samples, n_features)
|
||
|
The data.
|
||
|
|
||
|
Returns
|
||
|
-------
|
||
|
ll : ndarray of shape (n_samples,)
|
||
|
Log-likelihood of each sample under the current model.
|
||
|
"""
|
||
|
check_is_fitted(self)
|
||
|
|
||
|
X = self._validate_data(X, dtype=[np.float64, np.float32], reset=False)
|
||
|
Xr = X - self.mean_
|
||
|
n_features = X.shape[1]
|
||
|
precision = self.get_precision()
|
||
|
log_like = -.5 * (Xr * (np.dot(Xr, precision))).sum(axis=1)
|
||
|
log_like -= .5 * (n_features * log(2. * np.pi) -
|
||
|
fast_logdet(precision))
|
||
|
return log_like
|
||
|
|
||
|
def score(self, X, y=None):
|
||
|
"""Return the average log-likelihood of all samples.
|
||
|
|
||
|
See. "Pattern Recognition and Machine Learning"
|
||
|
by C. Bishop, 12.2.1 p. 574
|
||
|
or http://www.miketipping.com/papers/met-mppca.pdf
|
||
|
|
||
|
Parameters
|
||
|
----------
|
||
|
X : array-like of shape (n_samples, n_features)
|
||
|
The data.
|
||
|
|
||
|
y : Ignored
|
||
|
|
||
|
Returns
|
||
|
-------
|
||
|
ll : float
|
||
|
Average log-likelihood of the samples under the current model.
|
||
|
"""
|
||
|
return np.mean(self.score_samples(X))
|
||
|
|
||
|
def _more_tags(self):
|
||
|
return {'preserves_dtype': [np.float64, np.float32]}
|