[WIP] Add CSD whitening to DICS beamformer #7021

wmvanvliet · 2019-11-05T08:31:33Z

Working on adding whitening support to the DICS beamformer, so we can mix different channel types.

Todo:

Implement whitening of CSD and Leadfield
Modify DICS example to use both grads and mags
Do proper check to see if data and noise CSDs are compatible (working on part of this in [WIP] Make equalize_channels capable of operating on Forward, Covariance and CrossSpectralDensity #7029)
Add new unit tests for DICS whitening
Have only one whitening matrix for all frequencies
Update What's New

larsoner

Looks like a great start, love to see a small diff bring a nice enhancement. (I know it's not done yet but it looks close!)

larsoner · 2019-11-05T15:07:25Z

examples/inverse/plot_dics_source_power.py


-# Set picks, use a single sensor type
-picks = mne.pick_types(raw.info, meg='grad', exclude='bads')
-


mne/beamformer/_dics.py

mne/time_frequency/csd.py

larsoner · 2019-11-05T15:10:05Z

mne/time_frequency/csd.py

+            from ..cov import Covariance  # to avoid circular import
+            return Covariance(data, self.ch_names, bads=[], projs=self.projs,
+                              nfree=self.n_fft)
+            return None


debugging cruft?

probably :)

britta-wstnr

Hooray for combining sensors for DICS!
I fear I came to question leadfield whitening and its effect on DICS with multiple frequencies - see my reviewing comments ... What is your opinions on this @larsoner @wmvanvliet ?

britta-wstnr · 2019-11-06T12:21:40Z

mne/beamformer/_dics.py

        computing the inverse.
+    noise_csd : instance of CrossSpectralDensity | None
+        Noise cross-spectral density (CSD) matrices. If provided, whitening
+        will be done. The noise CSDs need to have been computed for the same


The noise CSD matrices need to be computed

My reasoning with the grammar was that as you are supplying the CSDs are a parameter to this function, they already have been computed.

britta-wstnr · 2019-11-06T12:24:11Z

mne/beamformer/_dics.py

+
+        _, _, proj, vertices, G, _, nn, orient_std = _prepare_beamformer_input(
+            info, forward, label, pick_ori, noise_cov=noise_csd_freq,
+            combine_xyz=combine_xyz, exp=exp)


Does this really need to be done per frequency now? It was outside the loop before if I see it right.

Oh wait ... we now whiten each lead field differently, is that it? Hummmmm, is that something we want (I do understand why we have it)? Different lead fields per frequency ... ? It strikes me as an odd thing to have. One step further: do we actually have to whiten the lead field or would it be sufficient to whiten the data covariance /CSD matrix.

We absolutely have to whiten per frequency. We have a different CSD at each frequency, and thus a different beamformer at each frequency. It's just a convenience feature that we pack multiple CSDs into one object and multiple beamformers in a single datastructure.

However, there is now a lot of redundant work happening within _prepare_beamformer_input. I plan to refactor that function to unpack the different things it does.

Well I do understand we get different weights at each frequency. The question is whether we want to have different lead fields, too. Probably not. And that raises the general question why we whiten lead fields - we do not invert the lead field in beamforming as compared to MNE solutions.

I don't know what happens if the lead field is not whitened - if it mathematically would still pan out. Maybe worth a try! AFAIK there is no literature on lead field whitening with beamformers, only covariance matrices - but I might be missing that paper ;)

And that raises the general question why we whiten lead fields - we do not invert the lead field in beamforming as compared to MNE solutions.

At least for MNE any linear operator you apply to the data you should apply to the lead fields because otherwise the results will be biased. I assume the same would hold for beamformers, no?

For example if you apply a projector to the data, you apply it to the lead fields. If you know that this is the correct thing mathematically to do for beamformers when you apply a projection operation to the data, then it is probably mathematically also the correct thing to do in the case of whitening as well. (The whitener implicitly also applies the projection matrix FWIW.)

To me it's not so much about having "different lead fields" as it is about applying the same operators to data and the lead fields before doing computations on them. But maybe I'm thinking about it the wrong way...

Well, it is def. the case that whatever is done on the original data needs to be done on the lead field as well. However - does the data need to be whitened to take care of scaling or is it enough to do the covariance matrix?
I quickly discussed things with @sarangnemo and he shares my intuition about this maybe not being ideal. I will try to run some data through differently whitened beamformers and report back!

in the DICS code, there is no data, there is only the CSD matrix.

Here is a comparison between whitening both leadfield and the CSD (left) and only whitening the CSD (right):

Co-Authored-By: Eric Larson <larson.eric.d@gmail.com>

codecov · 2019-11-06T13:50:14Z

Codecov Report

❗ No coverage uploaded for pull request base (master@1ffebaf). Click here to learn what that means.
The diff coverage is n/a.

@@            Coverage Diff            @@
##             master    #7021   +/-   ##
=========================================
  Coverage          ?   89.76%           
=========================================
  Files             ?      444           
  Lines             ?    79430           
  Branches          ?    12699           
=========================================
  Hits              ?    71297           
  Misses            ?     5349           
  Partials          ?     2784

Co-Authored-By: Eric Larson <larson.eric.d@gmail.com>

wmvanvliet · 2019-11-07T14:59:41Z

Oops, completely forgot that _prepare_beamformer does not whiten the CSD. So we've been looking at results with only a whitened leadfield. Fixed now. Results look GOOD! :)

larsoner · 2019-11-07T17:29:18Z

Azure style is real:

https://dev.azure.com/mne-tools/mne-python/_build/results?buildId=4348&view=logs&jobId=fabb5f62-f417-5e9e-e129-2b63f80d6540&j=fabb5f62-f417-5e9e-e129-2b63f80d6540&t=dbab9c63-5656-5fa8-4c44-6c5548889c04

Looks like all beamformer tests fail:

https://travis-ci.org/mne-tools/mne-python/jobs/608773104?utm_medium=notification&utm_source=github_status

So far the example looks perhaps a tiny bit better than master (?):

https://16649-1301584-gh.circle-artifacts.com/0/dev/auto_examples/inverse/plot_dics_source_power.html

https://mne.tools/dev/auto_examples/inverse/plot_dics_source_power.html

Great to see this progress @wmvanvliet !

britta-wstnr · 2019-11-19T09:21:22Z

So there is no big difference in this case whether we whiten the lead field or not? I guess then the sparse thing to do, also in accordance with beamformer literature (Sekihara et al., 2008 https://ieeexplore.ieee.org/document/4454055 ), is to drop whitening of the lead field and data.
I will try to check this out with the LCMV as soon as I have time.

cc @sarangnemo

wmvanvliet · 2019-11-19T09:36:22Z

I have tested the whitening also on some other data we have internally (cannot share) that is more noisy and with smaller effect sizes. For that case, whitening the leadfield "explodes" the result. Not whitening the leadfield gives sensible results. Evidence is accumulating against whitening the leadfield...

wmvanvliet · 2019-11-19T09:38:36Z

... I still don't completely understand why whitening the leadfield would be bad in beamforming but good in MNE. Could it be because MNE consideres the leadfield "as a whole", while a beamformer only operates on 3 columns at a time?

agramfort · 2019-11-19T13:15:08Z

is the whitener consistently applied to every object? data cov, CSD matrices, noise cov and leadfield? if you do a simulation with some random coloring the pipeline with whitening should be invariant. I would do something like this to check the implementation as I don't understand why it may be a bad idea to whiten

…

wmvanvliet · 2019-11-19T14:16:58Z

What do you mean by "simulation with some random coloring"? I agree that this implementation needs to be tested. I think I apply whitening constently to all objects, but we need to make sure.

agramfort · 2019-11-19T21:12:16Z

apply a full rank random matrix to the data and the leadfield and see if it affects the results. It should not if whitening is properly implemented.

…

wmvanvliet · 2019-11-20T07:31:07Z

ah, I see. Good idea!

wmvanvliet · 2020-01-06T08:28:44Z

I think I finally know where things went wrong! I forgot to whiten the CSD in the apply_dics_csd function. (We also forgot to do it in the apply_lcmv_cov function). Now, the whitening gives much more sane results.

Check out this gist that multiplied the first 50 channels by 100 to make them artificially super loud. The whitening neatly undoes this:

https://gist.github.com/wmvanvliet/bf27cbf56a176651205789d1507870b6

agramfort · 2020-01-06T13:55:42Z

you want to MRG this one then? good for review?

britta-wstnr · 2020-01-06T14:03:11Z

So, @wmvanvliet, @ckiefer0, and me have been discussing this a few times over the last weeks off Github ... Our conclusion from today: we should take the same whitener (i.e., same noise cov matrix) over the different frequencies to prevent having different forward models per frequency (the motivation that kicked this whole discussion off).

agramfort · 2020-01-06T14:04:55Z

sounds reasonable.

larsoner · 2020-01-06T14:29:33Z

Agreed this seems reasonable.

So instead of having noise_csd as in this PR, you would just have noise_cov as an argument, right?

wmvanvliet · 2020-01-07T17:10:57Z

Definitely not ready for merging yet. Lots of work still to be done to whip this into shape API wise. Of course, comments are always welcome!

Regarding the noise CSD, it is probably best if it were computed using exactly the same method as the data CSD. API-wise, I was thinking of putting in a check that the same frequencies are defined for both CSDs. Internally, we will average across the frequencies to obtain a single noise CSD.

Still some things to think about:

What happens when noise CSD is computed using less data than data CSD? I know COV matrices are scaled according to the number of averages, and CSD has n_fft field for this purpose.
Could it occur that one of the noise CSDs has a different rank than the data CSDs? Would this be cause for alarm?

larsoner · 2020-06-01T16:11:04Z

We also forgot to do it in the apply_lcmv_cov function

This is now addressed in #7656 with a test that you can manually left- and right-apply_lcmv and obtain the same result

wmvanvliet · 2020-07-09T11:05:15Z

Closing in favor of #7656

wmvanvliet added 3 commits November 1, 2019 12:49

WIP: DICS whitening

7c53b5b

Hacky way of making things work

34a1578

Undo broken equalize_channels implementation

56320f4

larsoner reviewed Nov 5, 2019

View reviewed changes

britta-wstnr reviewed Nov 6, 2019

View reviewed changes

Update mne/beamformer/_dics.py

1e65b19

Co-Authored-By: Eric Larson <larson.eric.d@gmail.com>

wmvanvliet and others added 2 commits November 6, 2019 15:50

Update mne/time_frequency/csd.py

275f68b

Co-Authored-By: Eric Larson <larson.eric.d@gmail.com>

remove stray line

46efd8d

wmvanvliet mentioned this pull request Nov 7, 2019

[WIP] Make equalize_channels capable of operating on Forward, Covariance and CrossSpectralDensity #7029

Merged

6 tasks

wmvanvliet added 3 commits November 7, 2019 15:11

Make DICS example work again, now with whitening

da95615

Also whiten the CSD

9965058

Also whiten the CSD

444ec4c

britta-wstnr mentioned this pull request Nov 8, 2019

ENH: agenda for beamformer module #5365

Open

15 tasks

wmvanvliet added 4 commits November 11, 2019 09:58

Make unit tests pass again

0dc4b00

Undo

b112081

Be careful not to use a 32-bit value as dtype

6556ff6

PEP8

a132bb3

Fix the whitening!

60066c8

wmvanvliet added 3 commits January 8, 2020 13:05

Only have one noise CSD

ce8bd08

Bit of cleanup

1ffebaf

Merge branch 'master' into dics_whitening

d7d792d

larsoner mentioned this pull request Jun 3, 2020

MRG, BUG: Refactor beamformer code #7656

Merged

13 tasks

wmvanvliet closed this Jul 9, 2020

wmvanvliet deleted the dics_whitening branch October 18, 2021 06:06


		# Set picks, use a single sensor type
		picks = mne.pick_types(raw.info, meg='grad', exclude='bads')

Uh oh!

[WIP] Add CSD whitening to DICS beamformer #7021

[WIP] Add CSD whitening to DICS beamformer #7021

Uh oh!

Conversation

wmvanvliet commented Nov 5, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

larsoner left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

britta-wstnr left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

wmvanvliet Nov 6, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

codecov bot commented Nov 6, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

wmvanvliet commented Nov 7, 2019

Uh oh!

larsoner commented Nov 7, 2019

Uh oh!

britta-wstnr commented Nov 19, 2019

Uh oh!

wmvanvliet commented Nov 19, 2019

Uh oh!

wmvanvliet commented Nov 19, 2019 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

agramfort commented Nov 19, 2019 via email

Uh oh!

wmvanvliet commented Nov 19, 2019

Uh oh!

agramfort commented Nov 19, 2019 via email

Uh oh!

wmvanvliet commented Nov 20, 2019

Uh oh!

wmvanvliet commented Jan 6, 2020

Uh oh!

agramfort commented Jan 6, 2020 via email

Uh oh!

britta-wstnr commented Jan 6, 2020

Uh oh!

agramfort commented Jan 6, 2020 via email

Uh oh!

larsoner commented Jan 6, 2020

Uh oh!

wmvanvliet commented Jan 7, 2020

wmvanvliet commented Nov 5, 2019 •

edited

Loading

wmvanvliet Nov 6, 2019 •

edited

Loading

codecov bot commented Nov 6, 2019 •

edited

Loading

wmvanvliet commented Nov 19, 2019 •

edited

Loading

larsoner commented Jun 1, 2020 •

edited

Loading