MNT use scipy.sparse.csgraph.shortest_path in Isomap (scikit-learn#20531)

Co-authored-by: Guillaume Lemaitre <[email protected]>

Author: TomDLT

doc/whats_new/v1.0.rst

@@ -43,6 +43,10 @@ random sampling procedures.
 - |Fix| :class:`manifold.TSNE` now avoids numerical underflow issues during
   affinity matrix computation.
 
+- |Fix| :class:`manifold.Isomap` now connects disconnected components of the
+  neighbors graph along some minimum distance pairs, instead of changing
+  every infinite distances to zero.
+
 - |Fix| The splitting criterion of :class:`tree.DecisionTreeClassifier` and
   :class:`tree.DecisionTreeRegressor` can be impacted by a fix in the handling
   of rounding errors. Previously some extra spurious splits could occur.
@@ -505,6 +509,12 @@ Changelog
   be scaled to have standard deviation 1e-4 in 1.2.
   :pr:`19491` by :user:`Dmitry Kobak <dkobak>`.
 
+- |Fix| :class:`manifold.Isomap` now uses `scipy.sparse.csgraph.shortest_path`
+  to compute the graph shortest path. It also connects disconnected components
+  of the neighbors graph along some minimum distance pairs, instead of changing
+  every infinite distances to zero. :pr:`20531` by `Roman Yurchak`_ and `Tom
+  Dupre la Tour`_.
+
 :mod:`sklearn.metrics`
 ......................
 
@@ -721,6 +731,10 @@ Changelog
   these functions were not documented and part from the public API.
   :pr:`20521` by :user:`Olivier Grisel <ogrisel>`.
 
+- |API| Fixed several bugs in :func:`utils.graph.graph_shortest_path`, which is
+  now deprecated. Use `scipy.sparse.csgraph.shortest_path` instead. :pr:`20531`
+  by `Tom Dupre la Tour`_.
+
 :mod:`sklearn.validation`
 .........................
 

sklearn/cluster/_agglomerative.py

@@ -15,12 +15,13 @@
 from scipy.sparse.csgraph import connected_components
 
 from ..base import BaseEstimator, ClusterMixin
-from ..metrics.pairwise import paired_distances, pairwise_distances
+from ..metrics.pairwise import paired_distances
 from ..neighbors import DistanceMetric
 from ..neighbors._dist_metrics import METRIC_MAPPING
 from ..utils import check_array
 from ..utils._fast_dict import IntFloatDict
 from ..utils.fixes import _astype_copy_false
+from ..utils.graph import _fix_connected_components
 from ..utils.validation import check_memory
 
 # mypy error: Module 'sklearn.cluster' has no attribute '_hierarchical_fast'
@@ -68,21 +69,14 @@ def _fix_connectivity(X, connectivity, affinity):
             stacklevel=2,
         )
         # XXX: Can we do without completing the matrix?
-        for i in range(n_connected_components):
-            idx_i = np.where(labels == i)[0]
-            Xi = X[idx_i]
-            for j in range(i):
-                idx_j = np.where(labels == j)[0]
-                Xj = X[idx_j]
-                if affinity == "precomputed":
-                    D = X[np.ix_(idx_i, idx_j)]
-                else:
-                    D = pairwise_distances(Xi, Xj, metric=affinity)
-                ii, jj = np.where(D == np.min(D))
-                ii = ii[0]
-                jj = jj[0]
-                connectivity[idx_i[ii], idx_j[jj]] = True
-                connectivity[idx_j[jj], idx_i[ii]] = True
+        connectivity = _fix_connected_components(
+            X=X,
+            graph=connectivity,
+            n_connected_components=n_connected_components,
+            component_labels=labels,
+            metric=affinity,
+            mode="connectivity",
+        )
 
     return connectivity, n_connected_components
 
@@ -661,7 +655,6 @@ def _single_linkage(*args, **kwargs):
     single=_single_linkage,
 )
 
-
 ###############################################################################
 # Functions for cutting hierarchical clustering tree
 

sklearn/manifold/_isomap.py

@@ -2,14 +2,20 @@
 
 # Author: Jake Vanderplas  -- <[email protected]>
 # License: BSD 3 clause (C) 2011
+import warnings
 
 import numpy as np
+import scipy
+from scipy.sparse.csgraph import shortest_path
+from scipy.sparse.csgraph import connected_components
+
 from ..base import BaseEstimator, TransformerMixin
 from ..neighbors import NearestNeighbors, kneighbors_graph
 from ..utils.validation import check_is_fitted
-from ..utils.graph import graph_shortest_path
 from ..decomposition import KernelPCA
 from ..preprocessing import KernelCenterer
+from ..utils.graph import _fix_connected_components
+from ..externals._packaging.version import parse as parse_version
 
 
 class Isomap(TransformerMixin, BaseEstimator):
@@ -186,9 +192,46 @@ def _fit_transform(self, X):
             n_jobs=self.n_jobs,
         )
 
-        self.dist_matrix_ = graph_shortest_path(
-            kng, method=self.path_method, directed=False
-        )
+        # Compute the number of connected components, and connect the different
+        # components to be able to compute a shortest path between all pairs
+        # of samples in the graph.
+        # Similar fix to cluster._agglomerative._fix_connectivity.
+        n_connected_components, labels = connected_components(kng)
+        if n_connected_components > 1:
+            if self.metric == "precomputed":
+                raise RuntimeError(
+                    "The number of connected components of the neighbors graph"
+                    f" is {n_connected_components} > 1. The graph cannot be "
+                    "completed with metric='precomputed', and Isomap cannot be"
+                    "fitted. Increase the number of neighbors to avoid this "
+                    "issue."
+                )
+            warnings.warn(
+                "The number of connected components of the neighbors graph "
+                f"is {n_connected_components} > 1. Completing the graph to fit"
+                " Isomap might be slow. Increase the number of neighbors to "
+                "avoid this issue.",
+                stacklevel=2,
+            )
+
+            # use array validated by NearestNeighbors
+            kng = _fix_connected_components(
+                X=self.nbrs_._fit_X,
+                graph=kng,
+                n_connected_components=n_connected_components,
+                component_labels=labels,
+                mode="distance",
+                metric=self.nbrs_.effective_metric_,
+                **self.nbrs_.effective_metric_params_,
+            )
+
+        if parse_version(scipy.__version__) < parse_version("1.3.2"):
+            # make identical samples have a nonzero distance, to account for
+            # issues in old scipy Floyd-Warshall implementation.
+            kng.data += 1e-15
+
+        self.dist_matrix_ = shortest_path(kng, method=self.path_method, directed=False)
+
         G = self.dist_matrix_ ** 2
         G *= -0.5
 

sklearn/manifold/tests/test_isomap.py

@@ -94,7 +94,7 @@ def test_transform():
     X, y = datasets.make_s_curve(n_samples, random_state=0)
 
     # Compute isomap embedding
-    iso = manifold.Isomap(n_components=n_components, n_neighbors=2)
+    iso = manifold.Isomap(n_components=n_components)
     X_iso = iso.fit_transform(X)
 
     # Re-embed a noisy version of the points
@@ -190,6 +190,25 @@ def test_sparse_input():
     for eigen_solver in eigen_solvers:
         for path_method in path_methods:
             clf = manifold.Isomap(
-                n_components=2, eigen_solver=eigen_solver, path_method=path_method
+                n_components=2,
+                eigen_solver=eigen_solver,
+                path_method=path_method,
+                n_neighbors=8,
             )
             clf.fit(X)
+
+
+def test_multiple_connected_components():
+    # Test that a warning is raised when the graph has multiple components
+    X = np.array([0, 1, 2, 5, 6, 7])[:, None]
+    with pytest.warns(UserWarning, match="number of connected components"):
+        manifold.Isomap(n_neighbors=2).fit(X)
+
+
+def test_multiple_connected_components_metric_precomputed():
+    # Test that an error is raised when the graph has multiple components
+    # and when the metric is "precomputed".
+    X = np.array([0, 1, 2, 5, 6, 7])[:, None]
+    X_graph = neighbors.kneighbors_graph(X, n_neighbors=2, mode="distance")
+    with pytest.raises(RuntimeError, match="number of connected components"):
+        manifold.Isomap(n_neighbors=1, metric="precomputed").fit(X_graph)

sklearn/utils/graph.py

@@ -10,9 +10,11 @@
 #          Jake Vanderplas <[email protected]>
 # License: BSD 3 clause
 
+import numpy as np
 from scipy import sparse
 
-from .graph_shortest_path import graph_shortest_path  # noqa
+from .deprecation import deprecated
+from ..metrics.pairwise import pairwise_distances
 
 
 ###############################################################################
@@ -67,3 +69,130 @@ def single_source_shortest_path_length(graph, source, *, cutoff=None):
             break
         level += 1
     return seen  # return all path lengths as dictionary
+
+
+@deprecated(
+    "`graph_shortest_path` is deprecated in 1.0 (renaming of 0.25) and will "
+    "be removed in 1.2. Use `scipy.sparse.csgraph.shortest_path` instead."
+)
+def graph_shortest_path(dist_matrix, directed=True, method="auto"):
+    """Shortest-path graph search on a positive directed or undirected graph.
+
+    Parameters
+    ----------
+    dist_matrix : arraylike or sparse matrix, shape = (N,N)
+        Array of positive distances.
+        If vertex i is connected to vertex j, then dist_matrix[i,j] gives
+        the distance between the vertices.
+        If vertex i is not connected to vertex j, then dist_matrix[i,j] = 0
+
+    directed : boolean
+        if True, then find the shortest path on a directed graph: only
+        progress from a point to its neighbors, not the other way around.
+        if False, then find the shortest path on an undirected graph: the
+        algorithm can progress from a point to its neighbors and vice versa.
+
+    method : string ['auto'|'FW'|'D']
+        method to use.  Options are
+        'auto' : attempt to choose the best method for the current problem
+        'FW' : Floyd-Warshall algorithm.  O[N^3]
+        'D' : Dijkstra's algorithm with Fibonacci stacks.  O[(k+log(N))N^2]
+
+    Returns
+    -------
+    G : np.ndarray, float, shape = [N,N]
+        G[i,j] gives the shortest distance from point i to point j
+        along the graph.
+
+    Notes
+    -----
+    As currently implemented, Dijkstra's algorithm does not work for
+    graphs with direction-dependent distances when directed == False.
+    i.e., if dist_matrix[i,j] and dist_matrix[j,i] are not equal and
+    both are nonzero, method='D' will not necessarily yield the correct
+    result.
+    Also, these routines have not been tested for graphs with negative
+    distances.  Negative distances can lead to infinite cycles that must
+    be handled by specialized algorithms.
+    """
+    return sparse.csgraph.shortest_path(dist_matrix, method=method, directed=directed)
+
+
+def _fix_connected_components(
+    X,
+    graph,
+    n_connected_components,
+    component_labels,
+    mode="distance",
+    metric="euclidean",
+    **kwargs,
+):
+    """Add connections to sparse graph to connect unconnected components.
+
+    For each pair of unconnected components, compute all pairwise distances
+    from one component to the other, and add a connection on the closest pair
+    of samples. This is a hacky way to get a graph with a single connected
+    component, which is necessary for example to compute a shortest path
+    between all pairs of samples in the graph.
+
+    Parameters
+    ----------
+    X : array of shape (n_samples, n_features) or (n_samples, n_samples)
+        Features to compute the pairwise distances. If `metric =
+        "precomputed"`, X is the matrix of pairwise distances.
+
+    graph : sparse matrix of shape (n_samples, n_samples)
+        Graph of connection between samples.
+
+    n_connected_components : int
+        Number of connected components, as computed by
+        `scipy.sparse.csgraph.connected_components`.
+
+    component_labels : array of shape (n_samples)
+        Labels of connected components, as computed by
+        `scipy.sparse.csgraph.connected_components`.
+
+    mode : {'connectivity', 'distance'}, default='distance'
+        Type of graph matrix: 'connectivity' corresponds to the connectivity
+        matrix with ones and zeros, and 'distance' corresponds to the distances
+        between neighbors according to the given metric.
+
+    metric : str
+        Metric used in `sklearn.metrics.pairwise.pairwise_distances`.
+
+    kwargs : kwargs
+        Keyword arguments passed to
+        `sklearn.metrics.pairwise.pairwise_distances`.
+
+    Returns
+    -------
+    graph : sparse matrix of shape (n_samples, n_samples)
+        Graph of connection between samples, with a single connected component.
+    """
+
+    for i in range(n_connected_components):
+        idx_i = np.flatnonzero(component_labels == i)
+        Xi = X[idx_i]
+        for j in range(i):
+            idx_j = np.flatnonzero(component_labels == j)
+            Xj = X[idx_j]
+
+            if metric == "precomputed":
+                D = X[np.ix_(idx_i, idx_j)]
+            else:
+                D = pairwise_distances(Xi, Xj, metric=metric, **kwargs)
+
+            ii, jj = np.unravel_index(D.argmin(axis=None), D.shape)
+            if mode == "connectivity":
+                graph[idx_i[ii], idx_j[jj]] = 1
+                graph[idx_j[jj], idx_i[ii]] = 1
+            elif mode == "distance":
+                graph[idx_i[ii], idx_j[jj]] = D[ii, jj]
+                graph[idx_j[jj], idx_i[ii]] = D[ii, jj]
+            else:
+                raise ValueError(
+                    "Unknown mode=%r, should be one of ['connectivity', 'distance']."
+                    % mode
+                )
+
+    return graph