search

mne-tools / mne-python

MNE: Magnetoencephalography (MEG) and Electroencephalography (EEG) in Python
https://mne.tools
BSD 3-Clause "New" or "Revised" License
2.69k stars 1.31k forks source link

Regularisation of rank-deficient matrices for beamformers incorrect? #5520

Closed olafhauk closed 5 years ago

olafhauk commented 6 years ago

The function _reg_pinv() in mne.beamformer._lcmv uses Tikhonov regularisation by adding the identity to the covariance matrix. This is only correct for full-rank matrices. It is equivalent to inverting via SVD, but filtering the inverted singular values 1/s by s.^2/(s.^2 + lambda). This obviously shouldn't be applied to zero singular values, i.e. in the case of a rank-deficient matrix. Adding the identity renders the matrix full rank, and re-introduces the singular vectors associated with zero singular values into the inverse. This can potentially produce problems for example for SSSed (MEG) or average-referenced (EEG) data, even with small lambdas. Maybe this is already taken care off somewhere, but I couldn't see it.

agramfort commented 6 years ago

@olafhauk are you talking about https://github.com/mne-tools/mne-python/blob/master/mne/beamformer/_compute_beamformer.py#L35?

I see that the rank is estimated first so it should be correct no?

britta-wstnr commented 5 years ago

We discussed with @olafhauk and @vlitvak yesterday. Let's look into it at the sprint together with @vlitvak

olafhauk commented 5 years ago

We discussed with @olafhauk and @vlitvak yesterday. Let's look into it at the sprint together with @vlitvak

Looking at the code now, it seems to do the right thing (taking rank into account for Tikhonov regularisation), assuming that the rank is determined correctly.

larsoner commented 5 years ago

Okay I'll close for now but we'll keep it in mind as we try to unify minimum-norm and beamformer treatment of rank and reg