Implementation of Block Covariance estimation

[gabelstein]

This includes block diagonal covariance estimation.
Useful for tasks where multiple band-pass filters are applied to the original data and covariance between frequencies do not hold information.

[sylvchev]

Thanks for the PR, this is a nice code. There is some linting errors in CI, that you coud see by running flake8 on your files.

Could you update what's new page in doc/whatsnew.rst?
Could it be interesting to update the examples/SSVEP/plot_classify_ssvep_mdm.py with this new estimator?

[gabelstein]

Ah yes, sorry for not checking before.
I'll check the examples and see if it makes sense.

[qbarthelemy]

Obviously, block covariance processing is interesting after filter-bank, like in SSVEP paradigm. But I think that the suggested implementation is sub-optimal.

After a filter-bank with non-overlapped frequency bands, you usually build an extended version of the input signal with filter outputs concatenated along channel axis.
The resulting matrix is block-diagonal, because the off-diagonal blocks are almost null as there is a null inter-frequency covariance.

S = [[A, 0, 0],
     [0, B, 0],
     [0, 0, C]]

However, this type of matrices has a uselessly too high dimension, leading to:
ill-conditioned covariance matrices,
larger memory storage,
and longer computations (complexity of EVD is cubic with the number of channels).

If you consider two block-diagonal matrices:

S1 = [[A1, 0, 0],            and            S2 = [[A2, 0, 0],
      [0, B1, 0],                                 [0, B2, 0],
      [0, 0, C1]]                                 [0, 0, C2]]

then, you can see that:

• the squared distance between S1 and S2 is simply the sum of squared distance between diagonal blocks:

d^2(S1,S2) = d^2(A1, A2) + d^2(B1, B2) + d^2(C1, C2)

processing only low-dimensional matrices.

• the mean is the concatenation of means computed on each blocks.

My two cents:

• conceptually, block covariance matrix S could be simply stored as [A, B, C], and processed accordingly;
• concretly, I have no idea how to implement in pyRiemann a kind of "block-diagonal matrix frontend view", being actually processed by a "by-block optimized backend computation".

I hope that’s clear...

[gabelstein]

I had thought about that as well, but I couldn't find a reasonable solution that didn't include writing and own BlockCovariance class and rewriting all things like eigenvalue decomposition etc for that.

Interestingly, although it is stored as a n*n matrix, the runtime for EVD and such reduces to about a half from the full covariances (at least in some toy examples that I wrote up).

So although this is definitely far from optimal code, it could have some usefulness. But if you think that this shouldn't be included like this, I could try to figure out a way to make it more efficient.

[sylvchev]

Yes, apart from the heavy solution to write dedicated classes for SPD and blockSPD matrices, which increase code complexity for an unclear efficiency gain, I see another possibility. This is to use scipy sparse matrices and there is even block diagonal matrix structure. But the efficiency gain is still unclear as most of our codebase rely on numpy's linear algebra, so the whole codebase requires a major rewriting.

I agree with your general remark @qbarthelemy but I think we could use the code of this PR as it is, especially if this could save some runtime when performing EVD.

[sylvchev]

@gabelstein This seems good to me one you fixed the line length for the CI; ok to merge after that?

[gabelstein]

@gabelstein This seems good to me; ok to merge?

Jup, looks good.

[sylvchev]

Could you correct the CI error?

[qbarthelemy]

Ok for me to integrate the current implementation.

[qbarthelemy]

@gabelstein , can you check last modifications?

[gabelstein]

@gabelstein , can you check last modifications?

Everything looks good.

[qbarthelemy]

Thx @gabelstein for the contribution!