TRCA variation

[ludovicdmt]

Hello,

In addition to the original method of TRCA, I propose in this PR a variation that use regularization and riemannian mean for the computation of the matrix S (instead of euclidean mean as in original implementation). It provides improvements when data are noisy and have access to large number of calibration data.

This variation is inspired by a similar work from A. Barachant on CSP: https://hal.archives-ouvertes.fr/hal-00602686/document

Ideas are :

• Empirical estimator for covariance matrix used in original paper is unbiased but with high variance especially with a high number of channel. It could benefit of some regularization.
• The matrix S in TRCA method represent the average (euclidean) inter-trial covariance matrix. Covariance matrices are semi-definite positive and lie on a riemannian manifold. A measure of distance (geodesic distance) is defined on this curvated manifold: the riemannian mean. Euclidean mean does not respect the curvature of the space.
  However the quality of the estimation of the riemannian mean depends on the number of covariance matrices used.

It does not provide an improvement of performance for the specific dataset used in ./testd/data/trcadata.mat but we notice at least 10% of improvement on our data that are more noisy and use more calibration data (about 10 trial per class).

This variation rely on Pyriemann toolbox which is already in requirements.

[codecov]

Codecov Report

Merging #39 (40735c4) into master (4aa4ba4) will increase coverage by 0.33%.
The diff coverage is 95.52%.

@@            Coverage Diff             @@
##           master      #39      +/-   ##
==========================================
+ Coverage   77.94%   78.27%   +0.33%     
==========================================
  Files          20       20              
  Lines        2158     2200      +42     
==========================================
+ Hits         1682     1722      +40     
- Misses        476      478       +2     

Impacted Files	Coverage Δ
meegkit/utils/trca.py	93.75% <93.33%> (+2.57%)	⬆️
meegkit/trca.py	95.87% <96.15%> (-1.23%)	⬇️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 4aa4ba4...40735c4. Read the comment docs.

[nbara]

Hi @ludovicdmt, that's great thanks a lot! I've left a few comments if you could please address them?

If you don't have time I'll do it myself later this week. I'll also try to improve the example_trca.py a bit

[ludovicdmt]

Thanks for the feedback, I'll have a look and do that tomorrow.

[ludovicdmt]

Should be good now

[nbara]

I'm a bit concerned by the very poor performance compared to the original method.

A performance drop can be expected, but in this cas we sometimes get ~30% accuracy on the example dataset.

The tests fail because I am testing all combinations of ensemble, method and regularization.

edit: Actually the accuracy drop is also visible even when always the euclidean mean (mean_covariance(S, metric='logeuclid')) with scm estimator, which in principle should be close to the original implementation.

edit2: actually metric='euclid with scm works ok-ish

[ludovicdmt]

Regarding performance, as explained in 1 geodesic mean is more robust to outliers than the euclidean mean. But it is also really sensible to ill-conditioned covariance matrices.
In this dataset, data are especially clean (an accuracy of 100% on some blocks). In addition, they are only 4 samples per class to estimate the $6\times6$ S matrix, so probably ill-conditioned estimations of covariance matrix S.

On data with smaller SNR and with more training samples, this TRCA variation gave us ~10% of performance improvement. I think the two methods are complementary in terms of use cases.

But maybe that's not the point of Meegkit toolbox which is more about collecting standard methods and I can't definitely understand that.

[nbara]

The fact that the riemann method doesn't work as well on this dataset is a problem in itself, it's just that it makes it difficult to unit-test. For the original implementation, I can be reasonably sure that my code is correct because it yields good classification results. Here it's a bit more complex, and it will make catching potential future bugs (due to dependency changes, etc.) more difficult to track

[nbara]

OK I restricted the test scope, and added an illustration to the example. Thanks @ludovicdmt

[ludovicdmt]

Ok thank you for the help and work @nbara !