Inconsistent results with same random seed

### Describe the bug

I'm not sure this is a bug and I couldn't find anything in the issues archive, but I'm seeing an inconsistency between consecutive runs of kmeans.fit, even when setting the same random seed via np.random.seed or random_state. I pinned it down to multithreading: if I force self._n_threads to 1 inside _kmeans.py/fit the inconsistency goes away. See test code below. Is this expected behavior, and is there a way to eliminate multithreading during fitting?

```python
from sklearn.cluster import KMeans
import numpy as np

np.random.seed(0)
matrix = np.random.randn(1000, 100)

for i in range(100):
    k1 = KMeans(n_clusters=5, max_iter=1, random_state=0, n_init=1)
    k1.fit(matrix)
    k2 = KMeans(n_clusters=5, max_iter=1, random_state=0, n_init=1)
    k2.fit(matrix)
    print('min diff: %s, max diff: %s' % (str(np.min(k1.cluster_centers_ - k2.cluster_centers_)), str(np.max(k1.cluster_centers_ - k2.cluster_centers_))))
```

### Steps/Code to Reproduce

```python
from sklearn.cluster import KMeans
import numpy as np

np.random.seed(0)
matrix = np.random.randn(1000, 100)

for i in range(100):
    k1 = KMeans(n_clusters=5, max_iter=1, random_state=0, n_init=1)
    k1.fit(matrix)
    k2 = KMeans(n_clusters=5, max_iter=1, random_state=0, n_init=1)
    k2.fit(matrix)
    print('min diff: %s, max diff: %s' % (str(np.min(k1.cluster_centers_ - k2.cluster_centers_)), str(np.max(k1.cluster_centers_ - k2.cluster_centers_))))
```

### Expected Results

min diff: 0.0, max diff: 0.0
min diff: 0.0, max diff: 0.0
min diff: 0.0, max diff: 0.0
min diff: 0.0, max diff: 0.0
min diff: 0.0, max diff: 0.0
min diff: 0.0, max diff: 0.0
...

### Actual Results

min diff: -1.1102230246251565e-16, max diff: 2.220446049250313e-16
min diff: -2.220446049250313e-16, max diff: 1.1102230246251565e-16
min diff: -1.1102230246251565e-16, max diff: 1.1102230246251565e-16
min diff: 0.0, max diff: 0.0
min diff: -1.1102230246251565e-16, max diff: 1.1102230246251565e-16
min diff: -1.1102230246251565e-16, max diff: 2.220446049250313e-16
,,,

### Versions

```shell
System:
    python: 3.8.10 (tags/v3.8.10:3d8993a, May  3 2021, 11:48:03) [MSC v.1928 64 bit (AMD64)]
executable: C:\Python38\python.exe
   machine: Windows-10-10.0.19045-SP0
Python dependencies:
      sklearn: 1.3.1
          pip: 21.1.1
   setuptools: 56.0.0
        numpy: 1.24.2
        scipy: 1.8.0
       Cython: 0.29.28
       pandas: 1.3.4
   matplotlib: 3.4.3
       joblib: 1.3.2
threadpoolctl: 3.1.0
Built with OpenMP: True
threadpoolctl info:
       user_api: openmp
   internal_api: openmp
         prefix: vcomp
       filepath: C:\Python38\Lib\site-packages\sklearn\.libs\vcomp140.dll
        version: None
    num_threads: 16
       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: C:\Python38\Lib\site-packages\numpy\.libs\libopenblas64__v0.3.21-gcc_10_3_0.dll
        version: 0.3.21
threading_layer: pthreads
   architecture: SkylakeX
    num_threads: 16
       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: C:\Python38\Lib\site-packages\scipy\.libs\libopenblas.XWYDX2IKJW2NMTWSFYNGFUWKQU3LYTCZ.gfortran-win_amd64.dll
        version: 0.3.17
threading_layer: pthreads
   architecture: Prescott
    num_threads: 16
```

COMMENT:
Workaround for consistent results using a modified version of KMeans fit where self._n_threads is set to 1 (replaces line 1488 in _kmeans.py, scikit-learn 1.3.1):

```
import warnings
import numpy as np
import scipy.sparse as sp
from sklearn.cluster import KMeans
from sklearn.utils import check_array, check_random_state
from sklearn.utils.validation import _check_sample_weight, _is_arraylike_not_scalar
from sklearn.utils.extmath import row_norms
from sklearn.cluster._kmeans import _kmeans_single_lloyd, _kmeans_single_elkan
from sklearn.cluster._k_means_common import _is_same_clustering
from sklearn.exceptions import ConvergenceWarning

def kmeans_fit(kmeans_obj, X, sample_weight=None):
    """Compute k-means clustering.

    Parameters
    ----------
    X : {array-like, sparse matrix} of shape (n_samples, n_features)
        Training instances to cluster. It must be noted that the data
        will be converted to C ordering, which will cause a memory
        copy if the given data is not C-contiguous.
        If a sparse matrix is passed, a copy will be made if it's not in
        CSR format.

    sample_weight : array-like of shape (n_samples,), default=None
        The weights for each observation in X. If None, all observations
        are assigned equal weight.

        .. versionadded:: 0.20

    Returns
    -------
    self : object
        Fitted estimator.
    """
    X = kmeans_obj._validate_data(
        X,
        accept_sparse="csr",
        dtype=[np.float64, np.float32],
        order="C",
        copy=kmeans_obj.copy_x,
        accept_large_sparse=False,
    )

    kmeans_obj._check_params(X)
    random_state = check_random_state(kmeans_obj.random_state)
    sample_weight = _check_sample_weight(sample_weight, X, dtype=X.dtype)
    kmeans_obj._n_threads = 1 # disable parallelization for consistent clustering

    # Validate init array
    init = kmeans_obj.init
    init_is_array_like = _is_arraylike_not_scalar(init)
    if init_is_array_like:
        init = check_array(init, dtype=X.dtype, copy=True, order="C")
        kmeans_obj._validate_center_shape(X, init)

    # subtract of mean of x for more accurate distance computations
    if not sp.issparse(X):
        X_mean = X.mean(axis=0)
        # The copy was already done above
        X -= X_mean

        if init_is_array_like:
            init -= X_mean

    # precompute squared norms of data points
    x_squared_norms = row_norms(X, squared=True)

    if kmeans_obj._algorithm == "elkan":
        kmeans_single = _kmeans_single_elkan
    else:
        kmeans_single = _kmeans_single_lloyd
        kmeans_obj._check_mkl_vcomp(X, X.shape[0])

    best_inertia, best_labels = None, None

    for i in range(kmeans_obj._n_init):
        # Initialize centers
        centers_init = kmeans_obj._init_centroids(
            X, x_squared_norms=x_squared_norms, init=init, random_state=random_state
        )
        if kmeans_obj.verbose:
            print("Initialization complete")

        # run a k-means once
        labels, inertia, centers, n_iter_ = kmeans_single(
            X,
            sample_weight,
            centers_init,
            max_iter=kmeans_obj.max_iter,
            verbose=kmeans_obj.verbose,
            tol=kmeans_obj._tol,
            x_squared_norms=x_squared_norms,
            n_threads=kmeans_obj._n_threads,
        )

        # determine if these results are the best so far
        # we chose a new run if it has a better inertia and the clustering is
        # different from the best so far (it's possible that the inertia is
        # slightly better even if the clustering is the same with potentially
        # permuted labels, due to rounding errors)
        if best_inertia is None or (
            inertia < best_inertia
            and not _is_same_clustering(labels, best_labels, kmeans_obj.n_clusters)
        ):
            best_labels = labels
            best_centers = centers
            best_inertia = inertia
            best_n_iter = n_iter_

    if not sp.issparse(X):
        if not kmeans_obj.copy_x:
            X += X_mean
        best_centers += X_mean

    distinct_clusters = len(set(best_labels))
    if distinct_clusters < kmeans_obj.n_clusters:
        warnings.warn(
            "Number of distinct clusters ({}) found smaller than "
            "n_clusters ({}). Possibly due to duplicate points "
            "in X.".format(distinct_clusters, kmeans_obj.n_clusters),
            ConvergenceWarning,
            stacklevel=2,
        )

    kmeans_obj.cluster_centers_ = best_centers
    kmeans_obj._n_features_out = kmeans_obj.cluster_centers_.shape[0]
    kmeans_obj.labels_ = best_labels
    kmeans_obj.inertia_ = best_inertia
    kmeans_obj.n_iter_ = best_n_iter
    return kmeans_obj
```

Labels: Bug