CSP error: LinAlgError: The leading minor of order 64 of B is not positive definite

Hi,

I face the following error when running the CSP function using the EEG BCI dataset:

LinAlgError : The leading minor of order 64 of B is not positive definite. The factorization of B could not be completed and no eigenvalues or eigenvectors were computed.

There seems to be no structure in when this error occurs: sometimes I receive an error, sometimes not (with the same data). I tested the code on subjects with indeces 0, 1, and 2 on the EEGBCI dataset on two PCs (with sometimes the same result, sometimes different). In response to a previous Forum post (‘CSP error’, 2017), we changed the reg parameter from 1e-10 to oas, but it still gives us the error sometimes. We provide the code below. Do you have a suggestion on how to solve this issue?

Best regards,
Friederike

  • MNE-Python version: 0.22.0
  • operating system: Windows 10

Code:

from mne.channels import make_standard_montage
from mne.io import concatenate_raws, read_raw_edf
from mne.datasets import eegbci
from mne.preprocessing import ICA, corrmap
import numpy as np
from mne import Epochs, pick_types, events_from_annotations
from sklearn.model_selection import RepeatedStratifiedKFold
from mne.decoding import CSP

# define load_data function:
def load_data(subjects, runs=[6, 10, 14]):
    raws = [None] * len(subjects)

    for i, subject in enumerate(subjects):
        raw_fnames = eegbci.load_data(subject+1, runs)
        raw = concatenate_raws([read_raw_edf(f, preload=True) for f in raw_fnames])
        eegbci.standardize(raw) 
        standard_1005 = make_standard_montage('standard_1005')
        raw.set_montage(standard_1005)
        raw.rename_channels(lambda x: x.strip('.'))
        raws[i] = raw

    return raws

# define ICA_analysis function
def ICA_analysis(raws):
    raw_template = load_data(subjects=[0])
    raw_template[0].filter(7, 30, fir_design='firwin', skip_by_annotation='edge')
    ica_template = ICA(n_components=20, random_state=0, max_iter=800).fit(raw_template[0])
    template_eog_component = ica_template.get_components()[:, 7]
    template_ecg_component = ica_template.get_components()[:, 0]
    icas = [None] * len(raws)

    for i, raw in enumerate(raws):
        icas[i] = ICA(n_components=20, random_state=0, max_iter=800).fit(raw)

    try:
        corrmap(icas, template=template_eog_component, threshold=0.8, label='blink', plot=False)
    except:
        None

    try:
        corrmap(icas, template=template_ecg_component, threshold=0.8, label='QRS', plot=False)
    except:
        None

    for subject in range(len(icas)):
        remove_list = np.empty([0], dtype='int')

        for k in icas[subject].labels_.values():
            remove_list = np.append(remove_list, k)

        print("remove_list: {}".format(remove_list))
        icas[subject].exclude = remove_list
        raws[subject] = icas[subject].apply(raws[subject])

    return raws

# define preprocessing function:
def preprocessing(raws):

    # bandpass filter raw object
    for subject in range(len(raws)):
        raws[subject].filter(7, 30, fir_design='firwin', skip_by_annotation='edge')

    # # ICA filtering raw object
    raws = ICA_analysis(raws)

    # Repairing artifacts with average reference as a projector
    for subject in range(len(raws)):
        raws[subject].set_eeg_reference('average', projection=True).apply_proj()

    # Crop and epoch data 
    epochs = [None] * len(raws)
    tmin, tmax = 0, 521/160
    n_samples = round(raws[0].info['sfreq'] * (tmax - tmin) + 1)
    n_channels = len(raws[0].info['chs'])
    X = np.empty([0,n_channels,n_samples])
    y = np.empty([0], dtype='int')

    for subject, raw in enumerate(raws):
        events, _ = events_from_annotations(raw, event_id=dict(T1=2, T2=3))
        picks = pick_types(raw.info, eeg=True, exclude='bads')
        event_id = dict(hands=2, feet=3)
        epochs[subject] = Epochs(raw, events, event_id, tmin, tmax, proj=True, picks=picks, baseline=No 
        ne, preload=True)

        X = np.append(X, epochs[subject].get_data(),0) 
        y = np.append(y, epochs[subject].events[:, -1] - 2)
    return X, y    


raws_train = load_data(subjects=[0]) # we tried [0,1,2]
X_train, y_train = preprocessing(raws_train)

cv_inner = RepeatedStratifiedKFold(n_splits=3, n_repeats=1, random_state=0)
for train_index, test_index in cv_inner.split(X_train, y_train):
    X1, X2 = X_train[train_index], X_train[test_index]
    y1, y2 = y_train[train_index], y_train[test_index]
    a = CSP(n_components=5, reg='oas').fit(X1, y1) # reg = 1e-10 gives same error

Hello,

could you please post the output of the command

mne sys_info

?

I’ve also seen you’ve run into a similar problem before:

Did you ever manage to follow up and resolve that?

Thanks for your reply. mne.sys_info() gives the following info:

mne.sys_info()
Platform:      Windows-10-10.0.19041-SP0
Python:        3.8.3 (default, Jul  2 2020, 17:30:36) [MSC v.1916 64 bit (AMD64)]
Executable:    C:\Users\fbehrens\Anaconda3\python.exe
CPU:           Intel64 Family 6 Model 142 Stepping 10, GenuineIntel: 8 cores
Memory:        7.9 GB

mne:           0.22.0
numpy:         1.18.5 {blas=mkl_rt, lapack=mkl_rt}
scipy:         1.5.4
matplotlib:    3.3.3 {backend=module://ipykernel.pylab.backend_inline}

sklearn:       0.23.2
numba:         0.50.1
nibabel:       Not found
nilearn:       Not found
dipy:          Not found
cupy:          Not found
pandas:        1.0.5
mayavi:        Not found
pyvista:       Not found
vtk:           Not found

And yes, the issue posted before about the SVD converge error has been resolved. Unfortunately, I do not know how… I did not change anything, but now a few weeks later the plot function seems to work fine. Interestingly, the same holds for my colleague.

Do you have any idea what could be the issue?

Hello, thanks for posting the output of mne.sys_info(). I’m not an expert with this stuff, but I hope @larsoner might have an idea here! Fingers crossed :crossed_fingers:

@larsoner do you have an idea about what might cause this issue?

There was a bug with Windows that was fixed:

https://developercommunity.visualstudio.com/t/fmod-after-an-update-to-windows-2004-is-causing-a/1207405

I think the update rolled out in January, so if things are magically fixed then this might be why.

Are things still broken, or are you just curious why it’s working now?

Thanks @larsoner for the update. This might be indeed the reason why our previous issue Richard referred to above is resolved. However, we still run into the issue that is outlined in this issue where the CSP function somtimes raises LinAlgError: The leading minor of order 64 of B is not positive definite. Any thought on that?

If you set the regularization very high, say to 1e-4, then do you still get the problem?

When you hit the error, if you do post-mortem debugging (import pdb; pdb.pm()), what are the first and last singular values of the matrix being decomposed like (print(np.linalg.svd(x)[1][[0, -1]])?

Changing the regulation parameter does not help unfortunately. Regarding the debugging, import pdb; pdb.pm() kept running for more than an hour without any result. Am I doing something wrong here?

Once you run that, it should drop you to an interactive debugging session with a prompt like:

(pdb) 

And you should be able to do debugging commands like “print” (“p”):

(pdb) p np.linalg.svd(x)[1][[0, -1]]

But if you’re using IPython or something you might want to use %debug magic to get there instead

ah thanks a lot for the explanation. I got the following values for the print statement: [5.14015052e-09 1.12520815e-24]

This seems like it must be a bug with CSP, since that looks un-regularized. Can you np.savez_compressed one set of X1, X2, y1, y2 that cause the error, upload them somewhere and open a GitHub issue with a minimal script that shows the problem (probably just needs to np.load the data and call mne.decoding.CSP on it)?

1 Like

ok, I just opened a new bug issue here. Thanks a lot for the support.

Hi, I am facing the same problem. Wondering if you have solved it already? Cheers

Hey, I’m also facing the problem! Any Updates on this?

Hello @dk76 and @chicanino, unfortunately I don’t think there have been any new developments here. @larsoner, have you had any luck figuring out what’s going on?

It looks like @agramfort was unable to reproduce, so it’s probably a Windows+MKL issue. I’ll have to try it on my Windows machine next week. In the meantime I’d recommend making sure you have the latest versions of NumPy and SciPy (maybe make a new/clean env) to see if the problem persists.

Thanks for updating @richard. It seems that my problem is due to the “set_eeg_reference” function. I removed the average re-reference step and the error was gone. Don’t know if that could give someone a clue.

when you apply the average reference you reduce the rank of the data. CSP requires
you to pass full rank data to have positive definite matrices. You need to regularize the
covariance by adding a constant to the diagonal or apply a PCA before CSP

this can be done using
https://mne.tools/stable/generated/mne.decoding.UnsupervisedSpatialFilter.html#mne.decoding.UnsupervisedSpatialFilter

in a pipeline

1 Like

Hi, just thought I would drop a link to my post from last week. I ran into the same issue (I was able to avoid the mistake, however, by performing PCA on my data first and removing some noisy channels). Maybe my example can also help in pinning down the origin of the error :slight_smile:

1 Like