Gaussian Mixture: Diagonal covariance vectors might contain unreasonably negative values when the input datatype is np.float32

### Describe the bug

The Gaussian Mixture implementation shows numerical instabilities on single-precision floating point input numbers, that even large values of the regularization parameter reg_covar (like 0.1) cannot mitigate.

More specifically, diagonal covariance elements must not be negative. However, due to the numerical instabilities intrinsic to floating point arithmetic, they might end up being tiny negative numbers that reg_covar must compensate.
It turns out that, for some input float32 , the covariance can reach the unreasonable value of -0.99999979.
This is because squaring float32 numbers significantly magnifies their precision errors.

The proposed solution consists in converting float32 values to float64 before squaring them.
Care must be taken to not increase memory consumption in the overall process.
Hence, as avgX_means is equal to avg_means2, the return value can be simplified.

### Steps/Code to Reproduce

```python
import numpy as np
from sklearn.mixture import GaussianMixture

model = GaussianMixture(n_components=2, covariance_type="spherical", reg_covar=0.1)
model.fit(np.array([[9999.0], [0.0]], dtype=np.float32))
model.covariances_
```

### Expected Results

```python
array([0.1, 0.1])
```

### Actual Results

```python
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Input In [132], in <cell line: 49>()
     45 skgm._estimate_gaussian_covariances_diag = _optimized_estimate_gaussian_covariances_diag
     48 model = GaussianMixture(n_components=2,covariance_type="spherical", reg_covar=0.1)
---> 49 model.fit(np.array([[9999.0], [0.0]], dtype=np.float32))
     50 model.covariances_

File [~\AppData\Local\Programs\Python\Python39\lib\site-packages\sklearn\mixture\_base.py:200](http://localhost:8888/lab/tree/~/AppData/Local/Programs/Python/Python39/lib/site-packages/sklearn/mixture/_base.py#line=199), in BaseMixture.fit(self, X, y)
    174 def fit(self, X, y=None):
    175     """Estimate model parameters with the EM algorithm.
    176 
    177     The method fits the model ``n_init`` times and sets the parameters with
   (...)
    198         The fitted mixture.
    199     """
--> 200     self.fit_predict(X, y)
    201     return self

File [~\AppData\Local\Programs\Python\Python39\lib\site-packages\sklearn\mixture\_base.py:253](http://localhost:8888/lab/tree/~/AppData/Local/Programs/Python/Python39/lib/site-packages/sklearn/mixture/_base.py#line=252), in BaseMixture.fit_predict(self, X, y)
    250 self._print_verbose_msg_init_beg(init)
    252 if do_init:
--> 253     self._initialize_parameters(X, random_state)
    255 lower_bound = -np.inf if do_init else self.lower_bound_
    257 if self.max_iter == 0:

File [~\AppData\Local\Programs\Python\Python39\lib\site-packages\sklearn\mixture\_base.py:160](http://localhost:8888/lab/tree/~/AppData/Local/Programs/Python/Python39/lib/site-packages/sklearn/mixture/_base.py#line=159), in BaseMixture._initialize_parameters(self, X, random_state)
    155 else:
    156     raise ValueError(
    157         "Unimplemented initialization method '%s'" % self.init_params
    158     )
--> 160 self._initialize(X, resp)

File [~\AppData\Local\Programs\Python\Python39\lib\site-packages\sklearn\mixture\_gaussian_mixture.py:723](http://localhost:8888/lab/tree/~/AppData/Local/Programs/Python/Python39/lib/site-packages/sklearn/mixture/_gaussian_mixture.py#line=722), in GaussianMixture._initialize(self, X, resp)
    721 if self.precisions_init is None:
    722     self.covariances_ = covariances
--> 723     self.precisions_cholesky_ = _compute_precision_cholesky(
    724         covariances, self.covariance_type
    725     )
    726 elif self.covariance_type == "full":
    727     self.precisions_cholesky_ = np.array(
    728         [
    729             linalg.cholesky(prec_init, lower=True)
    730             for prec_init in self.precisions_init
    731         ]
    732     )

File [~\AppData\Local\Programs\Python\Python39\lib\site-packages\sklearn\mixture\_gaussian_mixture.py:347](http://localhost:8888/lab/tree/~/AppData/Local/Programs/Python/Python39/lib/site-packages/sklearn/mixture/_gaussian_mixture.py#line=346), in _compute_precision_cholesky(covariances, covariance_type)
    345 else:
    346     if np.any(np.less_equal(covariances, 0.0)):
--> 347         raise ValueError(estimate_precision_error_message)
    348     precisions_chol = 1.0 / np.sqrt(covariances)
    349 return precisions_chol

ValueError: Fitting the mixture model failed because some components have ill-defined empirical covariance (for instance caused by singleton or collapsed samples). Try to decrease the number of components, or increase reg_covar.
```

### Versions

```shell
System:
    python: 3.9.13 (tags/v3.9.13:6de2ca5, May 17 2022, 16:36:42) [MSC v.1929 64 bit (AMD64)]
executable: C:\Users\leonce.mekinda\AppData\Local\Programs\Python\Python39\python.exe
   machine: Windows-10-10.0.19043-SP0

Python dependencies:
      sklearn: 1.1.1
          pip: 22.0.4
   setuptools: 58.1.0
        numpy: 1.22.4
        scipy: 1.8.1
       Cython: 0.29.30
       pandas: 1.4.3
   matplotlib: 3.5.2
       joblib: 1.1.0
threadpoolctl: 3.1.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: C:\Users\leonce.mekinda\AppData\Local\Programs\Python\Python39\Lib\site-packages\numpy\.libs\libopenblas.EL2C6PLE4ZYW3ECEVIV3OXXGRN2NRFM2.gfortran-win_amd64.dll
        version: 0.3.18
threading_layer: pthreads
   architecture: Haswell
    num_threads: 12

       user_api: openmp
   internal_api: openmp
         prefix: vcomp
       filepath: C:\Users\leonce.mekinda\AppData\Local\Programs\Python\Python39\Lib\site-packages\sklearn\.libs\vcomp140.dll
        version: None
    num_threads: 12

       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: C:\Users\leonce.mekinda\AppData\Local\Programs\Python\Python39\Lib\site-packages\scipy\.libs\libopenblas.XWYDX2IKJW2NMTWSFYNGFUWKQU3LYTCZ.gfortran-win_amd64.dll
        version: 0.3.17
threading_layer: pthreads
   architecture: Haswell
    num_threads: 12
```

COMMENT:
Thanks for the report.

For information, I updated the description is use the markdown delimiter for the verbatim code blocks to make it more readable.

It seems that in this case, we use `_estimate_gaussian_covariances_spherical` that directly calls into `_estimate_gaussian_covariances_diag` which does the following:

https://github.com/scikit-learn/scikit-learn/blob/8ded7f43f8141959ac267450c126014d4c935907/sklearn/mixture/_gaussian_mixture.py#L209-L232

So there is indeed a possibility for [catastrophic cancellation](https://en.wikipedia.org/wiki/Catastrophic_cancellation).

Converting to float64 might help, but if one we want to adapt this code to work with array API on GPU devices, some devices might not support float64 operation (e.g. the built-in GPUs of macOS machines).

However the code looks like:

```
(w a - b)^2 = w^2 a^2 - 2 w a b + b^
```

and I was wondering if computing the left hand side instead of the right hand size wouldn't be more numerically stable. But this does not actually match. In particular this is not exactly the same as:

```python
return ((resp.T @ X) / nk[:, np.newaxis] - means) ** 2 + reg_covar
```

because the responsibility weights are not simple scalar weights. So maybe it's not possible to do the factored computation.

Labels: Bug, Numerical Stability