scikit-learn · lorentzenchr · Dec 22, 2021 · Dec 22, 2021 · Dec 22, 2021 · Dec 22, 2021
diff --git a/doc/whats_new/v1.1.rst b/doc/whats_new/v1.1.rst
@@ -311,6 +311,14 @@ Changelog
    now validate input parameters in `fit` instead of `__init__`.
    :pr:`21880` by :user:`Mrinal Tyagi <MrinalTyagi>`.
 
+:mod:`sklearn.mixture`
+......................
+
+- |Fix| Fix a bug that correctly initialize `precisions_cholesky_` in
+  :class:`mixture.GaussianMixture` when providing `precisions_init` by taking
+  its square root.
+  :pr:`22058` by :user:`Guillaume Lemaitre <glemaitre>`.
+
 :mod:`sklearn.pipeline`
 .......................
 

diff --git a/sklearn/mixture/_gaussian_mixture.py b/sklearn/mixture/_gaussian_mixture.py
@@ -728,7 +728,7 @@ def _initialize(self, X, resp):
                 self.precisions_init, lower=True
             )
         else:
-            self.precisions_cholesky_ = self.precisions_init
+            self.precisions_cholesky_ = np.sqrt(self.precisions_init)
 
     def _m_step(self, X, log_resp):
         """M step.

diff --git a/sklearn/mixture/tests/test_gaussian_mixture.py b/sklearn/mixture/tests/test_gaussian_mixture.py
@@ -11,6 +11,7 @@
 import numpy as np
 from scipy import stats, linalg
 
+from sklearn.cluster import KMeans
 from sklearn.covariance import EmpiricalCovariance
 from sklearn.datasets import make_spd_matrix
 from io import StringIO
@@ -21,6 +22,7 @@
     _estimate_gaussian_covariances_tied,
     _estimate_gaussian_covariances_diag,
     _estimate_gaussian_covariances_spherical,
+    _estimate_gaussian_parameters,
     _compute_precision_cholesky,
     _compute_log_det_cholesky,
 )
@@ -1241,6 +1243,7 @@ def test_gaussian_mixture_setting_best_params():
         random_state=rnd,
         n_components=len(weights_init),
         precisions_init=precisions_init,
+        max_iter=1,
     )
     # ensure that no error is thrown during fit
     gmm.fit(X)
@@ -1258,3 +1261,64 @@ def test_gaussian_mixture_setting_best_params():
         "lower_bound_",
     ]:
         assert hasattr(gmm, attr)
+
+
+def test_gaussian_mixture_precisions_init_diag():
+    """Check that we properly initialize `precision_cholesky_` when we manually
+    provide the precision matrix.
+
+    In this regard, we check the consistency between estimating the precision
+    matrix and providing the same precision matrix as initialization. It should
+    lead to the same results with the same number of iterations.
+
+    If the initialization is wrong then the number of iterations will increase.
+
+    Non-regression test for:
+    https://github.com/scikit-learn/scikit-learn/issues/16944
+    """
+    # generate a toy dataset
+    n_samples = 300
+    rng = np.random.RandomState(0)
+    shifted_gaussian = rng.randn(n_samples, 2) + np.array([20, 20])
+    C = np.array([[0.0, -0.7], [3.5, 0.7]])
+    stretched_gaussian = np.dot(rng.randn(n_samples, 2), C)
+    X = np.vstack([shifted_gaussian, stretched_gaussian])
+
+    # common parameters to check the consistency of precision initialization
+    n_components, covariance_type, reg_covar, random_state = 2, "diag", 1e-6, 0
+
+    # execute the manual initialization to compute the precision matrix:
+    # - run KMeans to have an initial guess
+    # - estimate the covariance
+    # - compute the precision matrix from the estimated covariance
+    resp = np.zeros((X.shape[0], n_components))
+    label = (
+        KMeans(n_clusters=n_components, n_init=1, random_state=random_state)
+        .fit(X)
+        .labels_
+    )
+    resp[np.arange(X.shape[0]), label] = 1
+    _, _, covariance = _estimate_gaussian_parameters(
+        X, resp, reg_covar=reg_covar, covariance_type=covariance_type
+    )
+    precisions_init = 1 / covariance
+
+    gm_with_init = GaussianMixture(
+        n_components=n_components,
+        covariance_type=covariance_type,
+        reg_covar=reg_covar,
+        precisions_init=precisions_init,
+        random_state=random_state,
+    ).fit(X)
+
+    gm_without_init = GaussianMixture(
+        n_components=n_components,
+        covariance_type=covariance_type,
+        reg_covar=reg_covar,
+        random_state=random_state,
+    ).fit(X)
+
+    assert gm_without_init.n_iter_ == gm_with_init.n_iter_
+    assert_allclose(
+        gm_with_init.precisions_cholesky_, gm_without_init.precisions_cholesky_
+    )