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:
@yigalagam You can use the environment variables `OMP_NUM_THREADS=1` before executing the script. It will use a single thread.

I assume what we observe here are rounding errors that are not reproducible due to the parallelism. @jeremiedbb would know more in details if this is the issue.

Labels: Bug