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:
> > In this example, you still use np.float64. Just use float
> 
> I would rather not use implicit precision levels in the source code. But, in this particular case, this does not change anything: Python `float` are promoted to `float64` both by NumPy and Pytorch:
> 
> ```python
> >>> import torch
> >>> torch.ones(5, dtype=float, device="mps")
> Traceback (most recent call last):
>   Cell In[5], line 1
>     torch.ones(5, dtype=float, device="mps")
> TypeError: Cannot convert a MPS Tensor to float64 dtype as the MPS framework doesn't support float64. Please use float32 instead.
> ```
> 
> And even putting potential array API / torch support aside, I would rather avoid duplicating a datastructure of the size of X to change it's dtype from float32 to float64 to avoid triping the memory usage when fitting GMM on large data if we can.

I guess on torch you will use torch.float that maps to torch.float32.

Regarding the array duplication, the original code does X * X, which already creates a new array. What the proposed code does is to create a new array and then square values **in-place**

Labels: Bug, Numerical Stability