I am running into problems when trying to use LCMV beamformers on evoked MEG data (mags & grads) with a volume source space.
using an empty-room noise covariance matrix and a data covariance matrix computed from the epochs leads to confetti-like source activation values (upper right plot)
running the same pipeline with a noise covariance computed on the epochsâ baseline works fine (lower right)
Discussing this with @agramfort, he pointed out the 2 bumps in the covariance matrices computed on epochs (both baseline and post-stimulus time windows) â @larsoner what do you think about those?
Can you see if this problem occurs with minimum-norm localization of the same data? This will tell us if itâs a problem with beamforming specifically or the whitening that comes from the noise covariance. My guess is that itâs the noise cov but weâll see!
Did you process the task data and the empty-room data the same way, looks like using MaxFilter? I have some code for processing empty-room data the same way as run data, for example by hacking in the info['dev_head_t'] and setting destination the same way and actually using the head position from an actual run with movement compensation, etc. I havenât implemented this in MNE because I doubt it affects other people very much but if it turns out to be a not-identical-MaxFilter-processing issue we can figure out how to incorporate it. Also FIX: Improve covariance code · Issue #4676 · mne-tools/mne-python · GitHub might be to blame.
If you have a surface source space for this subject (recon w/FreeSurfer), can you try using a surface source space instead? I imagine the same issue will come up but if itâs easy to test we might as well.
We hacked something like this into the MNE Study Template because this is something we had to deal with for numerous datasets. Just saying this to let you know there is, in fact, a demand, and if you believe this could somehow be added to MNE directly, Iâd greatly appreciate it (and Iâd be willing to help test)
I tried different combinations and the stc looks similarly scattered with:
MNE/ dSPM on same volume source space
lcmv on surface source space
It looks normal however when using dSPM on the surface source space only, with the same noise-cov.
I used the study-template code by @richard to manually change the dev_head_t and dig of the empty room to the one from a reference run, and otherwise the same maxfilter parameters and filters as for the subjectsâ data.
The bad STCs seem to depend on the rank of the covariance matrix: when I compute the noise cov with rank=âinfoâ, it estimates a rank of 72, while the data has a rank of 73 (multiple runs combined).
When I instead compute the noise cov with rank = {meg=â73â} (estimated from the data) the resulting stc looks fine (volume source space, LCMV).
For make_lcmv I use the following rank parameters: rank = {meg=â73â} and reduce_rank=True; without explicitly setting the rank it gave me an error saying the rank of the noise and data covariance doesnât match. Could that be the issue finally?
how does this happen? when you pass rank=âinfoâ to compute_covariance and make_lcmv and the info in the same it should be 72 everywhere. Itâs not the case? Any explanation based on your code?
@agramfort I think this happens because maxfilter is applied on each of 8 raw runs separately, so they end up having slightly different ranks. Then the runs are combined to epochs with one resulting rank.
The empty room gets the dev_head_t from one reference run, so it might end up having a different rank than the epochs (all ranks range between 71-74).
For compute_covariance I use the info from the empty room raw, and for the lcmv the info from the epochs.
make_lcmv did issue a warning about the mismatch in ranks between noise and data cov, when rank=âinfoâ, which I thought I had overcome by explicitly providing the rank from the data (no warning).
There is no message about the mismatch if the rank is set explicitly.
I still donât fully grasp why this mismatch of ranks by 1 causes such trouble
since the rank of the epochs does anyways not reflect the rank of every run.
I would use the rank of the matrix that is doing the whitening. So @sophie in your case compute rank from the empty-room raw, not from the epoched / evoked task data, or process the empty-room in an identical way and you will probably be safe enough â but even then itâs probably safest to compute rank from the empty-room data.
Letâs say you have two matrices, noise = np.diag([1, 1, 1e-14, 0.9e-14]) and data=np.diag([1, 1, 1, 1e-14]). If we want to whiten data using the matrix noise, the correct behavior is to use rank=2 and your resulting data matrix will have one component thrown away because it is truly of rank 3. You lose a bit of information this way, but itâs probably the best you can do with a mismatch of spatial subspaces. And at least with MaxFiltered data, the mismatch will occur in the highest (spatial) orders, so youâll be losing some of your lowest-SNR components anyway (hopefully no huge loss).
OTOH if you whiten with rank=3 your whitened data matrix will be totally broken, because your whitening matrix will be something like np.diag([1, 1, 1e14, 0]), and 1e14 will dominate even though this was a valid component in your data. And worse, sometimes the subspace mismatch will mean that this component wasnât even present in your data so you can be literally amplifying noise rather than signal.
Processing your empty room the same way as generally raw helps with these situations because it helps ensure (but does not guarantee, at least in the case of movement compensation) that the high-order components that are thrown away during regularization are the same in both cases. And something like the proposals here would allow us to explicitly keep track of which spatial components are present and maybe do a better job of reconciling them:
Concretely though I would say that even though in general rank problems can be tricky, if we special-case MEG channels (and maybe even MEG+MaxFilter?) without reference channels we can probably detect this sort of problem. @sophie Can you share your noise and data covariances somewhere?
Hi @larsoner, thank you for the elaborate response!
Does your recommendation to use the rank of the noise covariance matrix also hold when it has a higher rank than the data covariance matrix (by 1 or 2)?
Yes it should still apply in that case. Think about the reverse case where you have data = np.diag([1, 1, 1e-14, 0.9e-14]) and noise = np.diag([1, 1, 1, 1e-14]). Your whitener with rank=3 based on noise will be np.diag([1, 1, 1, 0] and your whitened data will be np.diag([1, 1, 1e-14, 0], i.e., the component that was essentially zero in your original data will stay that way.
Hello @sophie, whatâs the current status of this? Did you manage to resolve the issue? I would be very curious to learn if and how you got it resolved, so we can follow the same approach in the Study Template. Thanks!
Hi @richard, I did follow Ericâs suggestion, to use the rank of the matrix used for whitening and the stcâs look ok now. I think this point would merit to be made explicit somewhere as a recommendation.
I will clean up my code and we can discuss for the study template.
Not sure what the âcorrectâ approach would be!