Repository: scikit-learn/scikit-learn
Pull request: ENH propagate `eigen_tol` to all eigen solver
#### Reference Issues/PRs
closes https://github.com/scikit-learn/scikit-learn/issues/21243
Resolves #11968 (stalled)

#### What does this implement/fix? Explain your changes.
PR #11968: Propagate the `eigen_tol` for the 'arpack', 'amg', and 'lobpcg' solvers.

This PR: 

- Adds tests

#### Any other comments?
Follow up PRs to: 
- Add `"auto"` option for `eigen_tol`
- Begins deprecation where appropriate to update default `eigen_tol=0` to `eigen_tol="auto"`


Problem statement:
eigen_tol in _spectral_embedding.py does not propagate to solvers other than arpack
### Describe the bug

In `_spectral_embedding.py` the option 
```
    eigen_tol : float, default=0.0
        Stopping criterion for eigendecomposition of the Laplacian matrix
        when ``eigen_solver='arpack'``.
```
does not appear to propagate to solvers 'amg' and 'lobpcg', e.g., adding `eigen_tol ` and changing its value in `examples/cluster/plot_coin_segmentation.py` affects the outcome only if `eigen_solver='arpack'`, which is the default but not for the other two: 'lobpcg' and 'amg'.

### Steps/Code to Reproduce

See above

### Expected Results

`eigen_tol ` works for all supported solvers

### Actual Results

`eigen_tol ` works for `eigen_solver='arpack' only and is ignored otherwise

### While at it

Also introduce and pass `max_iter` param to all eigen solvers, although maybe it should be named `eigen_max_iter` to avoid confusion.

### Versions

System:
    python: 3.7.3 (v3.7.3:ef4ec6ed12, Mar 25 2019, 22:22:05) [MSC v.1916 64 bit (AMD64)]
executable: C:\Users\Name\AppData\Local\Programs\Python\Python37\python.exe    
   machine: Windows-10-10.0.19041-SP0

Python dependencies:
          pip: 21.2.4
   setuptools: 58.1.0
      sklearn: 1.0
        numpy: 1.21.2
        scipy: 1.7.1
       Cython: 0.29.24
       pandas: 1.3.3
   matplotlib: 3.4.3
       joblib: 1.0.1
threadpoolctl: 2.2.0

Hints from issue discussion:
Currently, `eigen_tol` is different between solvers i.e. `logpcg` uses `1e-5` and arpack uses `0.0`. Should the solvers have different default tolerances?
1. The solvers should ideally have the same default accuracy, but `arpack` and `lobpcg` likely have different stopping criteria, so it may need some testing to figure out what the default numeric values of tolerances should be to make the accuracy the same. 
2. `arpack` formally uses `0.0`, but mathematically the tolerance cannot be zero - the iterations would never stop due to numeric errors. Thus,  `arpack` probably treats `0.0` as `None`, i.e., sets up some tolerance >0 internally.
3. At least `amg` is just `lobpcg` with some extra tricks, so tolerances in `amg` and `lobpcg` should be the same. 
4. Figuring out and setting default tolerances is a different issue, unrelated to the reported bug; e.g., see #21194
To be backward compatible, we can have `eigen_tol=None` and pass in 0 to arpack, which is machine precision internally. For `lobpcg` and `amg`, we continue to pass in `1e-5` when `eigen_tol=None`. If `eigen_tol` is a float, we can pass it in directly.

Does this sound reasonable?
> To be backward compatible, we can have `eigen_tol=None` and pass in 0 to arpack, which is machine precision internally. For `lobpcg` and `amg`, we continue to pass in `1e-5` when `eigen_tol=None`. If `eigen_tol` is a float, we can pass it in directly.
> 
> Does this sound reasonable?

@thomasjpfan  Yes. 

In a future PR, may be someone should also introduce parameter `maxit`, pass it to all solvers to overwrite the internal defaults, and set some defaults in `spectral_clustering`, just like for `eigen_tol`. That would be a desirable enhancement.

On a side note, I am unsure that 0 in `arpack` is indeed machine precision internally. But if it is 0, I am confident that it's not a good default for `arpack` here. E.g., something like `1e-8` seems to work just fine in #21148 `plot_coin_segmentation.py` while cutting 2x the compute time. But that's also not relevant to the present issue and should be addressed separately.
It would indeed be interesting to be able to pass a `max_iter` param to the eigen solvers, although maybe it should be named `eigen_max_iter` to avoid confusion with any outer-level `max_iter` loop in the estimator itself.

Having `eigen_max_iter` exposed would probably better help study the influence of the eigensolver convergence on the outer statistical estimation problem.

> On a side note, I am unsure that 0 in arpack is indeed machine precision internally. But if it is 0, I am confident that it's not a good default for arpack here. E.g., something like 1e-8 seems to work just fine in #21148 plot_coin_segmentation.py while cutting 2x the compute time. But that's also not relevant to the present issue and should be addressed separately.

That would deserve more investigation in tracking down the meaning of `tol=0` in ARPACK. Also to move forward I like the plan suggested by @thomasjpfan in https://github.com/scikit-learn/scikit-learn/issues/21243#issuecomment-939076860 to move backward compat problems out of the way and decouple the discussion of the choice of the default value for the arpack solver from the fix of the original problem reported here (namely `eigen_tol` propagation for other solvers).
In addition to `spectral_embedding`, the same changes need to be done in `SpectralEmbedding`, `spectral_clustering`, and `SpectralClustering`.

The PR already did those changes: https://github.com/scikit-learn/scikit-learn/pull/11968
It did most of the backward compatibility things that are discussed here. The tricky one is the `max_iter=2000` in the second call to `lobpcg`.

My concern is about the tests. I don't really know how to trigger some convergence issue to be sure that we probably handle the exception.
> I don't really know how to trigger some convergence issue

@glemaitre The theory tells us that convergence of any eigensolver will be slow for a random initiation if there is a small separation of wanted and unwanted eigenvalues, e.g., finding one, the smallest eigenvalue, in the case where the second smallest is very close. In spectral clustering that would correspond to determining two clusters in the case where there are three or more nearly identical clusters. E.g., make a Gaussian point cloud (a blob) in the positive quadrant in 2D far enough from the origin, then flip the signs of the 2 components of all the points to create symmetric copies of the blob in every quadrant. Now you have 4 identical clusters, due to symmetry, just perturb them slightly and use any distance to generate the graph adjacency matrix. As you try to partition this union of 4 nearly identical clusters into two clusters, where do you partition, vertically or horizontally?  This ambiguity should manifest itself in slow convergence in the solver.

I have not tried such an example, but the theory predicts that it should work. Could be actually a good new demo...
While at it, may be also fix https://github.com/scikit-learn/scikit-learn/issues/14713 ?

Commit message: Address review feedback

Patch to review:
diff --git a/doc/whats_new/v1.2.rst b/doc/whats_new/v1.2.rst
index 5851100a358c5..10f4451afcbfb 100644
--- a/doc/whats_new/v1.2.rst
+++ b/doc/whats_new/v1.2.rst
@@ -19,6 +19,12 @@ parameters, may produce different models from the previous version. This often
 occurs due to changes in the modelling logic (bug fixes or enhancements), or in
 random sampling procedures.
 
+- |Enhancement| The default `eigen_tol` for :class:`cluster.SpectralClustering`,
+  :class:`manifold.SpectralEmbedding`, :func:`cluster.spectral_clustering`,
+  and :func:`manifold.spectral_embedding` is now `None` when using the `'amg'`
+  or `'lobpcg'` solvers. This change improves numerical stability of the
+  solver, but may result in a different model.
+
 Changelog
 ---------
 
@@ -44,6 +50,13 @@ Changelog
 - |Enhancement| :class:`cluster.Birch` now preserves dtype for `numpy.float32`
   inputs. :pr:`22968` by `Meekail Zain <micky774>`.
 
+- |Enhancement| :class:`cluster.SpectralClustering` and
+  :func:`cluster.spectral_clustering` now propogates the `eigen_tol` parameter
+  to all choices of `eigen_solver`. Includes a new option `eigen_tol="auto"`
+  and begins deprecation to change the default from `eigen_tol=0` to
+  `eigen_tol="auto"` in version 1.3.
+  :pr:`23210` by :user:`Meekail Zain <micky774>`.
+
 - |Enhancement| :class:`cluster.KMeans` and :class:`cluster.MiniBatchKMeans`
   now accept a new `'auto'` option for `n_init` which changes the number of
   random initializations to one when using `init='k-means++'` for efficiency.
@@ -108,13 +121,17 @@ Changelog
   deterministic SVD used by the randomized SVD algorithm.
   :pr:`20617` by :user:`Srinath Kailasa <skailasa>`
 
-Code and Documentation Contributors
------------------------------------
 
-Thanks to everyone who has contributed to the maintenance and improvement of
-the project since version 1.1, including:
+:mod:`sklearn.manifold`
+.......................
 
-TODO: update at the time of the release.
+- |Enhancement| Adds `eigen_tol` parameter to
+  :class:`manifold.SpectralEmbedding`. Both :func:`manifold.spectral_embedding`
+  and :class:`manifold.SpectralEmbedding` now propogate `eigen_tol` to all
+  choices of `eigen_solver`. Includes a new option `eigen_tol="auto"`
+  and begins deprecation to change the default from `eigen_tol=0` to
+  `eigen_tol="auto"` in version 1.3.
+  :pr:`23210` by :user:`Meekail Zain <micky774>`.
 
 :mod:`sklearn.svm`
 ..................
@@ -122,3 +139,11 @@ TODO: update at the time of the release.
 - |API| The `class_weight_` attribute is now deprecated for
   :class:`svm.NuSVR`, :class:`svm.SVR`, :class:`svm.OneClassSVM`.
   :pr:`22898` by :user:`Meekail Zain <micky774>`.
+
+Code and Documentation Contributors
+-----------------------------------
+
+Thanks to everyone who has contributed to the maintenance and improvement of
+the project since version 1.1, including:
+
+TODO: update at the time of the release.
diff --git a/sklearn/cluster/_spectral.py b/sklearn/cluster/_spectral.py
index 390b567c0d0bb..967f37e5bb621 100644
--- a/sklearn/cluster/_spectral.py
+++ b/sklearn/cluster/_spectral.py
@@ -198,7 +198,7 @@ def spectral_clustering(
     eigen_solver=None,
     random_state=None,
     n_init=10,
-    eigen_tol=0.0,
+    eigen_tol="auto",
     assign_labels="kmeans",
     verbose=False,
 ):
@@ -259,9 +259,23 @@ def spectral_clustering(
         consecutive runs in terms of inertia. Only used if
         ``assign_labels='kmeans'``.
 
-    eigen_tol : float, default=0.0
-        Stopping criterion for eigendecomposition of the Laplacian matrix
-        when using arpack eigen_solver.
+    eigen_tol : float, default="auto"
+        Stopping criterion for eigendecomposition of the Laplacian matrix.
+        If `eigen_tol="auto"` then the passed tolerance will depend on the
+        `eigen_solver`:
+
+        - If `eigen_solver="arpack"`, then `eigen_tol=0.0`;
+        - If `eigen_solver="lobpcg"` or `eigen_solver="amg"`, then
+         `eigen_tol=None` which configures the underlying `lobpcg` solver to
+          automatically resolve the value according to their heuristics. See,
+          :func:`scipy.sparse.linalg.lobpcg` for details.
+
+        Note that when using `eigen_solver="lobpcg"` or `eigen_solver="amg"`
+        values of `tol<1e-5` may lead to convergence issues and should be
+        avoided.
+
+        .. versionadded:: 1.2
+           Added 'auto' option.
 
     assign_labels : {'kmeans', 'discretize', 'cluster_qr'}, default='kmeans'
         The strategy to use to assign labels in the embedding
@@ -461,9 +475,23 @@ class SpectralClustering(ClusterMixin, BaseEstimator):
         Number of neighbors to use when constructing the affinity matrix using
         the nearest neighbors method. Ignored for ``affinity='rbf'``.
 
-    eigen_tol : float, default=0.0
-        Stopping criterion for eigendecomposition of the Laplacian matrix
-        when ``eigen_solver='arpack'``.
+    eigen_tol : float, default="auto"
+        Stopping criterion for eigendecomposition of the Laplacian matrix.
+        If `eigen_tol="auto"` then the passed tolerance will depend on the
+        `eigen_solver`:
+
+        - If `eigen_solver="arpack"`, then `eigen_tol=0.0`;
+        - If `eigen_solver="lobpcg"` or `eigen_solver="amg"`, then
+         `eigen_tol=None` which configures the underlying `lobpcg` solver to
+          automatically resolve the value according to their heuristics. See,
+          :func:`scipy.sparse.linalg.lobpcg` for details.
+
+        Note that when using `eigen_solver="lobpcg"` or `eigen_solver="amg"`
+        values of `tol<1e-5` may lead to convergence issues and should be
+        avoided.
+
+        .. versionadded:: 1.2
+           Added 'auto' option.
 
     assign_labels : {'kmeans', 'discretize', 'cluster_qr'}, default='kmeans'
         The strategy for assigning labels in the embedding space. There are two
@@ -598,7 +626,7 @@ def __init__(
         gamma=1.0,
         affinity="rbf",
         n_neighbors=10,
-        eigen_tol=0.0,
+        eigen_tol="auto",
         assign_labels="kmeans",
         degree=3,
         coef0=1,
@@ -694,7 +722,7 @@ def fit(self, X, y=None):
             include_boundaries="left",
         )
 
-        if self.eigen_solver == "arpack":
+        if self.eigen_tol != "auto":
             check_scalar(
                 self.eigen_tol,
                 "eigen_tol",
diff --git a/sklearn/manifold/_spectral_embedding.py b/sklearn/manifold/_spectral_embedding.py
index 54a72313491ee..1170551e4faf8 100644
--- a/sklearn/manifold/_spectral_embedding.py
+++ b/sklearn/manifold/_spectral_embedding.py
@@ -145,7 +145,7 @@ def spectral_embedding(
     n_components=8,
     eigen_solver=None,
     random_state=None,
-    eigen_tol=0.0,
+    eigen_tol="auto",
     norm_laplacian=True,
     drop_first=True,
 ):
@@ -197,9 +197,22 @@ def spectral_embedding(
             https://github.com/pyamg/pyamg/issues/139 for further
             information.
 
-    eigen_tol : float, default=0.0
-        Stopping criterion for eigendecomposition of the Laplacian matrix
-        when using arpack eigen_solver.
+    eigen_tol : float, default="auto"
+        Stopping criterion for eigendecomposition of the Laplacian matrix.
+        If `eigen_tol="auto"` then the passed tolerance will depend on the
+        `eigen_solver`:
+
+        - If `eigen_solver="arpack"`, then `eigen_tol=0.0`;
+        - If `eigen_solver="lobpcg"` or `eigen_solver="amg"`, then
+         `eigen_tol=None` which configures the underlying `lobpcg` solver to
+          automatically resolve the value according to their heuristics. See,
+          :func:`scipy.sparse.linalg.lobpcg` for details.
+
+        Note that when using `eigen_solver="amg"` values of `tol<1e-5` may lead
+        to convergence issues and should be avoided.
+
+        .. versionadded:: 1.2
+           Added 'auto' option.
 
     norm_laplacian : bool, default=True
         If True, then compute symmetric normalized Laplacian.
@@ -293,10 +306,11 @@ def spectral_embedding(
         try:
             # We are computing the opposite of the laplacian inplace so as
             # to spare a memory allocation of a possibly very large array
+            tol = 0 if eigen_tol == "auto" else eigen_tol
             laplacian *= -1
             v0 = _init_arpack_v0(laplacian.shape[0], random_state)
             _, diffusion_map = eigsh(
-                laplacian, k=n_components, sigma=1.0, which="LM", tol=eigen_tol, v0=v0
+                laplacian, k=n_components, sigma=1.0, which="LM", tol=tol, v0=v0
             )
             embedding = diffusion_map.T[n_components::-1]
             if norm_laplacian:
@@ -338,7 +352,9 @@ def spectral_embedding(
         X = random_state.standard_normal(size=(laplacian.shape[0], n_components + 1))
         X[:, 0] = dd.ravel()
         X = X.astype(laplacian.dtype)
-        _, diffusion_map = lobpcg(laplacian, X, M=M, tol=1.0e-5, largest=False)
+
+        tol = None if eigen_tol == "auto" else eigen_tol
+        _, diffusion_map = lobpcg(laplacian, X, M=M, tol=tol, largest=False)
         embedding = diffusion_map.T
         if norm_laplacian:
             # recover u = D^-1/2 x from the eigenvector output x
@@ -371,8 +387,9 @@ def spectral_embedding(
             )
             X[:, 0] = dd.ravel()
             X = X.astype(laplacian.dtype)
+            tol = None if eigen_tol == "auto" else eigen_tol
             _, diffusion_map = lobpcg(
-                laplacian, X, tol=1e-5, largest=False, maxiter=2000
+                laplacian, X, tol=tol, largest=False, maxiter=2000
             )
             embedding = diffusion_map.T[:n_components]
             if norm_laplacian:
@@ -444,6 +461,23 @@ class SpectralEmbedding(BaseEstimator):
         to be installed. It can be faster on very large, sparse problems.
         If None, then ``'arpack'`` is used.
 
+    eigen_tol : float, default="auto"
+        Stopping criterion for eigendecomposition of the Laplacian matrix.
+        If `eigen_tol="auto"` then the passed tolerance will depend on the
+        `eigen_solver`:
+
+        - If `eigen_solver="arpack"`, then `eigen_tol=0.0`;
+        - If `eigen_solver="lobpcg"` or `eigen_solver="amg"`, then
+         `eigen_tol=None` which configures the underlying `lobpcg` solver to
+          automatically resolve the value according to their heuristics. See,
+          :func:`scipy.sparse.linalg.lobpcg` for details.
+
+        Note that when using `eigen_solver="lobpcg"` or `eigen_solver="amg"`
+        values of `tol<1e-5` may lead to convergence issues and should be
+        avoided.
+
+        .. versionadded:: 1.2
+
     n_neighbors : int, default=None
         Number of nearest neighbors for nearest_neighbors graph building.
         If None, n_neighbors will be set to max(n_samples/10, 1).
@@ -516,6 +550,7 @@ def __init__(
         gamma=None,
         random_state=None,
         eigen_solver=None,
+        eigen_tol="auto",
         n_neighbors=None,
         n_jobs=None,
     ):
@@ -524,6 +559,7 @@ def __init__(
         self.gamma = gamma
         self.random_state = random_state
         self.eigen_solver = eigen_solver
+        self.eigen_tol = eigen_tol
         self.n_neighbors = n_neighbors
         self.n_jobs = n_jobs
 
@@ -641,6 +677,7 @@ def fit(self, X, y=None):
             affinity_matrix,
             n_components=self.n_components,
             eigen_solver=self.eigen_solver,
+            eigen_tol=self.eigen_tol,
             random_state=random_state,
         )
         return self
diff --git a/sklearn/manifold/tests/test_spectral_embedding.py b/sklearn/manifold/tests/test_spectral_embedding.py
index 38e555643b034..f1d33d3a67ba7 100644
--- a/sklearn/manifold/tests/test_spectral_embedding.py
+++ b/sklearn/manifold/tests/test_spectral_embedding.py
@@ -1,3 +1,4 @@
+from unittest.mock import Mock
 import pytest
 
 import numpy as np
@@ -5,19 +6,21 @@
 from scipy import sparse
 from scipy.sparse import csgraph
 from scipy.linalg import eigh
+from scipy.sparse.linalg import eigsh
 
-from sklearn.manifold import SpectralEmbedding
+from sklearn.manifold import SpectralEmbedding, _spectral_embedding
 from sklearn.manifold._spectral_embedding import _graph_is_connected
 from sklearn.manifold._spectral_embedding import _graph_connected_component
 from sklearn.manifold import spectral_embedding
 from sklearn.metrics.pairwise import rbf_kernel
-from sklearn.metrics import normalized_mutual_info_score
+from sklearn.metrics import normalized_mutual_info_score, pairwise_distances
 from sklearn.neighbors import NearestNeighbors
 from sklearn.cluster import KMeans
 from sklearn.datasets import make_blobs
 from sklearn.utils.extmath import _deterministic_vector_sign_flip
 from sklearn.utils._testing import assert_array_almost_equal
 from sklearn.utils._testing import assert_array_equal
+from sklearn.utils.fixes import lobpcg
 
 try:
     from pyamg import smoothed_aggregation_solver  # noqa
@@ -480,3 +483,29 @@ def test_error_pyamg_not_available():
     err_msg = "The eigen_solver was set to 'amg', but pyamg is not available."
     with pytest.raises(ValueError, match=err_msg):
         se_precomp.fit_transform(S)
+
+
+@pytest.mark.parametrize("solver", ["arpack", "amg", "lobpcg"])
+def test_spectral_eigen_tol_auto(monkeypatch, solver):
+    """Test that `eigen_tol="auto"` is resolved correctly"""
+    if solver == "amg":
+        if not pyamg_available:
+            pytest.skip("PyAMG is not available.")
+    X = make_blobs(
+        n_samples=200, random_state=0, centers=[[1, 1], [-1, -1]], cluster_std=0.01
+    )[0]
+    D = pairwise_distances(X)  # Distance matrix
+    S = np.max(D) - D  # Similarity matrix
+
+    solver_func = eigsh if solver == "arpack" else lobpcg
+    default_value = 0 if solver == "arpack" else None
+    if solver == "amg":
+        S = sparse.csr_matrix(S)
+
+    mocked_solver = Mock(side_effect=solver_func)
+
+    monkeypatch.setattr(_spectral_embedding, solver_func.__qualname__, mocked_solver)
+
+    spectral_embedding(S, random_state=42, eigen_solver=solver, eigen_tol="auto")
+    mocked_solver.assert_called()
+    assert mocked_solver.call_args.kwargs["tol"] == default_value


Code review comments:
diff_hunk: @@ -480,3 +483,29 @@ def test_error_pyamg_not_available():
     err_msg = "The eigen_solver was set to 'amg', but pyamg is not available."
     with pytest.raises(ValueError, match=err_msg):
         se_precomp.fit_transform(S)
+
+
+@pytest.mark.parametrize("solver", ["arpack", "amg", "lobpcg"])
+def test_spectral_eigen_tol_auto(monkeypatch, solver):
+    """Test that `eigen_tol="auto"` is resolved correctly"""
+    if solver == "amg":
+        if not pyamg_available:
+            pytest.skip("PyAMG is not available.")
path: sklearn/manifold/tests/test_spectral_embedding.py
text: @user1:
```suggestion
    if solver == "amg" and not pyamg_available:
        pytest.skip("PyAMG is not available.")
```

diff_hunk: @@ -480,3 +483,29 @@ def test_error_pyamg_not_available():
     err_msg = "The eigen_solver was set to 'amg', but pyamg is not available."
     with pytest.raises(ValueError, match=err_msg):
         se_precomp.fit_transform(S)
+
+
+@pytest.mark.parametrize("solver", ["arpack", "amg", "lobpcg"])
+def test_spectral_eigen_tol_auto(monkeypatch, solver):
+    """Test that `eigen_tol="auto"` is resolved correctly"""
+    if solver == "amg":
+        if not pyamg_available:
+            pytest.skip("PyAMG is not available.")
+    X = make_blobs(
+        n_samples=200, random_state=0, centers=[[1, 1], [-1, -1]], cluster_std=0.01
+    )[0]
path: sklearn/manifold/tests/test_spectral_embedding.py
text: @user1:
Nit: The `[0]` is easier to miss:

```suggestion
    X, _ = make_blobs(
        n_samples=200, random_state=0, centers=[[1, 1], [-1, -1]], cluster_std=0.01
    )
```