mne.compute_covariance raises error with method="auto"

[barbaracassone]

Description of the problem

I am trying to follow the steps in this tutorial but the "compute_covariance" function results in an algorithmic convergence failure.

Steps to reproduce

I am running the Jupyter notebook "evoked_whitening.ipynb" downloadable at this link (https://mne.tools/stable/auto_examples/visualization/evoked_whitening.html) in a conda environment:

conda v26.1.1
jupyter v1.1.1
mne-base v1.11.0
python v3.13.9

The following lines raise an error:

method_params = dict(diagonal_fixed=dict(mag=0.01, grad=0.01, eeg=0.01))
noise_covs = compute_covariance(
    epochs,
    tmin=None,
    tmax=0,
    method="auto",
    return_estimators=True,
    n_jobs=None,
    projs=None,
    rank=None,
    method_params=method_params,
    verbose=True,
)

# With "return_estimator=True" all estimated covariances sorted
# by log-likelihood are returned.

print("Covariance estimates sorted from best to worst")
for c in noise_covs:
    print(f"{c['method']} : {c['loglik']}")

Link to data

No response

Expected results

As in this tutorial.

Actual results

Created an SSP operator (subspace dimension = 4)
    Setting small MEG eigenvalues to zero (without PCA)
    Setting small EEG eigenvalues to zero (without PCA)
Reducing data rank from 364 -> 360
Estimating covariance using SHRUNK
C:\Users\Utente\anaconda3\Lib\site-packages\sklearn\model_selection\_validation.py:953: UserWarning: Scoring failed. The score on this train-test partition for these parameters will be set to nan. Details: 
Traceback (most recent call last):
  File "C:\Users\Utente\anaconda3\Lib\site-packages\sklearn\model_selection\_validation.py", line 940, in _score
    scores = scorer(estimator, X_test, **score_params)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\sklearn\metrics\_scorer.py", line 492, in __call__
    return estimator.score(*args, **kwargs)
           ~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^
  File "C:\Users\Utente\anaconda3\Lib\site-packages\sklearn\covariance\_empirical_covariance.py", line 288, in score
    res = log_likelihood(test_cov, self.get_precision())
                                   ~~~~~~~~~~~~~~~~~~^^
  File "C:\Users\Utente\anaconda3\Lib\site-packages\sklearn\covariance\_empirical_covariance.py", line 230, in get_precision
    precision = linalg.pinvh(self.covariance_, check_finite=False)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\_lib\_util.py", line 1233, in wrapper
    return f(*arrays, *other_args, **kwargs)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\linalg\_basic.py", line 1726, in pinvh
    s, u = _decomp.eigh(a, lower=lower, check_finite=False, driver='ev')
           ~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\_lib\_util.py", line 1233, in wrapper
    return f(*arrays, *other_args, **kwargs)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\linalg\_decomp.py", line 620, in eigh
    raise LinAlgError(msg)
numpy.linalg.LinAlgError: The algorithm failed to converge; 202 off-diagonal elements of an intermediate tridiagonal form did not converge to zero.

  warnings.warn(
C:\Users\Utente\anaconda3\Lib\site-packages\sklearn\model_selection\_validation.py:953: UserWarning: Scoring failed. The score on this train-test partition for these parameters will be set to nan. Details: 
Traceback (most recent call last):
  File "C:\Users\Utente\anaconda3\Lib\site-packages\sklearn\model_selection\_validation.py", line 940, in _score
    scores = scorer(estimator, X_test, **score_params)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\sklearn\metrics\_scorer.py", line 492, in __call__
    return estimator.score(*args, **kwargs)
           ~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^
  File "C:\Users\Utente\anaconda3\Lib\site-packages\sklearn\covariance\_empirical_covariance.py", line 288, in score
    res = log_likelihood(test_cov, self.get_precision())
                                   ~~~~~~~~~~~~~~~~~~^^
  File "C:\Users\Utente\anaconda3\Lib\site-packages\sklearn\covariance\_empirical_covariance.py", line 230, in get_precision
    precision = linalg.pinvh(self.covariance_, check_finite=False)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\_lib\_util.py", line 1233, in wrapper
    return f(*arrays, *other_args, **kwargs)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\linalg\_basic.py", line 1726, in pinvh
    s, u = _decomp.eigh(a, lower=lower, check_finite=False, driver='ev')
           ~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\_lib\_util.py", line 1233, in wrapper
    return f(*arrays, *other_args, **kwargs)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\linalg\_decomp.py", line 620, in eigh
    raise LinAlgError(msg)
numpy.linalg.LinAlgError: The algorithm failed to converge; 202 off-diagonal elements of an intermediate tridiagonal form did not converge to zero.

 [...] 

    warnings.warn(
C:\Users\Utente\anaconda3\Lib\site-packages\sklearn\model_selection\_search.py:1135: UserWarning: One or more of the test scores are non-finite: [           nan            nan            nan            nan
            nan            nan            nan            nan
            nan            nan           -inf            nan
            nan            nan            nan           -inf
            nan           -inf            nan            nan
 -1032.04034359 -1039.33041563 -1048.84139831 -1060.94429543
            nan            nan            nan            nan
 -1158.92270153 -1246.99239465]
  warnings.warn(
C:\Users\Utente\anaconda3\Lib\site-packages\sklearn\model_selection\_search.py:1146: RuntimeWarning: invalid value encountered in subtract
  (array - array_means[:, np.newaxis]) ** 2, axis=1, weights=weights
Done.
Estimating covariance using DIAGONAL_FIXED
    MAG regularization : 0.01
    GRAD regularization : 0.01
    EEG regularization : 0.01
Done.
Estimating covariance using EMPIRICAL
Done.
Using cross-validation to select the best estimator.
    MAG regularization : 0.01
    GRAD regularization : 0.01
    EEG regularization : 0.01
C:\Users\Utente\anaconda3\Lib\site-packages\sklearn\model_selection\_validation.py:971: UserWarning: Scoring failed. The score on this train-test partition for these parameters will be set to nan. Details: 
Traceback (most recent call last):
  File "C:\Users\Utente\anaconda3\Lib\site-packages\sklearn\metrics\_scorer.py", line 156, in __call__
    score = scorer(estimator, *args, **routed_params.get(name).score)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\mne\cov.py", line 1443, in _gaussian_loglik_scorer
    precision = est.get_precision()
  File "C:\Users\Utente\anaconda3\Lib\site-packages\mne\cov.py", line 1681, in get_precision
    return self.estimator_.get_precision()
           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^
  File "C:\Users\Utente\anaconda3\Lib\site-packages\mne\fixes.py", line 301, in get_precision
    precision = linalg.pinvh(self.covariance_)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\_lib\_util.py", line 1233, in wrapper
    return f(*arrays, *other_args, **kwargs)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\linalg\_basic.py", line 1726, in pinvh
    s, u = _decomp.eigh(a, lower=lower, check_finite=False, driver='ev')
           ~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\_lib\_util.py", line 1233, in wrapper
    return f(*arrays, *other_args, **kwargs)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\linalg\_decomp.py", line 620, in eigh
    raise LinAlgError(msg)
numpy.linalg.LinAlgError: The algorithm failed to converge; 359 off-diagonal elements of an intermediate tridiagonal form did not converge to zero.

  warnings.warn(
    MAG regularization : 0.01
    GRAD regularization : 0.01
    EEG regularization : 0.01
C:\Users\Utente\anaconda3\Lib\site-packages\sklearn\model_selection\_validation.py:971: UserWarning: Scoring failed. The score on this train-test partition for these parameters will be set to nan. Details: 
Traceback (most recent call last):
  File "C:\Users\Utente\anaconda3\Lib\site-packages\sklearn\metrics\_scorer.py", line 156, in __call__
    score = scorer(estimator, *args, **routed_params.get(name).score)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\mne\cov.py", line 1443, in _gaussian_loglik_scorer
    precision = est.get_precision()
  File "C:\Users\Utente\anaconda3\Lib\site-packages\mne\cov.py", line 1612, in get_precision
    return self.estimator_.get_precision()
           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^
  File "C:\Users\Utente\anaconda3\Lib\site-packages\mne\fixes.py", line 301, in get_precision
    precision = linalg.pinvh(self.covariance_)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\_lib\_util.py", line 1233, in wrapper
    return f(*arrays, *other_args, **kwargs)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\linalg\_basic.py", line 1726, in pinvh
    s, u = _decomp.eigh(a, lower=lower, check_finite=False, driver='ev')
           ~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\_lib\_util.py", line 1233, in wrapper
    return f(*arrays, *other_args, **kwargs)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\linalg\_decomp.py", line 620, in eigh
    raise LinAlgError(msg)
numpy.linalg.LinAlgError: The algorithm failed to converge; 358 off-diagonal elements of an intermediate tridiagonal form did not converge to zero.

[...] 
  
  warnings.warn(
Number of samples used : 1891
log-likelihood on unseen data (descending order):
   shrunk: nan
   diagonal_fixed: nan
   empirical: nan
[done]
Covariance estimates sorted from best to worst
shrunk : nan
diagonal_fixed : nan
empirical : nan
C:\Users\Utente\anaconda3\Lib\site-packages\sklearn\model_selection\_validation.py:971: UserWarning: Scoring failed. The score on this train-test partition for these parameters will be set to nan. Details: 
Traceback (most recent call last):
  File "C:\Users\Utente\anaconda3\Lib\site-packages\sklearn\metrics\_scorer.py", line 156, in __call__
    score = scorer(estimator, *args, **routed_params.get(name).score)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\mne\cov.py", line 1443, in _gaussian_loglik_scorer
    precision = est.get_precision()
  File "C:\Users\Utente\anaconda3\Lib\site-packages\mne\fixes.py", line 301, in get_precision
    precision = linalg.pinvh(self.covariance_)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\_lib\_util.py", line 1233, in wrapper
    return f(*arrays, *other_args, **kwargs)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\linalg\_basic.py", line 1726, in pinvh
    s, u = _decomp.eigh(a, lower=lower, check_finite=False, driver='ev')
           ~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\_lib\_util.py", line 1233, in wrapper
    return f(*arrays, *other_args, **kwargs)
  File "C:\Users\Utente\anaconda3\Lib\site-packages\scipy\linalg\_decomp.py", line 620, in eigh
    raise LinAlgError(msg)
numpy.linalg.LinAlgError: The algorithm failed to converge; 358 off-diagonal elements of an intermediate tridiagonal form did not converge to zero.

  warnings.warn(

Additional information

Platform             Windows-11-10.0.26200-SP0
Python               3.13.9 | packaged by Anaconda, Inc. | (main, Oct 21 2025, 19:09:58) [MSC v.1929 64 bit (AMD64)]
Executable           C:\Users\Utente\anaconda3\python.exe
CPU                  AMD Ryzen 5 5600H with Radeon Graphics (12 cores)
Memory               15.3 GiB

Core
 + mne               1.11.0 (latest release)
 + numpy             2.3.5 (MKL 2025.0-Product with 6 threads)
 + scipy             1.16.3
 + matplotlib        3.10.6 (backend=module://matplotlib_inline.backend_inline)

Numerical (optional)
 + sklearn           1.7.2
 + numba             0.62.1
 + nibabel           5.4.0
 + pandas            2.3.3
 + h5io              0.2.5
 + h5py              3.15.1
 - unavailable       nilearn, dipy, openmeeg, cupy

Visualization (optional)
 + pyvista           0.47.1 (OpenGL 4.6.0 Compatibility Profile Context 22.20.44.221025 via AMD Radeon(TM) Graphics)
 + pyvistaqt         0.11.3
 + vtk               9.5.2
 + qtpy              2.4.3 (PyQt5=5.15.2)
 + ipympl            0.9.7
 + ipywidgets        8.1.7
 - unavailable       pyqtgraph, mne-qt-browser, trame_client, trame_server, trame_vtk, trame_vuetify

Ecosystem (optional)
 + defusedxml        0.7.1
 - unavailable       mne-bids, mne-nirs, mne-features, mne-connectivity, mne-icalabel, mne-bids-pipeline, neo, eeglabio, edfio, curryreader, mffpy, pybv, antio

[tsbinns]

Hi @barbaracassone, thanks for opening this issue.

A very basic question, but just to confirm: you are running the example exactly as given in the file? You aren't changing the data used or any of the preprocessing steps?

The following clean install runs the example without issues for me:

Platform             Linux-6.17.0-14-generic-x86_64-with-glibc2.39
Python               3.13.11 (main, Jan 14 2026, 19:38:04) [Clang 21.1.4 ]
Executable           /home/tsbinns/github/mne-python/.venv/bin/python
CPU                  AMD Ryzen 5 1600 Six-Core Processor (12 cores)
Memory               15.5 GiB

Core
├☑ mne               1.11.0 (latest release)
├☑ numpy             2.4.4 (OpenBLAS 0.3.31.188.0 with 12 threads)
├☑ scipy             1.17.1
└☑ matplotlib        3.10.8 (backend=tkagg)

Numerical (optional)
├☑ sklearn           1.8.0
└☐ unavailable       numba, nibabel, nilearn, dipy, openmeeg, cupy, pandas, h5io, h5py

Visualization (optional)
└☐ unavailable       pyvista, pyvistaqt, vtk, qtpy, ipympl, pyqtgraph, mne-qt-browser, ipywidgets, trame_client, trame_server, trame_vtk, trame_vuetify

Ecosystem (optional)
└☐ unavailable       mne-bids, mne-nirs, mne-features, mne-connectivity, mne-icalabel, mne-bids-pipeline, neo, eeglabio, edfio, curryreader, mffpy, pybv, pymef, antio, defusedxml

Do you still see this problem if you upgrade your scipy and sklearn versions?

[barbaracassone]

Hello @tsbinns, thank you for your reply.

I confirm that I am running the code exactly as it is provided, without making any changes.
I upgraded all the libraries you mentioned, but nothing changed.
I am no expert, so I don't know how or why this is happening, but I suspect it may be a Windows-related issue, because the same script runs without any problems on my professor's Linux machine.

[tsbinns]

Thanks for confirming @barbaracassone! I'll try and replicate this on my Windows machine and work from there.

[larsoner]

 + numpy             2.3.5 (MKL 2025.0-Product with 6 threads)

I suspect this is perhaps some bug with MKL. Perhaps conda install nomkl (which will force a switch to OpenBLAS) would fix it...

[JiwaniZakir]

The error originates in sklearn's cross-validation loop used internally by the SHRUNK estimator (via sklearn.covariance.ShrunkCovariance) to select the optimal shrinkage parameter — the "Scoring failed" warning means the log-likelihood scorer is throwing an exception on at least one train/test fold, likely because the fold's covariance estimate is singular or near-singular. The rank reduction from 364 → 360 due to SSP projectors is the probable trigger: with a rank-deficient data matrix, some CV folds end up with a covariance that cannot be inverted cleanly for log-likelihood scoring, causing sklearn to set that fold's score to nan and ultimately raising downstream. This behavior is likely exacerbated by a stricter numerical handling or API change in newer sklearn versions compatible with Python 3.13 — worth checking the sklearn version in the environment and comparing against the _compute_covariance_auto path in mne/cov.py to see whether the rank/regularization is being passed through correctly to the scorer.

[tsbinns]

I suspect this is perhaps some bug with MKL. Perhaps conda install nomkl (which will force a switch to OpenBLAS) would fix it...

Thanks for the suggestion @larsoner. Just tried this on my Windows machine with MKL 2025.0 and I do see the error, but I also get a warning on install that MKL 2025.0 is not compatible with numpy 2.3.5.

@barbaracassone Upgrading to MKL 2025.3 allows me to run the example without issues, or you can switch to OpenBLAS.

[barbaracassone]

Upgrading to MKL 2025.3 solved the issue. Thank you all for your time and assistance!