""" Principal Component Analysis """ # Author: Alexandre Gramfort # Olivier Grisel # Mathieu Blondel # License: BSD Style. import numpy as np import warnings from scipy import linalg from math import log from ..base import BaseEstimator, TransformerMixin from ..utils import array2d, check_random_state, as_float_array from ..utils import atleast2d_or_csr from ..utils.extmath import fast_logdet from ..utils.extmath import safe_sparse_dot from ..utils.extmath import randomized_svd def _assess_dimension_(spectrum, rank, n_samples, n_features): """Compute the 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: array of shape (n) data spectrum rank: int, tested rank value n_samples: int, number of samples dim: int, embedding/empirical dimension 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` """ if rank > len(spectrum): raise ValueError("The tested rank cannot exceed the rank of the" " dataset") from scipy.special import gammaln pu = -rank * np.log(2) for i in range(rank): pu += (gammaln((n_features - i) / 2) - np.log(np.pi) * (n_features - i) / 2) pl = np.sum(np.log(spectrum[:rank])) pl = -pl * n_samples / 2 if rank == n_features: pv = 0 v = 1 else: v = 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 = np.log(2 * np.pi) * (m + rank + 1) / 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 += (np.log((spectrum[i] - spectrum[j]) * (1. / spectrum_[j] - 1. / spectrum_[i])) + np.log(n_samples)) ll = pu + pl + pv + pp - pa / 2 - rank * np.log(n_samples) / 2 return ll def _infer_dimension_(spectrum, n_samples, n_features): """Infers the dimension of a dataset of shape (n_samples, n_features) The dataset is described by its spectrum `spectrum`. """ ll = [] for rank in range(len(spectrum)): ll.append(_assess_dimension_(spectrum, rank, n_samples, n_features)) ll = np.array(ll) return ll.argmax() class PCA(BaseEstimator, TransformerMixin): """Principal component analysis (PCA) Linear dimensionality reduction using Singular Value Decomposition of the data and keeping only the most significant singular vectors to project the data to a lower dimensional space. This implementation uses the scipy.linalg implementation of the singular value decomposition. It only works for dense arrays and is not scalable to large dimensional data. The time complexity of this implementation is ``O(n ** 3)`` assuming n ~ n_samples ~ n_features. Parameters ---------- n_components : int, None or string 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', Minka\'s MLE is used to guess the dimension if ``0 < n_components < 1``, 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 copy : bool If False, data passed to fit are overwritten whiten : bool, optional When True (False by default) the `components_` vectors are divided by n_samples times 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 there data respect some hard-wired assumptions. Attributes ---------- `components_` : array, [n_components, n_features] Components with maximum variance. `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 Notes ----- For n_components='mle', this class uses the method of `Thomas P. Minka: Automatic Choice of Dimensionality for PCA. NIPS 2000: 598-604` Due to implementation subtleties of the Singular Value Decomposition (SVD), which is used in this implementation, running fit twice on the same matrix can lead to principal components with signs flipped (change in direction). For this reason, it is important to always use the same estimator object to transform data in a consistent fashion. 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(copy=True, n_components=2, whiten=False) >>> print(pca.explained_variance_ratio_) # doctest: +ELLIPSIS [ 0.99244... 0.00755...] See also -------- ProbabilisticPCA RandomizedPCA KernelPCA SparsePCA """ def __init__(self, n_components=None, copy=True, whiten=False): self.n_components = n_components self.copy = copy self.whiten = whiten def fit(self, X, y=None, **params): """Fit the model with X. Parameters ---------- X: array-like, shape (n_samples, n_features) Training data, where n_samples in the number of samples and n_features is the number of features. Returns ------- self : object Returns the instance itself. """ self._fit(X, **params) return self def fit_transform(self, X, y=None, **params): """Fit the model with X and apply the dimensionality reduction on X. Parameters ---------- X : array-like, shape (n_samples, n_features) Training data, where n_samples in the number of samples and n_features is the number of features. Returns ------- X_new : array-like, shape (n_samples, n_components) """ U, S, V = self._fit(X, **params) U = U[:, :self.n_components] if self.whiten: # X_new = X * V / S * sqrt(n_samples) = U * sqrt(n_samples) U *= np.sqrt(X.shape[0]) else: # X_new = X * V = U * S * V^T * V = U * S U *= S[:self.n_components] return U def _fit(self, X): X = array2d(X) n_samples, n_features = X.shape X = as_float_array(X, copy=self.copy) # Center data self.mean_ = np.mean(X, axis=0) X -= self.mean_ U, S, V = linalg.svd(X, full_matrices=False) self.explained_variance_ = (S ** 2) / n_samples self.explained_variance_ratio_ = (self.explained_variance_ / self.explained_variance_.sum()) if self.whiten: self.components_ = V / S[:, np.newaxis] * np.sqrt(n_samples) else: self.components_ = V if self.n_components == 'mle': if n_samples < n_features: raise ValueError("n_components='mle' is only supported " "if n_samples >= n_features") self.n_components = _infer_dimension_(self.explained_variance_, n_samples, n_features) elif (self.n_components is not None and 0 < self.n_components and self.n_components < 1.0): # number of components for which the cumulated explained variance # percentage is superior to the desired threshold ratio_cumsum = self.explained_variance_ratio_.cumsum() self.n_components = np.sum(ratio_cumsum < self.n_components) + 1 if self.n_components is not None: self.components_ = self.components_[:self.n_components, :] self.explained_variance_ = \ self.explained_variance_[:self.n_components] self.explained_variance_ratio_ = \ self.explained_variance_ratio_[:self.n_components] return (U, S, V) def transform(self, X): """Apply the dimensionality reduction on X. Parameters ---------- X : array-like, shape (n_samples, n_features) New data, where n_samples in the number of samples and n_features is the number of features. Returns ------- X_new : array-like, shape (n_samples, n_components) """ X = array2d(X) X_transformed = X - self.mean_ X_transformed = np.dot(X_transformed, self.components_.T) return X_transformed def inverse_transform(self, X): """Transform data back to its original space, i.e., return an input X_original whose transform would be X Parameters ---------- X : array-like, shape (n_samples, n_components) New data, where n_samples in the number of samples and n_components is the number of components. Returns ------- X_original array-like, shape (n_samples, n_features) Notes ----- If whitening is enabled, inverse_transform does not compute the exact inverse operation as transform. """ return np.dot(X, self.components_) + self.mean_ class ProbabilisticPCA(PCA): """Additional layer on top of PCA that adds a probabilistic evaluation""" __doc__ += PCA.__doc__ def fit(self, X, y=None, homoscedastic=True): """Additionally to PCA.fit, learns a covariance model Parameters ---------- X : array of shape(n_samples, n_features) The data to fit homoscedastic : bool, optional, If True, average variance across remaining dimensions """ PCA.fit(self, X) n_features = X.shape[1] self._dim = n_features Xr = X - self.mean_ Xr -= np.dot(np.dot(Xr, self.components_.T), self.components_) n_samples = X.shape[0] if n_features <= self.n_components: delta = np.zeros(n_features) elif homoscedastic: delta = ((Xr ** 2).sum() * np.ones(n_features) / (n_samples * n_features)) else: delta = (Xr ** 2).mean(0) / (n_features - self.n_components) self.covariance_ = np.diag(delta) n_components = self.n_components if n_components is None: n_components = n_features for k in range(n_components): add_cov = np.outer(self.components_[k], self.components_[k]) self.covariance_ += self.explained_variance_[k] * add_cov return self @property def dim(self): warnings.warn("Using dim is deprecated " "since version 0.12, and backward compatibility " "won't be maintained from version 0.14 onward. ", DeprecationWarning, stacklevel=2) return self._dim def score(self, X, y=None): """Return a score associated to new data Parameters ---------- X: array of shape(n_samples, n_features) The data to test Returns ------- ll: array of shape (n_samples), log-likelihood of each row of X under the current model """ Xr = X - self.mean_ n_features = X.shape[1] log_like = np.zeros(X.shape[0]) self.precision_ = linalg.inv(self.covariance_) log_like = -.5 * (Xr * (np.dot(Xr, self.precision_))).sum(axis=1) log_like -= .5 * (fast_logdet(self.covariance_) + n_features * log(2 * np.pi)) return log_like class RandomizedPCA(BaseEstimator, TransformerMixin): """Principal component analysis (PCA) using randomized SVD Linear dimensionality reduction using approximated Singular Value Decomposition of the data and keeping only the most significant singular vectors to project the data to a lower dimensional space. This implementation uses a randomized SVD implementation and can handle both scipy.sparse and numpy dense arrays as input. Parameters ---------- n_components : int Maximum number of components to keep: default is 50. copy : bool If False, data passed to fit are overwritten iterated_power : int, optional Number of iteration for the power method. 3 by default. whiten : bool, optional When True (False by default) the `components_` vectors are 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. random_state : int or RandomState instance or None (default) Pseudo Random Number generator seed control. If None, use the numpy.random singleton. Attributes ---------- `components_` : array, [n_components, n_features] Components with maximum variance. `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 Examples -------- >>> import numpy as np >>> from sklearn.decomposition import RandomizedPCA >>> X = np.array([[-1, -1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]]) >>> pca = RandomizedPCA(n_components=2) >>> pca.fit(X) # doctest: +ELLIPSIS +NORMALIZE_WHITESPACE RandomizedPCA(copy=True, iterated_power=3, n_components=2, random_state=None, whiten=False) >>> print(pca.explained_variance_ratio_) # doctest: +ELLIPSIS [ 0.99244... 0.00755...] See also -------- PCA ProbabilisticPCA References ---------- .. [Halko2009] `Finding structure with randomness: Stochastic algorithms for constructing approximate matrix decompositions Halko, et al., 2009 (arXiv:909)` .. [MRT] `A randomized algorithm for the decomposition of matrices Per-Gunnar Martinsson, Vladimir Rokhlin and Mark Tygert` """ def __init__(self, n_components=None, copy=True, iterated_power=3, whiten=False, random_state=None): self.n_components = n_components self.copy = copy self.iterated_power = iterated_power self.whiten = whiten self.mean_ = None self.random_state = random_state def fit(self, X, y=None): """Fit the model to the data X. Parameters ---------- X: array-like or scipy.sparse matrix, shape (n_samples, n_features) Training vector, where n_samples in the number of samples and n_features is the number of features. Returns ------- self : object Returns the instance itself. """ random_state = check_random_state(self.random_state) if not hasattr(X, 'todense'): # not a sparse matrix, ensure this is a 2D array X = np.atleast_2d(as_float_array(X, copy=self.copy)) n_samples = X.shape[0] if not hasattr(X, 'todense'): # Center data self.mean_ = np.mean(X, axis=0) X -= self.mean_ if self.n_components is None: n_components = X.shape[1] else: n_components = self.n_components U, S, V = randomized_svd(X, n_components, n_iter=self.iterated_power, random_state=random_state) self.explained_variance_ = exp_var = (S ** 2) / n_samples self.explained_variance_ratio_ = exp_var / exp_var.sum() if self.whiten: n = X.shape[0] self.components_ = V / S[:, np.newaxis] * np.sqrt(n) else: self.components_ = V return self def transform(self, X): """Apply dimensionality reduction on X. Parameters ---------- X : array-like or scipy.sparse matrix, shape (n_samples, n_features) New data, where n_samples in the number of samples and n_features is the number of features. Returns ------- X_new : array-like, shape (n_samples, n_components) """ X = atleast2d_or_csr(X) if self.mean_ is not None: X = X - self.mean_ X = safe_sparse_dot(X, self.components_.T) return X def inverse_transform(self, X): """Transform data back to its original space. Returns an array X_original whose transform would be X. Parameters ---------- X : array-like or scipy.sparse matrix, shape (n_samples, n_components) New data, where n_samples in the number of samples and n_components is the number of components. Returns ------- X_original array-like, shape (n_samples, n_features) Notes ----- If whitening is enabled, inverse_transform does not compute the exact inverse operation of transform. """ X_original = safe_sparse_dot(X, self.components_) if self.mean_ is not None: X_original = X_original + self.mean_ return X_original