Codecov Report

Merging #5253 into master will increase coverage by <.01%.
The diff coverage is 86.66%.

@@            Coverage Diff             @@
##           master    #5253      +/-   ##
==========================================
+ Coverage   88.13%   88.13%   +<.01%     
==========================================
  Files         358      358              
  Lines       65820    65875      +55     
  Branches    11206    11210       +4     
==========================================
+ Hits        58008    58061      +53     
- Misses       4985     4987       +2     
  Partials     2827     2827

@agramfort can you check and see if you agree this could be useful? If so I'll finish. If not then I'll try to make a more compelling case before getting into the rest of the cleanup steps.

can you show me quantified evidence of the gain you observed? sorry if you already did but I find no trace of it.

I only skimmed over this, but the idea is to take the ExG events and model them as "ERP" events with XDawn?

Yeah

Interesting. Have you benchmarked the 3+ different versions of the same algorithm we have here? You could presumably also do it with the encoder module, or linear_regression_raw, right?

I have not, my assumption was that they were roughly equivalent

@jona-sassenhagen do you expect rERP / XDawn to give a better estimate for removal? What is the best test case to show it? I've tried to come up with one here but it's not entirely convincing.

It depends on what you need - if you only need the encoding models, rERP might be faster. If you need the stuff beyond that, you got to use XDawn anyways.

My question is more: how do we show that, say, the ECG waveform estimated using Xdawn (or rERP) is better than that obtained by the simple trial average? Under what conditions is it better?

... If it only happens when events overlap -- like I suspect might be the case for rERP? -- then it's not so useful for EOG / ECG. I think Xdawn has an additional noise-suppression step (the eigh call) that might help, though, not sure...

Hm, I may still be missing what you're saying.

The regression-based approach and averaging are the same if no overlap is present. You actually helped me write the CI test where we test that :)

But I could imagine the overlap correcting actually helps?

The part of XDawn after the regression, i.e.,

But I could imagine the overlap correcting actually helps?

For blinks maybe, but likely only in cases of a lot of blinking (i.e., rare). For ECG, I hope the complexes do not typically overlap. :)

So then the noise cancellation eigh step might be the important one. I was only able to show a modest benefit w/regularization here. @agramfort do you have an idea where I should see the eigh step of Xdawn improve estimate of the spatial pattern? I tried adding some additional noise / artifacts that I hoped would dominate the spatial pattern without the "division" provided by the eigh call and wasn't very successful, but maybe I did it wrong.

For blinks maybe, but likely only in cases of a lot of blinking (i.e., rare). For ECG, I hope the complexes do not typically overlap. :)

Well then it depends on the length of the window, right? If one epoch captures both the beginning of the next and the end of the preceding one, ...

Btw small note, I used cupy's dot products and linalg.solve for the solver yesterday and it gave me approx a 5x speed increase for the ridge here. So it might in principle be worth investigating that for all linear regressions we use. I'm sure it wouldn't even be hard to use skcuda.linalg here.

There are many linalg operations that slow us down. It would be great to try replacing some of them. There was no suitable replacement when I looked a couple of years ago.

cupy is in many ways as drop-in as it works -- as evidenced by the fact that even I can use it. Check it out: https://cupy.chainer.org

A sparse X, dense Y solver is just:

import cupy as cp

def cuda_solver(X, y):
    X = cp.sparse.csc_matrix(X)
    return cp.asnumpy(cp.linalg.solve((X.T.dot(X) +  cp.sparse.eye(X.shape[1])).toarray(),
                                      X.T.dot(cp.array(y)))).T

And this is cool:

https://docs-cupy.chainer.org/en/stable/reference/generated/cupy.get_array_module.html#cupy.get_array_module

I.e., write generic functions that chose cupy/CUDA or numpy as the backend depending on the input.

do the projs look visually different in sample?

Not that different from average=True.

what was the motivation for this?

If it's better to use Xdawn than SSP to find signals of interest, it should be better at finding artifacts of interest, too.

did you observe in one dataset that plain SSP was failing?

Originally I did, but I was using average=False. With average=True the difference is not as large, and in general the results seem pretty similar in the cases I've tested. However, I'm assuming there are cases where Xdawn does better than SSP -- otherwise we probably shouldn't have it as a signal tool, either, and should just have SSP. So at the least it's a consistency issue at this point. But perhaps it's not worth adding Xdawn in more places (like artifact processing) until we find a really convincing use case? Or maybe it's not so useful for artifacts that have high "SNR" anyway?

A meta-argument: it's probably good to set up MNE highly modular - so that it is possible to exchange one encapsulated part of a pipeline with another. So ideally, given that Xdawn is already in, it should be very low cost to hook it in here.

Okay let's close this one and if there is some problematic use case in the future we can see if it works better