ValueError: Covariance matrices must be positive definite. Add regularization to avoid this error.

[FarnoodF]

I run into this error for some of my EEG inputs when I run the meegkit.asr.asr_process():

ValueError: Covariance matrices must be positive definite. Add regularization to avoid this error.

Complete error message:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-199-a17002d5a7b8> in <module>
     52             am = 0
     53             for i in range(len(X_t_w)):
---> 54                 x_t_w = asrC.transform(X_t_w[i].T).T
     55                 # x_t_w = octave.asr_process(X_t_w[i].T, sfreq, asr_state).T
     56                 X_t = avaDM.LagGenerator(x_t_w, lags)

c:\~\python\python38\lib\site-packages\meegkit\asr.py in transform(self, X, y, **kwargs)
    242 
    243         # Clean data, using covariances weighted by sample_weight
--> 244         out, self.state_ = asr_process(X, X_filt, self.state_,
    245                                        cov=np.stack(self.cov_),
    246                                        method=self.method,

c:\~\python\python38\lib\site-packages\meegkit\asr.py in asr_process(X, X_filt, state, cov, detrend, method, sample_weight)
    577     else:
    578         if cov.ndim == 3:
--> 579             cov = pyriemann.utils.mean.mean_covariance(
    580                 cov, metric='riemann', sample_weight=sample_weight)
    581 

c:\~\python\python38\lib\site-packages\pyriemann\utils\mean.py in mean_covariance(covmats, metric, sample_weight, *args)
    330         C = metric(covmats, sample_weight=sample_weight, *args)
    331     else:
--> 332         C = mean_methods[metric](covmats, sample_weight=sample_weight, *args)
    333     return C
    334 

c:\~\python\python38\lib\site-packages\pyriemann\utils\mean.py in mean_riemann(covmats, tol, maxiter, init, sample_weight)
     59         for index in range(Nt):
     60             tmp = numpy.dot(numpy.dot(Cm12, covmats[index, :, :]), Cm12)
---> 61             J += sample_weight[index] * logm(tmp)
     62 
     63         crit = numpy.linalg.norm(J, ord='fro')

c:\~\python\python38\lib\site-packages\pyriemann\utils\base.py in logm(Ci)
     44 
     45     """
---> 46     return _matrix_operator(Ci, numpy.log)
     47 
     48 

c:\~\python\python38\lib\site-packages\pyriemann\utils\base.py in _matrix_operator(Ci, operator)
      8     """matrix equivalent of an operator."""
      9     if Ci.dtype.char in typecodes['AllFloat'] and not numpy.isfinite(Ci).all():
---> 10         raise ValueError("Covariance matrices must be positive definite. Add regularization to avoid this error.")
     11     eigvals, eigvects = scipy.linalg.eigh(Ci, check_finite=False)
     12     eigvals = numpy.diag(operator(eigvals))

ValueError: Covariance matrices must be positive definite. Add regularization to avoid this error.

Running the same code on Matlab using the same EEG does not give me the error.

[nbara]

The traceback is quite explicit here. There's something wrong with the data that causes the covariance matrices to not be SPD.

This can happen if:

• you have more channels than time samples;
• for some reason your data is rank-deficient (if you applied common average reference or some component analysis beforehand);
• maybe you have included a portion of your data which is grossly artefactual (eg entirely flat or NaN?)

If you can not fix the root cause of the problem, then we will need to add regularization inside ASR.transform(). I'll try to provide an option to do that shortly.

Note that in general, Matlab is a bit more flexible that Python and handles these edge cases more gracefully. In Python we have to take care of each scenario explicitly...

[FarnoodF]

The EEG data looks normal: no flat or nan areas, and the number of samples is significantly more than number of channels. I also have not applied any component analysis beforehand. I think the issue arises from the CAR (Common Average Reference) as you mentioned. But CAR is a very common practice in BCI systems to lower the common noise in multiple channels.

At this point I don't think I can do anything on the input side. It would be better if we had the regularization option.

I also looked at your asr_process code and it is implemented differently from its EEGLAB counterpart, is there a reason for that?

[nbara]

At this point I don't think I can do anything on the input side. It would be better if we had the regularization option.

I will try to add an option to perform regularization.

I also looked at your asr_process code and it is implemented differently from its EEGLAB counterpart, is there a reason for that?

The original EEGLAB code is not the cleanest I've ever seen (that's an understatement).

And there are many utilities in Matlab that are not available in numpy or scipy natively so I had to come up with alternatives.

I've done my best to make a working implementation, but I've also taken the liberty to A) clean things up when I deemed it necessary B) deviate from the Matlab code when it made no sense to me.

[nbara]

Can you try https://github.com/nbara/python-meegkit/pull/18 to see if it fixes your problem? You can now specify a covariance estimator with regularization.

[FarnoodF]

@nbara Thanks, it works! But the ASR.transform() output is complex now with very small imaginary values. Probably add np.real()?

[nbara]

Can you check #18 again? This should be solved and I'm now properly testing that the output is real-valued.