"""Isomap for manifold learning""" # Author: Jake Vanderplas -- # License: BSD, (C) 2011 import numpy as np import warnings from ..base import BaseEstimator from ..neighbors import NearestNeighbors, kneighbors_graph from ..utils.graph import graph_shortest_path from ..decomposition import KernelPCA from ..preprocessing import KernelCenterer class Isomap(BaseEstimator): """Isomap Embedding Non-linear dimensionality reduction through Isometric Mapping Parameters ---------- n_neighbors : integer number of neighbors to consider for each point. n_components : integer number of coordinates for the manifold eigen_solver : ['auto'|'arpack'|'dense'] 'auto' : attempt to choose the most efficient solver for the given problem. 'arpack' : use Arnoldi decomposition to find the eigenvalues and eigenvectors. Note that arpack can handle both dense and sparse data efficiently 'dense' : use a direct solver (i.e. LAPACK) for the eigenvalue decomposition. tol : float convergence tolerance passed to arpack or lobpcg. not used if eigen_solver == 'dense' max_iter : integer maximum number of iterations for the arpack solver. not used if eigen_solver == 'dense' path_method : string ['auto'|'FW'|'D'] method to use in finding shortest path. 'auto' : attempt to choose the best algorithm automatically 'FW' : Floyd-Warshall algorithm 'D' : Dijkstra algorithm with Fibonacci Heaps neighbors_algorithm : string ['auto'|'brute'|'kd_tree'|'ball_tree'] algorithm to use for nearest neighbors search, passed to neighbors.NearestNeighbors instance Attributes ---------- `embedding_` : array-like, shape (n_samples, n_components) Stores the embedding vectors `kernel_pca_` : `KernelPCA` object used to implement the embedding `training_data_` : array-like, shape (n_samples, n_features) Stores the training data `nbrs_` : sklearn.neighbors.NearestNeighbors instance Stores nearest neighbors instance, including BallTree or KDtree if applicable. `dist_matrix_` : array-like, shape (n_samples, n_samples) Stores the geodesic distance matrix of training data References ---------- [1] Tenenbaum, J.B.; De Silva, V.; & Langford, J.C. A global geometric framework for nonlinear dimensionality reduction. Science 290 (5500) """ def __init__(self, n_neighbors=5, n_components=2, eigen_solver='auto', tol=0, max_iter=None, path_method='auto', neighbors_algorithm='auto', out_dim=None): if out_dim: warnings.warn("Parameter ``out_dim`` was renamed to " "``n_components`` and is now deprecated.", DeprecationWarning, stacklevel=2) self.out_dim = out_dim self.n_neighbors = n_neighbors self.n_components = n_components self.eigen_solver = eigen_solver self.tol = tol self.max_iter = max_iter self.path_method = path_method self.neighbors_algorithm = neighbors_algorithm self.nbrs_ = NearestNeighbors(n_neighbors=n_neighbors, algorithm=neighbors_algorithm) def _fit_transform(self, X): if self.out_dim: warnings.warn("Parameter ``out_dim`` was renamed to " "``n_components`` and is now deprecated.", DeprecationWarning, stacklevel=3) self.n_components = self.out_dim self.out_dim = None self.nbrs_.fit(X) self.training_data_ = self.nbrs_._fit_X self.kernel_pca_ = KernelPCA(n_components=self.n_components, kernel="precomputed", eigen_solver=self.eigen_solver, tol=self.tol, max_iter=self.max_iter) kng = kneighbors_graph(self.nbrs_, self.n_neighbors, mode='distance') self.dist_matrix_ = graph_shortest_path(kng, method=self.path_method, directed=False) G = self.dist_matrix_ ** 2 G *= -0.5 self.embedding_ = self.kernel_pca_.fit_transform(G) def reconstruction_error(self): """Compute the reconstruction error for the embedding. Returns ------- reconstruction_error : float Notes ------- The cost function of an isomap embedding is ``E = frobenius_norm[K(D) - K(D_fit)] / n_samples`` Where D is the matrix of distances for the input data X, D_fit is the matrix of distances for the output embedding X_fit, and K is the isomap kernel: ``K(D) = -0.5 * (I - 1/n_samples) * D^2 * (I - 1/n_samples)`` """ G = -0.5 * self.dist_matrix_ ** 2 G_center = KernelCenterer().fit_transform(G) evals = self.kernel_pca_.lambdas_ return np.sqrt(np.sum(G_center ** 2) - np.sum(evals ** 2)) / G.shape[0] def fit(self, X, y=None): """Compute the embedding vectors for data X Parameters ---------- X : {array-like, sparse matrix, BallTree, cKDTree, NearestNeighbors} Sample data, shape = (n_samples, n_features), in the form of a numpy array, sparse array, precomputed tree, or NearestNeighbors object. Returns ------- self : returns an instance of self. """ self._fit_transform(X) return self def fit_transform(self, X, y=None): """Fit the model from data in X and transform X. Parameters ---------- X: {array-like, sparse matrix, BallTree, cKDTree} 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_transform(X) return self.embedding_ def transform(self, X): """Transform X. This is implemented by linking the points X into the graph of geodesic distances of the training data. First the `n_neighbors` nearest neighbors of X are found in the training data, and from these the shortest geodesic distances from each point in X to each point in the training data are computed in order to construct the kernel. The embedding of X is the projection of this kernel onto the embedding vectors of the training set. Parameters ---------- X: array-like, shape (n_samples, n_features) Returns ------- X_new: array-like, shape (n_samples, n_components) """ distances, indices = self.nbrs_.kneighbors(X, return_distance=True) #Create the graph of shortest distances from X to self.training_data_ # via the nearest neighbors of X. #This can be done as a single array operation, but it potentially # takes a lot of memory. To avoid that, use a loop: G_X = np.zeros((X.shape[0], self.training_data_.shape[0])) for i in range(X.shape[0]): G_X[i] = np.min((self.dist_matrix_[indices[i]] + distances[i][:, None]), 0) G_X **= 2 G_X *= -0.5 return self.kernel_pca_.transform(G_X)