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:
Thank you for taking care of it.
Using float (alias to np.float64) instead of float64 would be supported on 32-bit platforms as well, and provide the floating point double-precision of all the other computations on 64-bit platforms:
```python
    # Convert to a new array of float or float64 depending on the platform
    X2 = X.astype(float)
    # In-place squaring
    X2 **= 2 
    avg_X2 = np.dot(resp.T, X2) / nk[:, np.newaxis]
    avg_means2 = means**2
    return np.clip(avg_X2 - avg_means2, 0, None) + reg_covar
```
Note that avg_X_means is exactly avg_means2. This is why this term can be removed and the returned value simplified 
to return np.clip(avg_X2 - avg_means2, 0, None )+ reg_covar

Returning  ((resp.T @ X) / nk[:, np.newaxis] - means) ** 2 would not be correct because the responsibility matrix must not be squared.

Labels: Bug, Numerical Stability