GaussianMixture log-probabilities are numerically inaccurate

[JohannesBuchner]

Describe the bug

While building an extension to (already fitted) GaussianMixture models (https://github.com/JohannesBuchner/askcarl/),
I was using https://hypothesis.readthedocs.io/ for testing and came across numerical inaccuracies in the log-probabilities computed by scikit-learn's GaussianMixture. These already occur with single components, where one can take scipy.stats.multivariate_normal.logpdf as the ground truth.

The inaccuracies appear in the function _estimate_log_gaussian_prob.

The log-probabilities can be off by 0.2 (see the very last example), which really is not small if those probabilities are used for likelihoods.

I uploaded the full script at https://github.com/JohannesBuchner/gmm-tests/blob/main/test_sklearn.py but the important excerpts are below. A unrelated issue I have is that I want to build sklearn.mixture.GaussianMixture from scratch without fitting it, but .score and .predict_proba always give 0 and 1, respectively. Is there a initialization step I am missing to be allowed to use these functions?

EDIT (by @ogrisel): here is a copy of the reproducer (to make this report self-contained):

Details

import numpy as np
from numpy import array
from scipy.stats import multivariate_normal
from scipy.special import logsumexp
from numpy.testing import assert_allclose
from hypothesis import given, strategies as st, example, settings
from hypothesis.extra.numpy import arrays
import sklearn.mixture
from  sklearn.mixture._gaussian_mixture import _estimate_log_gaussian_prob

def valid_QR(vectors):
    q, r = np.linalg.qr(vectors)
    return q.shape == vectors.shape and np.all(np.abs(np.diag(r)) > 1e-3) and np.all(np.abs(np.diag(r)) < 1000)

def make_covariance_matrix_via_QR(normalisations, vectors):
    q, r = np.linalg.qr(vectors)
    orthogonal_vectors = q @ np.diag(np.diag(r))
    cov = orthogonal_vectors @ np.diag(normalisations) @ orthogonal_vectors.T
    return cov

def valid_covariance_matrix(A, min_std=1e-6):
    if not np.isfinite(A).all():
        return False
    #if not np.std(A) > min_std:
    #    return False
    if (np.diag(A) <= min_std).any():
        return False

    try:
        np.linalg.inv(np.linalg.inv(A))
    except np.linalg.LinAlgError:
        return False

    try:
        multivariate_normal(mean=np.zeros(len(A)), cov=A)
    except ValueError:
        return False

    return True

# Strategy to generate arbitrary dimensionality mean and covariance
@st.composite
def mean_and_cov(draw):
    dim = draw(st.integers(min_value=1, max_value=10))  # Arbitrary dimensionality
    mu = draw(arrays(np.float64, (dim,), elements=st.floats(-10, 10)))  # Mean vector
    eigval = draw(arrays(np.float64, (dim,), elements=st.floats(1e-6, 10)))
    vectors = draw(arrays(np.float64, (dim,dim), elements=st.floats(-10, 10)).filter(valid_QR))
    cov = make_covariance_matrix_via_QR(eigval, vectors)
    x = draw(arrays(np.float64, (2, dim), elements=st.floats(-1e6, 1e6)))
    return dim, mu, cov, x


@settings(max_examples=1000)
@given(mean_and_cov())
@example(
    mean_cov=(2, array([0., 0.]), array([[1., 0.],
            [0., 1.]]), array([[0., 0.],
            [0., 0.]])),
).via('discovered failure')
@example(
    mean_cov=(1, array([0.]), array([[1.]]), array([[1.],
            [0.]])),
).via('discovered failure')
@example(
    mean_cov=(6,
     array([0., 0., 0., 0., 0., 0.]),
     array([[ 1.00000107e+00, -5.71254024e-08,  9.99999696e-01,
              9.99999583e-01,  9.99999675e-01,  9.99999979e-01],
            [-5.71254024e-08,  1.54534373e-06,  1.67651586e-08,
             -2.54656267e-07, -4.52842673e-09,  2.99544937e-07],
            [ 9.99999696e-01,  1.67651586e-08,  1.00000095e+00,
              9.99999657e-01,  9.99999707e-01,  9.99999994e-01],
            [ 9.99999583e-01, -2.54656267e-07,  9.99999657e-01,
              1.00000119e+00,  9.99999635e-01,  9.99999940e-01],
            [ 9.99999675e-01, -4.52842673e-09,  9.99999707e-01,
              9.99999635e-01,  1.00000099e+00,  9.99999990e-01],
            [ 9.99999979e-01,  2.99544937e-07,  9.99999994e-01,
              9.99999940e-01,  9.99999990e-01,  1.00000010e+00]]),
     array([[0., 0., 0., 0., 0., 0.],
            [0., 0., 0., 0., 0., 0.]])),
).via('discovered failure')
@example(
    mean_cov=(6,
     array([0., 0., 0., 0., 0., 0.]),
     array([[ 2.28000099, -4.12000012,  2.27999939, -1.59999953,  2.27999988,
              2.27999986],
            [-4.12000012, 21.48000008, -4.12000007,  6.39999967, -4.11999995,
             -4.11999995],
            [ 2.27999939, -4.12000007,  2.28000087, -1.59999973,  2.27999991,
              2.2799999 ],
            [-1.59999953,  6.39999967, -1.59999973,  2.0000013 , -1.60000021,
             -1.6000002 ],
            [ 2.27999988, -4.11999995,  2.27999991, -1.60000021,  2.28000032,
              2.27999984],
            [ 2.27999986, -4.11999995,  2.2799999 , -1.6000002 ,  2.27999984,
              2.28000035]]),
     array([[0., 0., 0., 0., 0., 0.],
            [0., 0., 0., 0., 0., 0.]])),
).via('discovered failure')
@example(
    mean_cov=(7,
     array([0., 0., 0., 0., 0., 0., 0.]),
     array([[156.25637361, -31.26969455,  12.40439959, -31.24564618,
              -31.22989933, -31.24225165, -31.2688804 ],
             [-31.26969455,   6.33006684,  -2.20457805,   6.2198964 ,
                6.16861271,   6.20884123,   6.34227887],
             [ 12.40439959,  -2.20457805,   2.43400622,  -2.56530348,
               -2.80150635,  -2.61622149,  -2.21679021],
             [-31.24564618,   6.2198964 ,  -2.56530348,   7.16206397,
                5.80788728,   5.8481158 ,   6.20768424],
             [-31.22989933,   6.16861271,  -2.80150635,   5.80788728,
                7.34003101,   5.75696927,   6.15640055],
             [-31.24225165,   6.20884123,  -2.61622149,   5.8481158 ,
                5.75696927,   7.23169779,   6.19662906],
             [-31.2688804 ,   6.34227887,  -2.21679021,   6.20768424,
                6.15640055,   6.19662906,   6.36588918]]),
     array([[0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0.]])),
).via('discovered failure')
@example(
    mean_cov=(1, array([0.]), array([[1.]]), array([[0.], [0.]])),
).via('discovered failure')
@example(
    mean_cov=(8,
     array([0., 0., 0., 0., 0., 0., 0., 0.]),
     array([[ 1.53365164e-06, -4.56057424e-07, -2.47157040e-07,
              -1.88188005e-07,  1.47416388e-07,  1.47416391e-07,
               4.49153299e-07,  1.47416391e-07],
             [-4.56057424e-07,  1.76562629e+00,  9.99999625e-01,
               9.99999679e-01, -3.26562499e+00,  1.98437501e+00,
               2.53124938e+00,  1.98437501e+00],
             [-2.47157040e-07,  9.99999625e-01,  1.00000099e+00,
               9.99999580e-01,  9.99999915e-01,  9.99999915e-01,
               1.00000006e+00,  9.99999915e-01],
             [-1.88188005e-07,  9.99999679e-01,  9.99999580e-01,
               1.00000091e+00,  9.99999927e-01,  9.99999927e-01,
               1.00000005e+00,  9.99999927e-01],
             [ 1.47416388e-07, -3.26562499e+00,  9.99999915e-01,
               9.99999927e-01,  2.47656250e+01, -4.48437497e+00,
              -7.53124995e+00, -4.48437497e+00],
             [ 1.47416391e-07,  1.98437501e+00,  9.99999915e-01,
               9.99999927e-01, -4.48437497e+00,  2.26562529e+00,
               2.96875004e+00,  2.26562479e+00],
             [ 4.49153299e-07,  2.53124938e+00,  1.00000006e+00,
               1.00000005e+00, -7.53124995e+00,  2.96875004e+00,
               4.06250038e+00,  2.96875004e+00],
             [ 1.47416391e-07,  1.98437501e+00,  9.99999915e-01,
               9.99999927e-01, -4.48437497e+00,  2.26562479e+00,
               2.96875004e+00,  2.26562529e+00]]),
     array([[0., 0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0., 0.]])),
).via('discovered failure')
@example(
    mean_cov=(7,
        array([0., 0., 0., 0., 0., 0., 0.]),
        array([[ 7.52684752, -4.9545883 , -0.31361835,  0.52269766,  0.52269929,
             6.92515709, -9.70648636],
           [-4.9545883 ,  3.26139136,  0.20644122, -0.3440688 , -0.34406688,
            -4.55852195,  6.38935566],
           [-0.31361835,  0.20644122,  0.01308347, -0.02177316, -0.0217791 ,
            -0.28854817,  0.40443697],
           [ 0.52269766, -0.3440688 , -0.02177316,  0.03630094,  0.0362984 ,
             0.48091352, -0.67406172],
           [ 0.52269929, -0.34406688, -0.0217791 ,  0.0362984 ,  0.03630217,
             0.48091543, -0.67405981],
           [ 6.92515709, -4.55852195, -0.28854817,  0.48091352,  0.48091543,
             6.37156541, -8.9305567 ],
           [-9.70648636,  6.38935566,  0.40443697, -0.67406172, -0.67405981,
            -8.9305567 , 12.51732134]]),
       array([[0., 0., 0., 0., 0., 0., 0.], [0., 0., 0., 0., 0., 0., 0.]]))
).via('discovered failure')
def test_single(mean_cov):
    # a askcarl with one component must behave the same as a single gaussian
    ndim, mu, cov, x = mean_cov
    ncomponents = 1
    if not valid_covariance_matrix(cov):
        return

    rv_truth = multivariate_normal(mu, cov)
    covs = [cov]
    means = [mu]
    weights = [1.0]
    assert x.shape == (2, ndim), (x.shape, 2, ndim)

    truth_logp = rv_truth.logpdf(x)
    assert truth_logp.shape == (len(x),)
    truth_p = rv_truth.pdf(x)

    precisions = [np.linalg.inv(cov) for cov in covs]

    # initialize GMM:
    skgmm = sklearn.mixture.GaussianMixture(
        n_components=ncomponents, weights_init=weights,
        means_init=means, precisions_init=precisions)
    skgmm._initialize(np.zeros((1, 1)), None)
    skgmm._set_parameters((weights, np.array(means), covs, skgmm.precisions_cholesky_))
    assert_allclose(skgmm.weights_, weights)
    assert_allclose(skgmm.means_, means)
    assert_allclose(skgmm.covariances_, covs)
    print('input:', x, skgmm.weights_, skgmm.means_, 'covs:', skgmm.covariances_, skgmm.precisions_cholesky_)
    print('expectation:', truth_p, 'log:', truth_logp)

    # compute probability:
    sk_logp = _estimate_log_gaussian_prob(x, skgmm.means_, skgmm.precisions_cholesky_, 'full').flatten()
    print('actual (_estimate_log_gaussian_prob):', sk_logp)
    assert sk_logp.shape == (len(x),), (sk_logp.shape, len(x))
    assert_allclose(truth_logp, sk_logp, atol=1e-4, rtol=1e-4)
    sk_logp1, sk_logp2 = skgmm._estimate_log_prob_resp(x)
    print('actual (_estimate_log_prob_resp):', sk_logp1, sk_logp2)
    assert_allclose(truth_logp, sk_logp1, atol=1e-4, rtol=1e-4)
    #assert_allclose(truth_logp, sk_logp2)

    sk_logp = skgmm.score(x)
    sk_p = skgmm.predict_proba(x)
    print('actual:', sk_p, 'log', sk_logp)

    # TODO: enable these two tests as well
    # assert_allclose(truth_logp, sk_logp, atol=1e-2, rtol=1e-2)  # -> sklearn-GMM gives score=0
    # assert_allclose(truth_p, sk_p, atol=1e-300, rtol=1e-4)  # -> sklearn-GMM gives p=1

Steps/Code to Reproduce

1. test setup, given means mu, covariances covs and its inversion precisions (values are listed further down).

    skgmm = sklearn.mixture.GaussianMixture(
        n_components=ncomponents, weights_init=weights,
        means_init=means, precisions_init=precisions)
    skgmm._initialize(np.zeros((1, 1)), None)
    skgmm._set_parameters((weights, np.array(means), covs, skgmm.precisions_cholesky_))

2. compute logprob:

    truth_logp = rv_truth.logpdf(x)   # expectation
    sk_logp = _estimate_log_gaussian_prob(x, skgmm.means_, skgmm.precisions_cholesky_, 'full').flatten()   # actual

3. test:

    print('actual (_estimate_log_gaussian_prob):', sk_logp)
    assert sk_logp.shape == (len(x),), (sk_logp.shape, len(x))
    assert_allclose(truth_logp, sk_logp, atol=1e-4, rtol=1e-4)
    sk_logp1, sk_logp2 = skgmm._estimate_log_prob_resp(x)
    print('actual (_estimate_log_prob_resp):', sk_logp1, sk_logp2)
    assert_allclose(truth_logp, sk_logp1, atol=1e-4, rtol=1e-4)

Expected Results

failing example:

    |     mean_cov=(6,
    |      array([0., 0., 0., 0., 0., 0.]),
    |      array([[ 2.28000099, -4.12000012,  2.27999939, -1.59999953,  2.27999988,
    |               2.27999986],
    |             [-4.12000012, 21.48000008, -4.12000007,  6.39999967, -4.11999995,
    |              -4.11999995],
    |             [ 2.27999939, -4.12000007,  2.28000087, -1.59999973,  2.27999991,
    |               2.2799999 ],
    |             [-1.59999953,  6.39999967, -1.59999973,  2.0000013 , -1.60000021,
    |              -1.6000002 ],
    |             [ 2.27999988, -4.11999995,  2.27999991, -1.60000021,  2.28000032,
    |               2.27999984],
    |             [ 2.27999986, -4.11999995,  2.2799999 , -1.6000002 ,  2.27999984,
    |               2.28000035]]),
    |      array([[0., 0., 0., 0., 0., 0.],
    |             [0., 0., 0., 0., 0., 0.]])),

where mean_cov is the tuple ndim, mu, cov, and points x where the density is evaluated.

fails at: assert_allclose(truth_logp, sk_logp, atol=1e-4, rtol=1e-4) with x=truth_logp, y=sk_logp.
| x: array([20.388033, 20.388033])

Actual Results

|  y: array([20.385592, 20.385592])

further failures:

    | Mismatched elements: 2 / 2 (100%)
    | Max absolute difference: 0.00596896
    | Max relative difference: 0.00098574
    |  x: array([6.061245, 6.061245])
    |  y: array([6.055276, 6.055276])
    | Falsifying explicit example: test_single(
    |     mean_cov=(7,
    |      array([0., 0., 0., 0., 0., 0., 0.]),
    |      array([[156.25637361, -31.26969455,  12.40439959, -31.24564618,
    |              -31.22989933, -31.24225165, -31.2688804 ],
    |             [-31.26969455,   6.33006684,  -2.20457805,   6.2198964 ,
    |                6.16861271,   6.20884123,   6.34227887],
    |             [ 12.40439959,  -2.20457805,   2.43400622,  -2.56530348,
    |               -2.80150635,  -2.61622149,  -2.21679021],
    |             [-31.24564618,   6.2198964 ,  -2.56530348,   7.16206397,
    |                5.80788728,   5.8481158 ,   6.20768424],
    |             [-31.22989933,   6.16861271,  -2.80150635,   5.80788728,
    |                7.34003101,   5.75696927,   6.15640055],
    |             [-31.24225165,   6.20884123,  -2.61622149,   5.8481158 ,
    |                5.75696927,   7.23169779,   6.19662906],
    |             [-31.2688804 ,   6.34227887,  -2.21679021,   6.20768424,
    |                6.15640055,   6.19662906,   6.36588918]]),
    |      array([[0., 0., 0., 0., 0., 0., 0.],
    |             [0., 0., 0., 0., 0., 0., 0.]])),

and

    | Mismatched elements: 2 / 2 (100%)
    | Max absolute difference: 0.04808096
    | Max relative difference: 0.00142038
    |  x: array([33.89883, 33.89883])
    |  y: array([33.850749, 33.850749])
    | Falsifying explicit example: test_single(
    |     mean_cov=(8,
    |      array([0., 0., 0., 0., 0., 0., 0., 0.]),
    |      array([[ 1.53365164e-06, -4.56057424e-07, -2.47157040e-07,
    |              -1.88188005e-07,  1.47416388e-07,  1.47416391e-07,
    |               4.49153299e-07,  1.47416391e-07],
    |             [-4.56057424e-07,  1.76562629e+00,  9.99999625e-01,
    |               9.99999679e-01, -3.26562499e+00,  1.98437501e+00,
    |               2.53124938e+00,  1.98437501e+00],
    |             [-2.47157040e-07,  9.99999625e-01,  1.00000099e+00,
    |               9.99999580e-01,  9.99999915e-01,  9.99999915e-01,
    |               1.00000006e+00,  9.99999915e-01],
    |             [-1.88188005e-07,  9.99999679e-01,  9.99999580e-01,
    |               1.00000091e+00,  9.99999927e-01,  9.99999927e-01,
    |               1.00000005e+00,  9.99999927e-01],
    |             [ 1.47416388e-07, -3.26562499e+00,  9.99999915e-01,
    |               9.99999927e-01,  2.47656250e+01, -4.48437497e+00,
    |              -7.53124995e+00, -4.48437497e+00],
    |             [ 1.47416391e-07,  1.98437501e+00,  9.99999915e-01,
    |               9.99999927e-01, -4.48437497e+00,  2.26562529e+00,
    |               2.96875004e+00,  2.26562479e+00],
    |             [ 4.49153299e-07,  2.53124938e+00,  1.00000006e+00,
    |               1.00000005e+00, -7.53124995e+00,  2.96875004e+00,
    |               4.06250038e+00,  2.96875004e+00],
    |             [ 1.47416391e-07,  1.98437501e+00,  9.99999915e-01,
    |               9.99999927e-01, -4.48437497e+00,  2.26562479e+00,
    |               2.96875004e+00,  2.26562529e+00]]),
    |      array([[0., 0., 0., 0., 0., 0., 0., 0.],
    |             [0., 0., 0., 0., 0., 0., 0., 0.]])),

and

    | Mismatched elements: 2 / 2 (100%)
    | Max absolute difference: 0.1614742
    | Max relative difference: 0.00484743
    |  x: array([33.149859, 33.149859])
    |  y: array([33.311333, 33.311333])
    | Falsifying explicit example: test_single(
    |     mean_cov=(7,
    |      array([0., 0., 0., 0., 0., 0., 0.]),
    |      array([[ 7.52684752, -4.9545883 , -0.31361835,  0.52269766,  0.52269929,
    |               6.92515709, -9.70648636],
    |             [-4.9545883 ,  3.26139136,  0.20644122, -0.3440688 , -0.34406688,
    |              -4.55852195,  6.38935566],
    |             [-0.31361835,  0.20644122,  0.01308347, -0.02177316, -0.0217791 ,
    |              -0.28854817,  0.40443697],
    |             [ 0.52269766, -0.3440688 , -0.02177316,  0.03630094,  0.0362984 ,
    |               0.48091352, -0.67406172],
    |             [ 0.52269929, -0.34406688, -0.0217791 ,  0.0362984 ,  0.03630217,
    |               0.48091543, -0.67405981],
    |             [ 6.92515709, -4.55852195, -0.28854817,  0.48091352,  0.48091543,
    |               6.37156541, -8.9305567 ],
    |             [-9.70648636,  6.38935566,  0.40443697, -0.67406172, -0.67405981,
    |              -8.9305567 , 12.51732134]]),
    |      array([[0., 0., 0., 0., 0., 0., 0.],
    |             [0., 0., 0., 0., 0., 0., 0.]])),

Versions

System:
    python: 3.12.3 (main, Sep 11 2024, 14:17:37) [GCC 13.2.0]
executable: /usr/bin/python3
   machine: Linux-6.8.0-45-generic-x86_64-with-glibc2.39

Python dependencies:
      sklearn: 1.5.0
          pip: 24.0
   setuptools: 68.1.2
        numpy: 1.26.3
        scipy: 1.14.1
       Cython: 3.0.11
       pandas: 2.2.2
   matplotlib: 3.9.2
       joblib: 1.4.2
threadpoolctl: 3.5.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: openblas
    num_threads: 20
         prefix: libopenblas
       filepath: /home/user/.local/lib/python3.12/site-packages/numpy.libs/libopenblas64_p-r0-0cf96a72.3.23.dev.so
        version: 0.3.23.dev
threading_layer: pthreads
   architecture: Prescott

       user_api: blas
   internal_api: openblas
    num_threads: 20
         prefix: libscipy_openblas
       filepath: /home/user/.local/lib/python3.12/site-packages/scipy.libs/libscipy_openblas-c128ec02.so
        version: 0.3.27.dev
threading_layer: pthreads
   architecture: Haswell

       user_api: openmp
   internal_api: openmp
    num_threads: 20
         prefix: libgomp
       filepath: /home/user/.local/lib/python3.12/site-packages/scikit_learn.libs/libgomp-a34b3233.so.1.0.0
        version: None

[ogrisel]

Thanks for the report @JohannesBuchner. I ran your script locally, and it does not fail on my machine.

I see that the version of the openblas library shipped with numpy has detected a "Prescott" CPU which is fishy (it's an ancient CPU architecture). Therefore, I suspect that you are running into an OpenBLAS bug such as OpenMathLib/OpenBLAS#4268.

Could you please upgrade your numpy version to see if:

• sklearn.show_versions() still report the "Prescott" architecture;
• if that solves your numerical stability problem.

If it still detects "Prescott" even with the latest version of numpy, please include the output of the lscpu command.

[ogrisel]

A unrelated issue I have is that I want to build sklearn.mixture.GaussianMixture from scratch without fitting it, but .score and .predict_proba always give 0 and 1, respectively. Is there a initialization step I am missing to be allowed to use these functions?

It wasn't really intended to be used that way. You'll have to look into the code of the fit method to know what is necessary to make the predict_proba function work as expected with your own fake "fitted" parameters.

[JohannesBuchner]

I uninstalled numpy and installed it again. I have numpy 1.26.4 (not ready for numpy-2). It now says Haswell:

$ python3 -c 'import sklearn; sklearn.show_versions()'

System:
    python: 3.12.3 (main, Sep 11 2024, 14:17:37) [GCC 13.2.0]
executable: /usr/bin/python3
   machine: Linux-6.8.0-45-generic-x86_64-with-glibc2.39

Python dependencies:
      sklearn: 1.5.0
          pip: 24.0
   setuptools: 68.1.2
        numpy: 1.26.4
        scipy: 1.14.1
       Cython: 3.0.11
       pandas: 2.2.2
   matplotlib: 3.9.2
       joblib: 1.4.2
threadpoolctl: 3.5.0

Built with OpenMP: True

threadpoolctl info:
       user_api: blas
   internal_api: openblas
    num_threads: 20
         prefix: libopenblas
       filepath: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so
        version: 0.3.26
threading_layer: pthreads
   architecture: Haswell

       user_api: blas
   internal_api: openblas
    num_threads: 20
         prefix: libscipy_openblas
       filepath: /home/user/.local/lib/python3.12/site-packages/scipy.libs/libscipy_openblas-c128ec02.so
        version: 0.3.27.dev
threading_layer: pthreads
   architecture: Haswell

       user_api: openmp
   internal_api: openmp
    num_threads: 20
         prefix: libgomp
       filepath: /home/user/.local/lib/python3.12/site-packages/scikit_learn.libs/libgomp-a34b3233.so.1.0.0
        version: None

I still get the same errors.

I was having these issues also when running on github actions.

lscpu:

⌂26.94 [user:~/Downloads/gmm-tests] main* 1 ± lscpu 
Architecture:             x86_64
  CPU op-mode(s):         32-bit, 64-bit
  Address sizes:          39 bits physical, 48 bits virtual
  Byte Order:             Little Endian
CPU(s):                   20
  On-line CPU(s) list:    0-19
Vendor ID:                GenuineIntel
  Model name:             13th Gen Intel(R) Core(TM) i7-13700H
    CPU family:           6
    Model:                186
    Thread(s) per core:   2
    Core(s) per socket:   14
    Socket(s):            1
    Stepping:             2
    CPU(s) scaling MHz:   36%
    CPU max MHz:          5000,0000
    CPU min MHz:          400,0000
    BogoMIPS:             5836,80
    Flags:                fpu vme de pse tsc msr pae mce cx8 apic sep mtrr pge mca cmov pat pse36 clflush dts acpi 
                          mmx fxsr sse sse2 ss ht tm pbe syscall nx pdpe1gb rdtscp lm constant_tsc art arch_perfmon
                           pebs bts rep_good nopl xtopology nonstop_tsc cpuid aperfmperf tsc_known_freq pni pclmulq
                          dq dtes64 monitor ds_cpl vmx smx est tm2 ssse3 sdbg fma cx16 xtpr pdcm pcid sse4_1 sse4_2
                           x2apic movbe popcnt tsc_deadline_timer aes xsave avx f16c rdrand lahf_lm abm 3dnowprefet
                          ch cpuid_fault epb ssbd ibrs ibpb stibp ibrs_enhanced tpr_shadow flexpriority ept vpid ep
                          t_ad fsgsbase tsc_adjust bmi1 avx2 smep bmi2 erms invpcid rdseed adx smap clflushopt clwb
                           intel_pt sha_ni xsaveopt xsavec xgetbv1 xsaves split_lock_detect user_shstk avx_vnni dth
                          erm ida arat pln pts hwp hwp_notify hwp_act_window hwp_epp hwp_pkg_req hfi vnmi umip pku 
                          ospke waitpkg gfni vaes vpclmulqdq rdpid movdiri movdir64b fsrm md_clear serialize arch_l
                          br ibt flush_l1d arch_capabilities
Virtualization features:  
  Virtualization:         VT-x
Caches (sum of all):      
  L1d:                    544 KiB (14 instances)
  L1i:                    704 KiB (14 instances)
  L2:                     11,5 MiB (8 instances)
  L3:                     24 MiB (1 instance)
NUMA:                     
  NUMA node(s):           1
  NUMA node0 CPU(s):      0-19
Vulnerabilities:          
  Gather data sampling:   Not affected
  Itlb multihit:          Not affected
  L1tf:                   Not affected
  Mds:                    Not affected
  Meltdown:               Not affected
  Mmio stale data:        Not affected
  Reg file data sampling: Mitigation; Clear Register File
  Retbleed:               Not affected
  Spec rstack overflow:   Not affected
  Spec store bypass:      Mitigation; Speculative Store Bypass disabled via prctl
  Spectre v1:             Mitigation; usercopy/swapgs barriers and __user pointer sanitization
  Spectre v2:             Mitigation; Enhanced / Automatic IBRS; IBPB conditional; RSB filling; PBRSB-eIBRS SW sequ
                          ence; BHI BHI_DIS_S
  Srbds:                  Not affected
  Tsx async abort:        Not affected

[JohannesBuchner]

A unrelated issue I have is that I want to build sklearn.mixture.GaussianMixture from scratch without fitting it, but .score and .predict_proba always give 0 and 1, respectively. Is there a initialization step I am missing to be allowed to use these functions?

It wasn't really intended to be used that way. You'll have to look into the code of the fit method to know what is necessary to make the predict_proba function work as expected with your own fake "fitted" parameters.

That's the confusing part, there is no return 1 if not fitted or similar in
https://github.com/scikit-learn/scikit-learn/blob/main/sklearn/mixture/_base.py#L153 and no error is raised.
In https://github.com/scikit-learn/scikit-learn/blob/main/sklearn/mixture/_gaussian_mixture.py I see no relevant code.

Somehow unpickling a model gets away with doing this though.

[ogrisel]

Please include the full traceback of the failure. It would help make the report less ambiguous.

Do you simply run the linked script unchanged, or do you do something else?

Is the reproducer fully deterministic (same values for the compared arrays at each run) or is there some random aspect involved through hypothesis?

It would help if you could post a minimal reproducer that trims out as many lines as possible in your script.

[JohannesBuchner]

Please see the full log at:
sklearn.log
The error says which arguments were given (ndiv, mu, cov, x). The code lines are already trimmed to a minimum, the relevant ones are https://github.com/JohannesBuchner/gmm-tests/blob/main/test_sklearn.py#L190-L208.

Everything above "def test_single(mean_cov):" is only for fuzzing and discovering other test cases.

With the explicitly declared test cases, (@example) the reproducer is already deterministic. I believe only if they all pass, hypothesis will try generating new ones.

[ogrisel]

Actually, I can reproduce the failure. I must have run the script without pytest the first time. 🤦

I will try to investigate tomorrow.