Repository: scikit-learn/scikit-learn
Pull request: [MRG] Fix KMeans convergence when tol==0
Fixes #17428

When ``tol==0``, convergence based on the norm of the difference of the centers between 2 iterations can fail due to rounding errors making the norm not exactly 0 (and leave iterations go until max_iter is reached).

This PR implements convergence based on the difference of the labels between 2 iterations when tol == 0.

It does not impact performances. Here's a timing comparison between setting tol=0 and tol=1e-16 (both converge after 60 iterations to the same clustering):
```py
from sklearn.cluster import KMeans
from sklearn.datasets import make_blobs

X, _ = make_blobs(n_samples=100000, centers=100, n_features=10)

km = KMeans(n_clusters=100, tol=0, n_init=1, init='random', random_state=0)
%timeit km.fit(X)
2.13 s ± 3.17 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

km = KMeans(n_clusters=100, tol=1e-16, n_init=1, init='random', random_state=0)
%timeit km.fit(X)
2.14 s ± 16.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
```

Problem statement:
locally failing kmeans convergence test (WSL)
```
> pytest -sv sklearn/cluster/tests/test_k_means.py -k test_kmeans_convergence
```
fails for me for both algorithms:
```pythontb
>       assert km.n_iter_ < 300
E       AssertionError: assert 300 < 300
E        +  where 300 = KMeans(algorithm='full', n_clusters=5, n_init=1, random_state=0, tol=0).n_iter_
```
show_versions:
```
System:
    python: 3.8.3 (default, May 19 2020, 18:47:26)  [GCC 7.3.0]
executable: /home/andy/anaconda3/envs/sklearndev/bin/python
   machine: Linux-4.4.0-18362-Microsoft-x86_64-with-glibc2.10

Python dependencies:
          pip: 20.0.2
   setuptools: 46.4.0.post20200518
      sklearn: 0.24.dev0
        numpy: 1.18.1
        scipy: 1.4.1
       Cython: 0.29.17
       pandas: 1.0.3
   matplotlib: 3.1.3
       joblib: 0.15.1
threadpoolctl: 2.1.0

Built with OpenMP: True
```

Hints from issue discussion:
cc @jeremiedbb 
Inertia basically stays the same after iteration 50, but there seems to be some floating point issues maybe?

```
Iteration 48, inertia 40826.39061567664
Iteration 49, inertia 40826.37236625354
Iteration 50, inertia 40826.37236625354
Iteration 51, inertia 40826.37236625354
Iteration 52, inertia 40826.372366253534
Iteration 53, inertia 40826.372366253534
Iteration 54, inertia 40826.372366253534
Iteration 55, inertia 40826.372366253534
Iteration 56, inertia 40826.372366253534
Iteration 57, inertia 40826.372366253534
Iteration 58, inertia 40826.37236625354
```

showing center shift we can see:
```
Iteration 49, inertia 40826.37236625354, center shift 2.000920841596814e-31
Iteration 50, inertia 40826.37236625354, center shift 2.272958446241511e-31
```
Somewhat related to #17428.
I wonder if we should *also* check that no cluster centers are relabeled? This is tolerance is meant to be a weaker convergence criterion in my understanding. If nothing is relabeled, we certainly don't want to keep going.
cc @kno10
I see two solutions: remove the test because it's unstable, or add an extra condition into the convergence check for ``tol=0`` to also check whether the labels changed.
It does not come from relabelling indefinitely. When `tol=0` we are checking that `norm(old_centers - new_centers) == 0`. Floating point errors can make the norm to never be exactly 0 but ~ machine precision. 

I already added a comment in the docstring that using tol=0 is not advised.
No sure if you saw my comment just before yours:
https://github.com/scikit-learn/scikit-learn/issues/17428#issuecomment-637829543

What I mean is that either we should make sure the test that's actually performed when using tol=0 is more stable than a floating point comparison, or we shouldn't have that setting in the tests.
Initially this test was added in a PR which goal was to make sure that when tol=0, the iteration loop don't run until max_iter (due to strict inequality check before).

Since we can't guarantee that due to floating point errors, I think we should just remove the test.
I'm fine with that, but then maybe we should also raise a warning if someone provides ``tol=0``? Asking for the classical behavior of no changes in assignments seems like a reasonable usecase but it would require carrying around the old labels. which is not too bad.
I don't have a preference. I'm fine with both. I'd also be fine with silently switch to tol=eps when user provides tol=0 :)
The standard approach is simply to set a flag when any point is reassigned when computing the new assignment; and converge only when no point is relabeled.
I understand that providing a "tol" parameter is more in line with other machine learning algorithms that cannot converge otherwise. But since k-means *finishes* with "no more changes", this definitely should be an option to have. This is a *feature* of k-means; and quite some people will expect the algorithm to have reached such a stable state, where every point *is* really assigned to the closest center - if you have tol>0, that means there probably still is some point that is *not* assigned to the closest center. This may be confusing to users if the algorithm stops too early.
If a user wants floating-point tolerance convergence, he can pass tol=2e-308 or similar. tol=0 should really be a "no more reassignment happening" to be compatible with the common definition of the standard k-means algorithm as "stop if no point is reassigned".

Commit message: address comments

Patch to review:
diff --git a/sklearn/cluster/_kmeans.py b/sklearn/cluster/_kmeans.py
index a42a3b309a57d..0fb75f26d9ad1 100644
--- a/sklearn/cluster/_kmeans.py
+++ b/sklearn/cluster/_kmeans.py
@@ -226,8 +226,8 @@ def k_means(X, n_clusters, *, sample_weight=None, init='k-means++',
         Relative tolerance with regards to Frobenius norm of the difference
         in the cluster centers of two consecutive iterations to declare
         convergence.
-        It's not advised to set `tol=0` since convergence might never be
-        declared due to rounding errors. Use a very small number instead.
+        If set to 0, convergence will be declared when the labels are the
+        same between two consecutive iterations.
 
     random_state : int, RandomState instance, default=None
         Determines random number generation for centroid initialization. Use
@@ -358,6 +358,7 @@ def _kmeans_single_elkan(X, sample_weight, centers_init, max_iter=300,
     centers_new = np.zeros_like(centers)
     weight_in_clusters = np.zeros(n_clusters, dtype=X.dtype)
     labels = np.full(n_samples, -1, dtype=np.int32)
+    labels_old = labels.copy()
     center_half_distances = euclidean_distances(centers) / 2
     distance_next_center = np.partition(np.asarray(center_half_distances),
                                         kth=1, axis=0)[1]
@@ -377,6 +378,8 @@ def _kmeans_single_elkan(X, sample_weight, centers_init, max_iter=300,
     init_bounds(X, centers, center_half_distances,
                 labels, upper_bounds, lower_bounds)
 
+    strict_convergence = False
+
     for i in range(max_iter):
         elkan_iter(X, sample_weight, centers, centers_new,
                    weight_in_clusters, center_half_distances,
@@ -395,14 +398,25 @@ def _kmeans_single_elkan(X, sample_weight, centers_init, max_iter=300,
 
         centers, centers_new = centers_new, centers
 
-        center_shift_tot = (center_shift**2).sum()
-        if center_shift_tot <= tol:
+        if np.array_equal(labels, labels_old):
+            # First check the labels for strict convergence because
+            # center_shift might not be exactly 0 due to rounding errors.
             if verbose:
-                print(f"Converged at iteration {i}: center shift "
-                      f"{center_shift_tot} within tolerance {tol}.")
+                print(f"Converged at iteration {i}: strict convergence.")
+            strict_convergence = True
             break
+        else:
+            # No strict convergence, check for tol based convergence.
+            center_shift_tot = (center_shift**2).sum()
+            if center_shift_tot <= tol:
+                if verbose:
+                    print(f"Converged at iteration {i}: center shift "
+                          f"{center_shift_tot} within tolerance {tol}.")
+                break
 
-    if center_shift_tot > 0:
+        labels_old[:] = labels
+
+    if not strict_convergence:
         # rerun E-step so that predicted labels match cluster centers
         elkan_iter(X, sample_weight, centers, centers, weight_in_clusters,
                    center_half_distances, distance_next_center,
@@ -473,6 +487,7 @@ def _kmeans_single_lloyd(X, sample_weight, centers_init, max_iter=300,
     centers = centers_init
     centers_new = np.zeros_like(centers)
     labels = np.full(X.shape[0], -1, dtype=np.int32)
+    labels_old = labels.copy()
     weight_in_clusters = np.zeros(n_clusters, dtype=X.dtype)
     center_shift = np.zeros(n_clusters, dtype=X.dtype)
 
@@ -483,6 +498,8 @@ def _kmeans_single_lloyd(X, sample_weight, centers_init, max_iter=300,
         lloyd_iter = lloyd_iter_chunked_dense
         _inertia = _inertia_dense
 
+    strict_convergence = False
+
     # Threadpoolctl context to limit the number of threads in second level of
     # nested parallelism (i.e. BLAS) to avoid oversubsciption.
     with threadpool_limits(limits=1, user_api="blas"):
@@ -496,15 +513,25 @@ def _kmeans_single_lloyd(X, sample_weight, centers_init, max_iter=300,
 
             centers, centers_new = centers_new, centers
 
-            center_shift_tot = (center_shift**2).sum()
-            if center_shift_tot <= tol:
+            if np.array_equal(labels, labels_old):
+                # First check the labels for strict convergence because
+                # center_shift might not be exactly 0 due to rounding errors.
                 if verbose:
-                    print("Converged at iteration {0}: "
-                          "center shift {1} within tolerance {2}"
-                          .format(i, center_shift_tot, tol))
+                    print(f"Converged at iteration {i}: strict convergence.")
+                strict_convergence = True
                 break
+            else:
+                # No strict convergence, check for tol based convergence.
+                center_shift_tot = (center_shift**2).sum()
+                if center_shift_tot <= tol:
+                    if verbose:
+                        print(f"Converged at iteration {i}: center shift "
+                              f"{center_shift_tot} within tolerance {tol}.")
+                    break
+
+            labels_old[:] = labels
 
-        if center_shift_tot > 0:
+        if not strict_convergence:
             # rerun E-step so that predicted labels match cluster centers
             lloyd_iter(X, sample_weight, x_squared_norms, centers, centers,
                        weight_in_clusters, labels, center_shift, n_threads,
@@ -617,8 +644,8 @@ class KMeans(TransformerMixin, ClusterMixin, BaseEstimator):
         Relative tolerance with regards to Frobenius norm of the difference
         in the cluster centers of two consecutive iterations to declare
         convergence.
-        It's not advised to set `tol=0` since convergence might never be
-        declared due to rounding errors. Use a very small number instead.
+        If set to 0, convergence will be declared when the labels are the
+        same between two consecutive iterations.
 
     precompute_distances : {'auto', True, False}, default='auto'
         Precompute distances (faster but takes more memory).
diff --git a/sklearn/cluster/tests/test_k_means.py b/sklearn/cluster/tests/test_k_means.py
index 945db12363bec..d8f7a6f7ab7f2 100644
--- a/sklearn/cluster/tests/test_k_means.py
+++ b/sklearn/cluster/tests/test_k_means.py
@@ -136,7 +136,7 @@ def test_relocate_empty_clusters(array_constr):
 @pytest.mark.parametrize("distribution", ["normal", "blobs"])
 @pytest.mark.parametrize("array_constr", [np.array, sp.csr_matrix],
                          ids=["dense", "sparse"])
-@pytest.mark.parametrize("tol", [1e-2, 1e-4, 1e-8])
+@pytest.mark.parametrize("tol", [1e-2, 1e-8, 1e-100, 0])
 def test_kmeans_elkan_results(distribution, array_constr, tol):
     # Check that results are identical between lloyd and elkan algorithms
     rnd = np.random.RandomState(0)
@@ -160,21 +160,18 @@ def test_kmeans_elkan_results(distribution, array_constr, tol):
     assert km_elkan.inertia_ == pytest.approx(km_full.inertia_, rel=1e-6)
 
 
-@pytest.mark.parametrize("algorithm", ["full", "elkan"])
-def test_kmeans_convergence(algorithm):
+def test_kmeans_convergence():
     # Check that KMeans stops when convergence is reached when tol=0. (#16075)
-    # We can only ensure that if the number of threads is not to large,
-    # otherwise the roundings errors coming from the unpredictability of
-    # the order in which chunks are processed make the convergence criterion
-    # to never be exactly 0.
     rnd = np.random.RandomState(0)
     X = rnd.normal(size=(5000, 10))
+    max_iter = 300
 
-    with threadpool_limits(limits=1, user_api="openmp"):
-        km = KMeans(algorithm=algorithm, n_clusters=5, random_state=0,
-                    n_init=1, tol=0, max_iter=300).fit(X)
+    km_full = KMeans(algorithm="full", n_clusters=5, random_state=0,
+                     n_init=1, tol=0, max_iter=max_iter).fit(X)
+    km_elkan = KMeans(algorithm="elkan", n_clusters=5, random_state=0,
+                      n_init=1, tol=0, max_iter=max_iter).fit(X)
 
-    assert km.n_iter_ < 300
+    assert km_full.n_iter_ == km_elkan.n_iter_ < max_iter
 
 
 def test_minibatch_update_consistency():


Code review comments:
diff_hunk: @@ -496,15 +513,25 @@ def _kmeans_single_lloyd(X, sample_weight, centers_init, max_iter=300,
 
             centers, centers_new = centers_new, centers
 
-            center_shift_tot = (center_shift**2).sum()
-            if center_shift_tot <= tol:
+            if np.array_equal(labels, labels_old):
+                # First check the labels for strict convergence because
+                # center_shift might not be exactly 0 due to rounding errors.
path: sklearn/cluster/_kmeans.py
text: @user1:
This comment needs to be updated since this branch is now called even if `tol > 0`.

diff_hunk: @@ -617,8 +644,8 @@ class KMeans(TransformerMixin, ClusterMixin, BaseEstimator):
         Relative tolerance with regards to Frobenius norm of the difference
         in the cluster centers of two consecutive iterations to declare
         convergence.
-        It's not advised to set `tol=0` since convergence might never be
-        declared due to rounding errors. Use a very small number instead.
+        If set to 0, convergence will be declared when the labels are the
path: sklearn/cluster/_kmeans.py
text: @user1:
If set to 0 or very close to zero, convergence will be declared ...

@author:
I think this part is just not relevant anymore. When labels are the same, it implies that centers did not move so the norm of the diff is 0 <= tol. The label check is just a technicality for us to work around rounding errors.

diff_hunk: @@ -395,14 +398,25 @@ def _kmeans_single_elkan(X, sample_weight, centers_init, max_iter=300,
 
         centers, centers_new = centers_new, centers
 
-        center_shift_tot = (center_shift**2).sum()
-        if center_shift_tot <= tol:
+        if np.array_equal(labels, labels_old):
+            # First check the labels for strict convergence because
+            # center_shift might not be exactly 0 due to rounding errors.
path: sklearn/cluster/_kmeans.py
text: @user1:
This comment needs to be updated since this branch is now called even if `tol > 0`.