"""Kernel Principal Components Analysis""" # Author: Mathieu Blondel # License: BSD Style. import numpy as np from scipy import linalg from ..utils.arpack import eigsh from ..base import BaseEstimator, TransformerMixin from ..preprocessing import KernelCenterer from ..metrics.pairwise import pairwise_kernels class KernelPCA(BaseEstimator, TransformerMixin): """Kernel Principal component analysis (KPCA) Non-linear dimensionality reduction through the use of kernels. Parameters ---------- n_components: int or None Number of components. If None, all non-zero components are kept. kernel: "linear" | "poly" | "rbf" | "sigmoid" | "cosine" | "precomputed" Kernel. Default: "linear" degree : int, optional Degree for poly, rbf and sigmoid kernels. Default: 3. gamma : float, optional Kernel coefficient for rbf and poly kernels. Default: 1/n_features. coef0 : float, optional Independent term in poly and sigmoid kernels. alpha: int Hyperparameter of the ridge regression that learns the inverse transform (when fit_inverse_transform=True). Default: 1.0 fit_inverse_transform: bool Learn the inverse transform for non-precomputed kernels. (i.e. learn to find the pre-image of a point) Default: False eigen_solver: string ['auto'|'dense'|'arpack'] Select eigensolver to use. If n_components is much less than the number of training samples, arpack may be more efficient than the dense eigensolver. tol: float convergence tolerance for arpack. Default: 0 (optimal value will be chosen by arpack) max_iter : int maximum number of iterations for arpack Default: None (optimal value will be chosen by arpack) Attributes ---------- `lambdas_`, `alphas_`: Eigenvalues and eigenvectors of the centered kernel matrix `dual_coef_`: Inverse transform matrix `X_transformed_fit_`: Projection of the fitted data on the kernel principal components References ---------- Kernel PCA was intoduced in: Bernhard Schoelkopf, Alexander J. Smola, and Klaus-Robert Mueller. 1999. Kernel principal component analysis. In Advances in kernel methods, MIT Press, Cambridge, MA, USA 327-352. """ def __init__(self, n_components=None, kernel="linear", gamma=None, degree=3, coef0=1, alpha=1.0, fit_inverse_transform=False, eigen_solver='auto', tol=0, max_iter=None): if fit_inverse_transform and kernel == 'precomputed': raise ValueError( "Cannot fit_inverse_transform with a precomputed kernel.") self.n_components = n_components self.kernel = kernel.lower() self.gamma = gamma self.degree = degree self.coef0 = coef0 self.alpha = alpha self.fit_inverse_transform = fit_inverse_transform self.eigen_solver = eigen_solver self.tol = tol self.max_iter = max_iter self._centerer = KernelCenterer() @property def _pairwise(self): return self.kernel == "precomputed" def _get_kernel(self, X, Y=None): params = {"gamma": self.gamma, "degree": self.degree, "coef0": self.coef0} try: return pairwise_kernels(X, Y, metric=self.kernel, filter_params=True, **params) except AttributeError: raise ValueError("%s is not a valid kernel. Valid kernels are: " "rbf, poly, sigmoid, linear and precomputed." % self.kernel) def _fit_transform(self, K): """ Fit's using kernel K""" # center kernel K = self._centerer.fit_transform(K) if self.n_components is None: n_components = K.shape[0] else: n_components = min(K.shape[0], self.n_components) # compute eigenvectors if self.eigen_solver == 'auto': if K.shape[0] > 200 and n_components < 10: eigen_solver = 'arpack' else: eigen_solver = 'dense' else: eigen_solver = self.eigen_solver if eigen_solver == 'dense': self.lambdas_, self.alphas_ = linalg.eigh( K, eigvals=(K.shape[0] - n_components, K.shape[0] - 1)) elif eigen_solver == 'arpack': self.lambdas_, self.alphas_ = eigsh(K, n_components, which="LA", tol=self.tol, maxiter=self.max_iter) # sort eignenvectors in descending order indices = self.lambdas_.argsort()[::-1] self.lambdas_ = self.lambdas_[indices] self.alphas_ = self.alphas_[:, indices] # remove eigenvectors with a zero eigenvalue self.alphas_ = self.alphas_[:, self.lambdas_ > 0] self.lambdas_ = self.lambdas_[self.lambdas_ > 0] return K def _fit_inverse_transform(self, X_transformed, X): if hasattr(X, "tocsr"): raise NotImplementedError("Inverse transform not implemented for " "sparse matrices!") n_samples = X_transformed.shape[0] K = self._get_kernel(X_transformed) K.flat[::n_samples + 1] += self.alpha self.dual_coef_ = linalg.solve(K, X, sym_pos=True, overwrite_a=True) self.X_transformed_fit_ = X_transformed def fit(self, X, y=None): """Fit the model from data in X. Parameters ---------- X: array-like, 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. """ K = self._get_kernel(X) self._fit_transform(K) if self.fit_inverse_transform: sqrt_lambdas = np.diag(np.sqrt(self.lambdas_)) X_transformed = np.dot(self.alphas_, sqrt_lambdas) self._fit_inverse_transform(X_transformed, X) self.X_fit_ = X return self def fit_transform(self, X, y=None, **params): """Fit the model from data in X and transform X. Parameters ---------- X: array-like, shape (n_samples, n_features) Training vector, 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) """ self.fit(X, **params) X_transformed = self.alphas_ * np.sqrt(self.lambdas_) if self.fit_inverse_transform: self._fit_inverse_transform(X_transformed, X) return X_transformed def transform(self, X): """Transform X. Parameters ---------- X: array-like, shape (n_samples, n_features) Returns ------- X_new: array-like, shape (n_samples, n_components) """ K = self._centerer.transform(self._get_kernel(X, self.X_fit_)) return np.dot(K, self.alphas_ / np.sqrt(self.lambdas_)) def inverse_transform(self, X): """Transform X back to original space. Parameters ---------- X: array-like, shape (n_samples, n_components) Returns ------- X_new: array-like, shape (n_samples, n_features) References ---------- "Learning to Find Pre-Images", G BakIr et al, 2004. """ if not self.fit_inverse_transform: raise ValueError("Inverse transform was not fitted!") K = self._get_kernel(X, self.X_transformed_fit_) return np.dot(K, self.dual_coef_)