CSP component order selection

[mbillingr]

Reference issue

Resolves #8074.

What does this implement/fix?

1. Add new parameter to CSP to select how to order components.
2. Add additional tests to make sure nothing broke.
3. (optional) Refactoring
   • Moved duplicated data type check into _check_Xy.
   • Split .fit() into smaller functions for improved readability.

Additional information

I left the refactoring in for now because I already made it.
You may to prefer not to have these changes in order to keep the code closer to the original from pyRiemann.
Let me know :)

[cbrnr]

I'll take a detailed look soon, but here are already three questions:

1. I don't quite understand the difference between 'new' and 'mutual_info' - aren't both methods sorting by mutual information? Would it make sense to combine these methods into one?
2. Why do we need a None value for the component_order parameter?
3. Does the behavior of CSP change with this PR? If so, we need a deprecation cycle, but I'd prefer if the default behavior doesn't change.

[mbillingr]

1. 'new' sorts by abs(eig - 0.5) and 'mutual_info' sorts by mutual information. In the two-class case they seem to result in the same order but in the multi-class case you can't use 'new'. I'd rather not combine them because pyRiemann distinguishes these two cases.

2. None is the current behavior, which remains the default. It chooses between 'new' and 'mutual_info' based on the number of classes. Maybe 'auto' would feel more consistent than None?

3. The default behavior of CSP does not change. This PR does not touch any of the existing tests and they still pass.

[larsoner]

'new'

I would maybe use extrema or extreme or deviation or something, because it's more about which are closer to the extremes of the possible values of the eigh result (0 and 1).

Also if mutual_info is expected to generally be better, we should just make it the default and call it a bug fix

[cbrnr]

There are too many choices for my taste. Previously, we didn't have a parameter to select the sorting order - the function used MI-based sorting for both two- and multi-class cases, even though the implementation differed. The problem with that was that the "classic" approach of sorting patterns by alternating sorted eigenvalues was not available. IMO, it would be sufficient to have two possible values for the parameter, "mutual_info" and "alternating" or something like that.

[agramfort]

I am not sure about the name alternating I started reading too. Looking at https://www.sciencedirect.com/topics/engineering/common-spatial-pattern See sentence: "Since λ1j+λ2j=1, a high value for λ1j means that the filter output based on filter *wj* yields a high variance for input signals in class 1 and a low variance for signals in class 2 (and vice versa); spatial filtering with such filters can thus significantly enhance discrimination ability. " this suggests that you should sort by the distance to 0.5. I checked the maths and I agree. I don't know where the alternating option comes from. So this suggest we need 2 options : 'variance' (dist to 0.5), 'mutual_information' and eventually a third option which is 'alternating'. I would like to see a ref in the docstring for the 3 strategies so it does not seem we came up with that in MNE... thx

…

[cbrnr]

@agramfort see my comment here.

Currently we have no options, but we're already doing two different things for the two- and multi-class case. Why do you want to introduce a new argument for something we already did without?

Also, am I correct in assuming that the distance to 0.5 is not equal (not even conceptually) to sorting by mutual information?

Even if these are not equal, I'd still summarize them as one method (the advanced), whereas the alternating one (the naive) should be there for practical reasons (it has been published, and many people in the BCI community have been using that - not sure if this is still the case though).

[agramfort]

Currently we have no options, but we're already doing two different things for the two- and multi-class case. Why do you want to introduce a new argument for something we already did without?

the question is do we want mutual info in the 2d case? if so we need a parameter.

Also, am I correct in assuming that the distance to 0.5 is not equal (not even conceptually) to sorting by mutual information?

correct.

Even if these are not equal, I'd still summarize them as one method (the advanced), whereas the alternating one (the naive) should be there for practical reasons (it has been published, and many people in the BCI community have been using that - not sure if this is still the case though).

so we still want to support the alternating option? It's your call if you think it's useful

[cbrnr]

the question is do we want mutual info in the 2d case? if so we need a parameter.

I don't think so, this PR just adds the classic alternating sorting. I think it's still useful because the question how we sort using the new method(s) comes up frequently. However, we could also go with @larsoner's suggestion and just implement the new approach. This would mean no new parameter, no change to the code, just add links to the literature to explain the current sorting. I'm fine with either approach, but I am 👎 on introducing a new parameter with three or four possible values.

[agramfort]

@cbrnr so what you suggest in the 2d case? what is called old or new here?

[cbrnr]

@agramfort I think your comment about the distance to 0.5 equalling variance is not quite correct. After all, CSP is designed to maximize the variance difference between/within conditions, so all criteria will somehow work on the variance (at least indirectly).

Therefore, I suggest that we keep things simple (at the risk of being inaccurate) and use the following values for the component_order parameter:

1. mutual_info: this is the current behavior (and the default); it uses the distance from 0.5 in the 2D case and MI in the ND case.
2. alternate: this is the classic Blankertz 2008 behavior that can be enabled for the 2D case. If this argument is used for more than 2 classes, we throw an error.

[cbrnr]

In addition, I would rename component_order to order_components or order_filters.

[agramfort]

I find the term mutual_info for the distance to 0.5 a bit misleading but I will not fight :)

…

[cbrnr]

That's the inaccuracy I was talking about 😄 - I think it's more useful to keep it simple, and people can always read the docstring to find the exact methods with references. Also, we still haven't ruled out that the distance to 0.5 is somehow related to the mutual information as a special case. @alexandrebarachant maybe?

[agramfort]

fair enough...

[mbillingr]

use the following values for the component_order parameter:

1. mutual_info: this is the current behavior (and the default); it uses the distance from 0.5 in the 2D case and MI in the ND case.
2. alternate: this is the classic Blankertz 2008 behavior that can be enabled for the 2D case. If this argument is used for more than 2 classes, we throw an error.

Does everybody agree now it should be done like that?

[cbrnr]

@mbillingr and @agramfort? Plus I suggest using order_components - also OK with you?

[mbillingr]

So we'll have order_components := 'mutual_info' | 'alternate'. OK with me.

[agramfort]

fine with me

…

[cbrnr]

I implemented the changes, let's see if tests still come back green.

[cbrnr]

@larsoner I've implemented all suggestions, please check if the examples are OK and if the footbibliography works.

[agramfort]

@larsoner merge if happy

[larsoner]

Scores changed slightly in this one but it seems okay:

• https://21998-1301584-gh.circle-artifacts.com/0/dev/auto_examples/decoding/plot_decoding_csp_timefreq.html
• https://mne.tools/dev/auto_examples/decoding/plot_decoding_csp_timefreq.html

[cbrnr]

The title "Decoding in time-frequency space data using the Common Spatial Pattern (CSP)" could be shortened to "Decoding in time-frequency space using Common Spatial Patterns (CSP)".

[larsoner]

Feel free to push this change if you want

[cbrnr]

By "scores changed slightly" I'm assuming you mean coverage - where is our coverage CI?

[larsoner]

By "scores changed slightly" I'm assuming you mean coverage

No I mean the decoding scores in the bar plots. master:

PR:

For coverage, codecov comments inline (eventually) when there are un-covered lines.

[cbrnr]

Got it. The scores changed quite a lot - @mbillingr do you have an idea what could have caused this change? I thought that this PR didn't change the previous (default) behavior?

[mbillingr]

Run-to-run variability due to cross-validation shuffling?

[larsoner]

Could be -- @cbrnr or @mbillingr please check and if so set the random state so it's reproducible

[larsoner]

Indeed there is a missing random_state. Setting it to 42 I get the same image on master and this PR:

[larsoner]

@mbillingr @cbrnr !