from warnings import warn import numpy as np import scipy.sparse import scipy.sparse.csgraph from sklearn.manifold import SpectralEmbedding from sklearn.metrics import pairwise_distances from sklearn.metrics.pairwise import _VALID_METRICS as SKLEARN_PAIRWISE_VALID_METRICS from umap.distances import pairwise_special_metric, SPECIAL_METRICS from umap.sparse import SPARSE_SPECIAL_METRICS, sparse_named_distances def component_layout( data, n_components, component_labels, dim, random_state, metric="euclidean", metric_kwds={}, ): """Provide a layout relating the separate connected components. This is done by taking the centroid of each component and then performing a spectral embedding of the centroids. Parameters ---------- data: array of shape (n_samples, n_features) The source data -- required so we can generate centroids for each connected component of the graph. n_components: int The number of distinct components to be layed out. component_labels: array of shape (n_samples) For each vertex in the graph the label of the component to which the vertex belongs. dim: int The chosen embedding dimension. metric: string or callable (optional, default 'euclidean') The metric used to measure distances among the source data points. metric_kwds: dict (optional, default {}) Keyword arguments to be passed to the metric function. If metric is 'precomputed', 'linkage' keyword can be used to specify 'average', 'complete', or 'single' linkage. Default is 'average' Returns ------- component_embedding: array of shape (n_components, dim) The ``dim``-dimensional embedding of the ``n_components``-many connected components. """ if data is None: # We don't have data to work with; just guess return np.random.random(size=(n_components, dim)) * 10.0 component_centroids = np.empty((n_components, data.shape[1]), dtype=np.float64) if metric == "precomputed": # cannot compute centroids from precomputed distances # instead, compute centroid distances using linkage distance_matrix = np.zeros((n_components, n_components), dtype=np.float64) linkage = metric_kwds.get("linkage", "average") if linkage == "average": linkage = np.mean elif linkage == "complete": linkage = np.max elif linkage == "single": linkage = np.min else: raise ValueError( "Unrecognized linkage '%s'. Please choose from " "'average', 'complete', or 'single'" % linkage ) for c_i in range(n_components): dm_i = data[component_labels == c_i] for c_j in range(c_i + 1, n_components): dist = linkage(dm_i[:, component_labels == c_j]) distance_matrix[c_i, c_j] = dist distance_matrix[c_j, c_i] = dist else: for label in range(n_components): component_centroids[label] = data[component_labels == label].mean(axis=0) if scipy.sparse.isspmatrix(component_centroids): warn( "Forcing component centroids to dense; if you are running out of " "memory then consider increasing n_neighbors." ) component_centroids = component_centroids.toarray() if metric in SPECIAL_METRICS: distance_matrix = pairwise_special_metric( component_centroids, metric=metric, kwds=metric_kwds, ) elif metric in SPARSE_SPECIAL_METRICS: distance_matrix = pairwise_special_metric( component_centroids, metric=SPARSE_SPECIAL_METRICS[metric], kwds=metric_kwds, ) else: if callable(metric) and scipy.sparse.isspmatrix(data): function_to_name_mapping = { sparse_named_distances[k]: k for k in set(SKLEARN_PAIRWISE_VALID_METRICS) & set(sparse_named_distances.keys()) } try: metric_name = function_to_name_mapping[metric] except KeyError: raise NotImplementedError( "Multicomponent layout for custom " "sparse metrics is not implemented at " "this time." ) distance_matrix = pairwise_distances( component_centroids, metric=metric_name, **metric_kwds ) else: distance_matrix = pairwise_distances( component_centroids, metric=metric, **metric_kwds ) affinity_matrix = np.exp(-(distance_matrix ** 2)) component_embedding = SpectralEmbedding( n_components=dim, affinity="precomputed", random_state=random_state ).fit_transform(affinity_matrix) component_embedding /= component_embedding.max() return component_embedding def multi_component_layout( data, graph, n_components, component_labels, dim, random_state, metric="euclidean", metric_kwds={}, ): """Specialised layout algorithm for dealing with graphs with many connected components. This will first fid relative positions for the components by spectrally embedding their centroids, then spectrally embed each individual connected component positioning them according to the centroid embeddings. This provides a decent embedding of each component while placing the components in good relative positions to one another. Parameters ---------- data: array of shape (n_samples, n_features) The source data -- required so we can generate centroids for each connected component of the graph. graph: sparse matrix The adjacency matrix of the graph to be emebdded. n_components: int The number of distinct components to be layed out. component_labels: array of shape (n_samples) For each vertex in the graph the label of the component to which the vertex belongs. dim: int The chosen embedding dimension. metric: string or callable (optional, default 'euclidean') The metric used to measure distances among the source data points. metric_kwds: dict (optional, default {}) Keyword arguments to be passed to the metric function. Returns ------- embedding: array of shape (n_samples, dim) The initial embedding of ``graph``. """ result = np.empty((graph.shape[0], dim), dtype=np.float32) if n_components > 2 * dim: meta_embedding = component_layout( data, n_components, component_labels, dim, random_state, metric=metric, metric_kwds=metric_kwds, ) else: k = int(np.ceil(n_components / 2.0)) base = np.hstack([np.eye(k), np.zeros((k, dim - k))]) meta_embedding = np.vstack([base, -base])[:n_components] for label in range(n_components): component_graph = graph.tocsr()[component_labels == label, :].tocsc() component_graph = component_graph[:, component_labels == label].tocoo() distances = pairwise_distances([meta_embedding[label]], meta_embedding) data_range = distances[distances > 0.0].min() / 2.0 if component_graph.shape[0] < 2 * dim or component_graph.shape[0] <= dim + 1: result[component_labels == label] = ( random_state.uniform( low=-data_range, high=data_range, size=(component_graph.shape[0], dim), ) + meta_embedding[label] ) continue diag_data = np.asarray(component_graph.sum(axis=0)) # standard Laplacian # D = scipy.sparse.spdiags(diag_data, 0, graph.shape[0], graph.shape[0]) # L = D - graph # Normalized Laplacian I = scipy.sparse.identity(component_graph.shape[0], dtype=np.float64) D = scipy.sparse.spdiags( 1.0 / np.sqrt(diag_data), 0, component_graph.shape[0], component_graph.shape[0], ) L = I - D * component_graph * D k = dim + 1 num_lanczos_vectors = max(2 * k + 1, int(np.sqrt(component_graph.shape[0]))) try: eigenvalues, eigenvectors = scipy.sparse.linalg.eigsh( L, k, which="SM", ncv=num_lanczos_vectors, tol=1e-4, v0=np.ones(L.shape[0]), maxiter=graph.shape[0] * 5, ) order = np.argsort(eigenvalues)[1:k] component_embedding = eigenvectors[:, order] expansion = data_range / np.max(np.abs(component_embedding)) component_embedding *= expansion result[component_labels == label] = ( component_embedding + meta_embedding[label] ) except scipy.sparse.linalg.ArpackError: warn( "WARNING: spectral initialisation failed! The eigenvector solver\n" "failed. This is likely due to too small an eigengap. Consider\n" "adding some noise or jitter to your data.\n\n" "Falling back to random initialisation!" ) result[component_labels == label] = ( random_state.uniform( low=-data_range, high=data_range, size=(component_graph.shape[0], dim), ) + meta_embedding[label] ) return result def spectral_layout(data, graph, dim, random_state, metric="euclidean", metric_kwds={}): """Given a graph compute the spectral embedding of the graph. This is simply the eigenvectors of the laplacian of the graph. Here we use the normalized laplacian. Parameters ---------- data: array of shape (n_samples, n_features) The source data graph: sparse matrix The (weighted) adjacency matrix of the graph as a sparse matrix. dim: int The dimension of the space into which to embed. random_state: numpy RandomState or equivalent A state capable being used as a numpy random state. Returns ------- embedding: array of shape (n_vertices, dim) The spectral embedding of the graph. """ n_samples = graph.shape[0] n_components, labels = scipy.sparse.csgraph.connected_components(graph) if n_components > 1: return multi_component_layout( data, graph, n_components, labels, dim, random_state, metric=metric, metric_kwds=metric_kwds, ) diag_data = np.asarray(graph.sum(axis=0)) # standard Laplacian # D = scipy.sparse.spdiags(diag_data, 0, graph.shape[0], graph.shape[0]) # L = D - graph # Normalized Laplacian I = scipy.sparse.identity(graph.shape[0], dtype=np.float64) D = scipy.sparse.spdiags( 1.0 / np.sqrt(diag_data), 0, graph.shape[0], graph.shape[0] ) L = I - D * graph * D k = dim + 1 num_lanczos_vectors = max(2 * k + 1, int(np.sqrt(graph.shape[0]))) try: if L.shape[0] < 2000000: eigenvalues, eigenvectors = scipy.sparse.linalg.eigsh( L, k, which="SM", ncv=num_lanczos_vectors, tol=1e-4, v0=np.ones(L.shape[0]), maxiter=graph.shape[0] * 5, ) else: eigenvalues, eigenvectors = scipy.sparse.linalg.lobpcg( L, random_state.normal(size=(L.shape[0], k)), largest=False, tol=1e-8 ) order = np.argsort(eigenvalues)[1:k] return eigenvectors[:, order] except scipy.sparse.linalg.ArpackError: warn( "WARNING: spectral initialisation failed! The eigenvector solver\n" "failed. This is likely due to too small an eigengap. Consider\n" "adding some noise or jitter to your data.\n\n" "Falling back to random initialisation!" ) return random_state.uniform(low=-10.0, high=10.0, size=(graph.shape[0], dim))