Igor Carrara
email: igor.carrara@inria.fr

Théodore Papadopoulo
email: theodore.papadopoulo@inria.fr

CLASSIFICATION OF BCI-EEG BASED ON THE AUGMENTED COVARIANCE MATRIX

Keywords: Augmented Covariance, Autoregressive Models, Brain-computer interface, Riemannian geometry, Takens's Theorem

Objective: Electroencephalography signals are recorded as multidimensional datasets. We propose a new framework based on the augmented covariance that stems from an autoregressive model to improve motor imagery classification. Methods: From the autoregressive model can be derived the Yule-Walker equations, which show the emergence of a symmetric positive definite matrix: the augmented covariance matrix. The state-of the art for classifying covariance matrices is based on Riemannian Geometry. A fairly natural idea is therefore to apply this Riemannian Geometry based approach to these augmented covariance matrices. The methodology for creating the augmented covariance matrix shows a natural connection with the delay embedding theorem proposed by Takens for dynamical systems. Such an embedding method is based on the knowledge of two parameters: the delay and the embedding dimension, respectively related to the lag and the order of the autoregressive model. This approach provides new methods to compute the hyper-parameters in addition to standard grid search. Results: The augmented covariance matrix performed ACMs better than any state-of-the-art methods. We will test our approach on several datasets and several subjects using the MOABB framework, using both within-session and cross-session evaluation. Conclusion: The improvement in results is due to the fact that the augmented covariance matrix incorporates not only spatial but also temporal information. As such, it contains information on the nonlinear components of the signal through the embedding procedure, which allows the leveraging of dynamical systems algorithms. Significance: These results extend the concepts and the results of the Riemannian distance based classification algorithm.

Introduction

Electroencephalography (EEG) is a modality that allows the passive measurement outside of the head of the electrical potential arising mainly from the post-synaptic electrical activity of the brain. The measurements obtained at sensors at time t are linear combinations of the activities at this same time t of a set of electrical "sources" located at the level of the cerebral cortex. These sources are time dependent on each other according to a complex biophysical model that is poorly understood (because it depends on many parameters) and which is therefore not easy to establish. There are, however, simple models which, even if they are not completely justified by biology, are general enough to be used with EEG measurements. An example of those are autoregressive models (AR). These models are based on the hypothesis that the signal can be explained by a combination of its past values plus some random factor called innovation. The AR model has been used for EEG on several occasions and in different ways. For example, they allow to model the interactions between the different zones of the brain during cognitive tasks [START_REF] Belaoucha | Structural connectivity to reconstruct brain activation and effective connectivity between brain regions[END_REF], to extract information from signals [START_REF] Ding | Short-window spectral analysis of cortical event-related potentials by adaptive multivariate autoregressive modeling: data preprocessing, model validation, and variability assessment[END_REF] or detect state changes [START_REF] Kirch | Detection of changes in multivariate time series with application to EEG data[END_REF]. Finally, autoregressive models are sometimes used in conjunction with machine learning methods for performing classification tasks [START_REF] Faradji | A self-paced two-state mental taskbased brain-computer interface with few EEG channels[END_REF].

One promising application of EEG signal analysis is Brain Computer Interfaces (BCI). BCI can be defined as a technology that measures brain activity and translates that signal into instructions or commands for a digital system. We focus on electroencephalography-based BCI (BCI-EEG), which is non-invasive, has a high temporal resolution and is relatively inexpensive. Because of these characteristics, BCI-EEG interfaces are candidates to becoming the main BCI technology in a mass deployment perspective for everyday use. In particular, we focus on motor imagery BCI paradigms, where the user changes his brain activity by imagining the movement of a body part.

Even though BCI is often considered for medical purposes, the use of this technology has expanded, and found its way into non-medical applications ranging from drone control [START_REF] Xu Duan | Quadcopter flight control using a non-invasive multi-modal brain computer interface[END_REF], rehabilitation in neurological diseases [START_REF] Lorach | Walking naturally after spinal cord injury using a brain-spine interface[END_REF] to applications in the context of Virtual and Augmented reality [START_REF] Wen | Combining brain-computer interface and virtual reality for rehabilitation in neurological diseases: A narrative review[END_REF].

Over the past 10 years, the Riemann distance-based classification paradigm has shown interesting performances for BCI-EEG classification [START_REF] Barachant | Riemannian geometry applied to BCI classification[END_REF]. The key idea behind this method is to map covariance matrices into an appropriate geometrical space: the space of Symmetric Positive Definite (SPD) matrices, which is actually a differentiable manifold with a natural Riemann structure [START_REF] Förstner | A metric for covariance matrices[END_REF]. It is therefore natural to try to extend this Riemannian approach with ideas from auto-regressive models in order to optimize the performance even for the most difficult subjects.

The goal of this work is to improve the performance of the classification algorithms, exploiting the nonlinear components of the EEG signal and combining this approach with Riemann distance-based classifiers. The idea is thus to push further the logic of the Riemannian approach by defining an appropriate SPD matrix that maximizes the amount of information related to a particular task, contained in the EEG signal.

The computation of the AR model is achieved by solving the Yule-Walker equations, which involve a SPD matrix that we call the augmented covariance matrix with lags. Thus, a fairly natural idea is to extend the standard Riemannian approach by using these augmented covariance matrices with lags.

This approach can also be seen from the point of view of dynamical systems: the AR Yule-Walker matrix can be seen as a standard covariance matrix computed on a particular embedding of the original system in a high dimensional space. The delay embedding theorem proposed by Takens [START_REF] Takens | Detecting strange attractors in turbulence[END_REF] in the context of non-linear dynamical systems establishes that such a matrix incorporates the nonlinear properties of the EEG signal provided the embedding dimension is high enough. Most EEG classification algorithms are based on the assumption that the signal can be described through a linear model. However, studies have shown that the EEG signal follows a nonlinear model [START_REF] Martin | Non-linearity in invasive EEG recordings from patients with temporal lobe epilepsy[END_REF] and because of this nonlinear component, the use of signal analysis techniques based on linearity assumptions could lose important characteristics. Only recently some classification algorithms have been proposed based on the use of quantities defined in the context of nonlinear systems theory such as correlation dimension, Lyapunov exponents, mutual information and the minimum embedding dimension [START_REF] Fang | Extracting features from phase space of EEG signals in brain-computer interfaces[END_REF][START_REF] Hosseinifard | Classifying depression patients and normal subjects using machine learning techniques and nonlinear features from EEG signal[END_REF]. Most of these non-linear features, however, are extremely computationally demanding.

The basic idea of our approach is to embed the EEG data in a high dimensional space using the theory of non-linear dynamics, and then to classify the EEG signal using the Riemannian approach. In this sense, we are thus combining the two strengths of both approaches, the use of the nonlinear component of the signal and the performance of the Riemannian distance-based algorithm.

This idea of using an augmented covariance matrix with delays was also recently introduced [START_REF] Sosulski | Introducing Block-Toeplitz Covariance Matrices to Remaster Linear Discriminant Analysis for Event-related Potential Brain-computer Interfaces[END_REF] where it was used in combination with a Linear Discriminant Analysis algorithm and applied on Event Related Potential datasets. Another use of the augmented covariance matrix methodology is outlined in a recently published paper [START_REF] Zhou | Multivariate phase space reconstruction and Riemannian manifold for sleep stage classification[END_REF]. Interestingly, this approach is specifically tailored for the domain of EEG sleep stage classification.

In order to validate the methods, we apply our approach both with a classification on the Riemann surface using the Minimum Distance to the Mean and another classification on the Tangent Space using a Support Vector Machine (SVM). We will test our approach against several pipelines of the state-of-the-art (both Machine Learning (ML) and Deep Learning (DL)), on several datasets and several subjects using the MOABB framework [START_REF] Aristimunha | Mother of all BCI Benchmarks[END_REF], with both withinsession and cross-session evaluation.

The article is structured as follows: in section 2, we first describe the theoretical approach that underlies the model, the datasets considered and the performed statistical analysis. In section 3, we list the results obtained using within-session and cross-session evaluation. Finally, in section 4, we analyze the implications of the method and its current limitations. Section 5 summarizes the results of our study.

Materials and Methods

BCI classification

The EEG signal can be represented as a multivariate time-series X ∈ R d×T where T is the number of sampled data points and d is the number of electrodes. X t ∈ R d×1 is the EEG signal at time t.

This paper tackles the problem of predicting a task achieved by a subject using its EEG signal, which is central in the field of BCI. To do so, the EEG signal is segmented into different short-time windows (epochs) corresponding to a single task. The goal is to predict the actual task from the corresponding EEG measurements.

As an example, this work is validated on motor imagery classification tasks, for which we are interested in predicting whether a subject is thinking about moving e.g. the left or the right hand (see Fig. 1). In this example, we will get some dataset consisting of epochs associated respectively to the left and right hand tasks. This dataset is split to create disjoint training and testing sub-datasets. The training sub-dataset is then used to train a classification algorithm, which is used to predict the outcome on the testing sub-dataset.

Augmented Covariance

Over the past 10 years, Riemannian approaches have shown interesting performances for SPD (Symmetric Positive Definite) matrix classification [START_REF] Barachant | Riemannian geometry applied to BCI classification[END_REF]. The key idea behind this method is to map the covariance matrix into an appropriate geometrical space: the space of SPD matrices, which is a differentiable manifold with a natural Riemann structure [START_REF] Förstner | A metric for covariance matrices[END_REF]. For more details on Riemannian geometry approach please refer to Appendix A.1.

Riemann distance based classification algorithms are traditionally using the sample spatial covariance matrix for each epoch, which is naturally SPD:

Cov(X) = 1 T ime -1 T ime i=1 X i X T i , (1)
where T denotes the transpose operator. This is not the only option for estimating the covariance, other estimators are possible, but this particular one has the advantage of being unbiased when the number of electrodes is much smaller than the temporal instants contained in an epoch.

As the covariance matrix belongs to the SPD space, it is possible to classify the different mental states using the Riemannian framework. Note that the classification can be performed either on the Riemann surface, using algorithms specifically defined on the Riemann manifold such as the Minimum Distance to the Mean (MDM), or by moving manifold points to a tangent space and using standard classification algorithms such as Support Vector Machine (SVM) or Logistic Regression [START_REF] Barachant | A Brain-Switch using Riemannian Geometry[END_REF].

The Augmented Covariance Method (ACM) is an extension of this framework with the introduction of a new SPD matrix. This new SPD matrix stems from the theory of autoregressive models. This approach consists in considering the recorded EEG signal as a realization of a multivariate random process whose law is unknown.

EEG time series are, in general, not stationary: the signal properties change over time according to the underlying mental processes. But for sufficiently short time windows where the mental task does not vary, we can consider the signal as weakly stationary, i.e. its first and second momentum are time invariant.

Under this assumption of weak stationarity, we can model the signal as a random process that follows an autoregressive model: this means that future values are expressed as a function of past ones plus a random component called innovation. Thus, the multivariate EEG signal can be modeled using an autoregressive model of order p [18]

X t = p i=1 A i X t-iτ + ε t (2)
where A i ∈ R d×d are the autoregressive coefficients , τ is a parameter called lag and ε t ∈ R d×1 is the innovation component. The order p of the process, the lag as well as the autoregressive coefficients A i are unknowns and must be estimated from the signal. The classic approach to estimate the autoregressive coefficients is using the Yule-Walker equation. In order to derive it, we simply post multiply Equation (2) by X t-k and compute its expectation value. Denoting by Γ(k) = E(X t , X T t-k ) ∈ R d×d the matrix of auto-covariance with lag k and by U the auto-covariance matrix of the innovation considered as independent from past values, we get the following equations for i ∈

[1, p] Γ(0) = p k=1 A k Γ(-kτ ) + U for i = 0 Γ(i) = p k=1 A k Γ(i -kτ ) for i ̸ = 0 (3)
The second Equation of (3) in matrix notation yields  . . . . . . . . .

Γ p-1 Γ p-2 Γ p-3 • • •               A 1 A 2 A 3 . . . A p        =        Γ 1 Γ 2 Γ 3 . . . Γ p        (4)
The autoregressive coefficients A i are thus solutions of a linear system. We call the matrix of this system Γ aug the Augmented Covariance Matrix (ACM):

Γ Aug =  Γ p-1 Γ p-2 Γ p-3 • • •        (5)
Γ aug is symmetric by construction since Γ i = Γ -i . A fairly natural idea is therefore to use these Γ aug matrices for classification algorithms using the Riemannian framework.

Γ aug combines spatial covariance with some temporal information on the signal. Γ aug is exactly the spatial covariance matrix obtained considering an augmented dataset obtained at each time t by considering p consecutive samples (separated by the fixed lag τ ). This is an embedding of our original EEG dataset in a dimension equal to d × p, with a fixed delay τ (Fig. 2). Segmenting this new augmented dataset in epochs and computing the epoch covariance matrices results in exactly the same Γ aug matrix. The sample covariance matrix is an efficient and unbiased estimator of the real covariance matrix [START_REF] Hastie | The elements of statistical learning: data mining, inference, and prediction[END_REF]. In the case of BCI applications, we are in a Finite Observation Large-Dimensional Limit condition (FOLDL) [START_REF] Bartz | Advances in high-dimensional covariance matrix estimation[END_REF]. In the case of Γ aug , the number of time sample T is kept fixed while the dimensionality increases since it is multiplied by the embedding dimension. However, in the parameter regime where we use the ACM, the FOLDL limit continues to apply.

Non Linear Dynamical approach

The formulation of the problem using the augmented dataset shows a natural connection with the Takens delay-embedding theorem [START_REF] Takens | Detecting strange attractors in turbulence[END_REF] in the context of non-linear dynamical systems.

While most EEG classification algorithms are based on the assumption that the signal is describable through a linear theory, recent studies have proven that the EEG signal follows a nonlinear model [START_REF] Martin | Non-linearity in invasive EEG recordings from patients with temporal lobe epilepsy[END_REF]. Because of this nonlinear component, the use of signal analysis techniques based on linearity assumptions could lose important characteristics, which stresses the importance of considering the theory of dynamical systems.

Let us try to understand the idea behind phase space reconstruction based on Takens's theorem. Usually what we observe in an experiment is not a phase space itself but a time series, most likely only a sequence of scalar measurements, so that not all the dynamic variables involved in the dynamics of the system under study are made available through the measurement process. So, in most cases, a time series can be seen as the values taken by the observed variables of a partially observable dynamical system. It is obvious that, even with a precise knowledge of the measurement process, it may be impossible to reconstruct the phase space of the original system from the data. Fortunately, a reconstruction of the original phase space is not really necessary for data analysis and sometimes not even desirable. It is sufficient to construct a new space that is dynamically equivalent to the original one. The process that we have just enunciated is very well summarized by Fig. 3.

Figure 3: Illustration of the concept of the Delay-Embedding theorem. Reproduced with permission from Ana González-Marcos [START_REF] De | Phase space reconstruction from a biological time series: A photoplethysmographic signal case study[END_REF] .

Considering that a process of measurement create a time series s(n), according to Takens's theorem, the geometric structure of the multi-variable dynamics of the system can be unfolded from the observable s(n) in a D-dimensional space constructed from the delay vector

s E (n) = [s(n), s(n -τ ), ..., s(n -(D -1)τ )] T ( 6
)
where τ is a positive integer called the embedding delay and D is the embedding dimension. s E (n) ∈ R D is an uniform embedding of the original phase space. As seen previously, the augmented covariance matrix is equivalent to using a standard covariance matrix on an embedding the EEG data in a high dimensional space. The approach is thus related to the theory of nonlinear dynamics but, instead of trying to identify the dynamics's parameters, we classify the mental state using the Riemann framework. In this sense, the two strengths of both approaches are combined: the use of the nonlinear component of the signal and the performance of the Riemann distance-based algorithms. Our approach turns out to be a natural extension of the Riemann distance-based classification that considers new SPD matrices containing a greater amount of information extracted from the signal based on dynamical systems theory. It is expected that the classification method based on the Riemann framework becomes even more effective by accounting for some non-linearities of the signal measured through the embedding of the data. For the sake of simplicity, we have opted for a uniform embedding approach to minimize the number of hyper-parameters. However, it's important to note that the current framework is highly flexible and can readily accommodate extensions to support non-uniform embeddings.

Considered Datasets

To test the performance of our method, we used open-access datasets available from the MOABB framework. We consider five motor imagery BCI datasets consisting of several subjects for each dataset and several sessions for each subject. On all datasets, we apply a standard band-pass filter with range [8; 35] Hz, which is classical for motor imagery tasks [START_REF] Chang S Nam | Brain-computer interfaces handbook: technological and theoretical advances[END_REF]. We consider epochs with lengths equal to the duration of task conditions which varies with the datasets under consideration (see last column of Table 1). In these datasets, all classes are balanced.

Evaluation Procedure

We consider two evaluation paradigms: Within-Session (WS) and Cross-Session (CS). The WS evaluation procedure involves evaluating performance directly within the same session of a certain subject. The classical evaluation method employed in MOABB utilized a 5-fold Cross Validation. When used with grid search hyper-parameter tuning, this was leading to data leakage. To solve this problem, we implemented a Nested Cross Validation [START_REF] Gavin | On over-fitting in model selection and subsequent selection bias in performance evaluation[END_REF] in MOABB. In our implementation, the inner loop uses a 3-fold Cross Validation strategy. The CS evaluation procedure, on the other hand, focuses on a single subject. In MOABB, this is done with a Leave One Out Cross Validation, wherein each iteration involves using a single session for testing while utilizing all other sessions for training. Figure 4 provides a schematic explanation of the experimental procedure used in this research for the grid search based algorithms. Note that the MDOP strategy, as it is not based on classification score, does not suffer of data leakage. The hyper-parameter selection is run over each training set of the outer loop and thus does not require an inner loop.

Our model depends on two hyper-parameters, the order of the autoregressive model (equivalently the embedding dimension D E ) and the delay τ . In order to compute these two parameters for every dataset and every subject, we developed two different strategies:

• Grid Search:

In order to find the hyper-parameters, we allow the grid search on two parameters to run on the domain defined by: D ∈ [1, 10] and τ ∈ [START_REF] Belaoucha | Structural connectivity to reconstruct brain activation and effective connectivity between brain regions[END_REF][START_REF] Takens | Detecting strange attractors in turbulence[END_REF]. As this is an exhaustive search, it will always provide the best possible optimum in the defined domain. The selection of hyper-parameters can be made via grid search using the nested approach or employing MDOP algorithm. Specifically, the MDOP methodology selects different sets of delays for each embedding and for each epoch and then obtains a unique value for delay and embedding dimension using an average across delays and epochs. Once these parameters are selected, we use the phase space reconstruction process to obtain a dynamic system equivalent to the original one (in this figure, we see an embedding corresponding to D E = 2 and τ = 10). From this high-dimensional space, the covariance matrix is extracted and classified using algorithms on the Riemannian manifold (MDM) or on the Tangent space (SVM) through projection.

• Non Linear Dynamical approach:

The parameters D and τ can also be estimated using the concepts of non-linear dynamical systems. There are several algorithms for estimating optimal parameters in noisy conditions that can be grouped in two categories: traditional and unified approaches. In this research, we focused on unified approaches that were implemented using the DynamicalSystem.jl library [START_REF] Datseris | DynamicalSystems.jl: A Julia software library for chaos and nonlinear dynamics[END_REF]. Specifically, we used the MDOP algorithm [START_REF] Nichkawde | Optimal state-space reconstruction using derivatives on projected manifold[END_REF], which works iteratively by finding simultaneously the best delay and embedding dimensions. The algorithm starts with the current dataset and adds a lagged version of that dataset to it by finding the correct value of τ such that the redundancy between coordinates is minimized, i.e. by maximizing the beta statistics [START_REF] Nichkawde | Optimal state-space reconstruction using derivatives on projected manifold[END_REF]. After each embedding cycle the FNN-statistic [START_REF] Matthew B Kennel | Determining embedding dimension for phase-space reconstruction using a geometrical construction[END_REF] is being checked and as soon as this statistic drops below a threshold, the algorithm terminates. It should be noted that this algorithm produces different values of τ for each embedding cycle. We decided to select the optimal value of τ as an average of these values, while we select D E as the number of cycles. This decision is rooted in our methodology's reliance on uniform embedding.

Studied Pipelines

We intend to evaluate how the augmented covariance method (ACM) as a general technique can increase the performance in different settings. We used the ACM with several Riemannian classification methods on the surface or on the tangent space as shown in Table 2. These pipelines will be denoted by the blue color in the result tables: a lighter blue is used to distinguish those working directly on the Riemannian manifold from those working in the tangent space. A visual representation of the methodology employed in this paper is available in Fig. 5. We compared our pipelines against eight state-of-the-art algorithms of the BCI literature:

• A reference method (colored in red in the result tables):

1. "CSP+LDA", a combination of the Common spatial pattern (CSP) algorithm followed by a classification performed on a shrinkage Linear Discriminant Analysis (LDA) [START_REF] Lotte | Learning from other Subjects Helps Reducing Brain-Computer Interface Calibration Time[END_REF].

• Four Riemannian methods (colored in yellow in result tables):

2. "FgMDM", a classification method by Minimum Distance to the Mean after having applied a geodesic filtering [START_REF] Yger | Riemannian approaches in brain-computer interfaces: a review[END_REF].

3. "MDM", a classification method by Minimum Distance to the Mean [START_REF] Barachant | Riemannian geometry applied to BCI classification[END_REF].

4. "COV+EN", a classification method in the tangent space using an Elastic Network as classifier [START_REF] Corsi | Functional connectivity ensemble method to enhance BCI performance (FUCONE)[END_REF].

5. "TGSP+SVM", a classification method in the tangent space using a Support Vector Machine [START_REF] Barachant | Riemannian geometry applied to BCI classification[END_REF].

• Three Deep Learning (DL) methods (colored in green in result tables):

6. "EEGNet", a DL architecture based on CNN, with standardized and re-sampled EEG signal at 128Hz [START_REF] Vernon J Lawhern | EEGNet: a compact convolutional neural network for EEG-based brain-computer interfaces[END_REF].

7. "DeepConvNet", a DL architecture based on CNN, with standardized and re-sampled EEG signal at 250Hz [START_REF] Tibor | Deep learning with convolutional neural networks for EEG decoding and visualization[END_REF].

8. "ShallowConvNet", a DL architecture based on convolutional layer (CNN), with standardized and re-sampled EEG signal at 250Hz [START_REF] Tibor | Deep learning with convolutional neural networks for EEG decoding and visualization[END_REF].

The pipeline 1 is the most basic pipeline and is given for reference. For pipelines 2-5, which are based on Riemannian geometry, all the hyper-parameters of the classifiers used in this study are optimized with grid search. Table A1 summarizes the domains over which these grid searches are performed.

The DL pipelines 6-8 were implemented using the TensorFlow [START_REF] Abadi | TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems[END_REF] and SciKeras libraries. To enhance result stability, each DL pipeline starts with a standardization step, which normalizes every channel to have a zero mean and unit standard deviation. Additionally, we employed a resampling procedure to ensure that each architecture integrates a temporal filter aligned with the state-of-the-art techniques' implementations. For specific details on the DL hyper-parameters, please refer to Table A2.

Statistical Analysis

In order to validate the best pipeline from a statistical point of view, we use the instrument provided by the MOABB framework. To remain fast, it uses a one-tailed permutation-based paired t-test [START_REF] Student | The probable error of a mean[END_REF] for datasets with less than 20 subjects, or a Wilcoxon non-parametric signedrank test [START_REF] Wilcoxon | Individual comparisons by ranking methods[END_REF] otherwise. This statistical test is used to generate a p-value that compares two pipelines for each pair of pipelines. Then, the p-values are combined using Stouffer's method [START_REF] Samuel | The american soldier: Adjustment during army life[END_REF] in order to get a final p-value for each hypothesis. In this way, we can validate the results using several subjects and several datasets.

The score are obtained using the Area Under Curve (AUC) of the Receiver Operating Characteristic (ROC) curve. AUC ranges in value from 0 to 1. A higher AUC value characterizes a better classifier. In case of a multi-class classification, the metric used is the accuracy.

Results

In this section, we report the results of the augmented covariance method over the classification on both the Riemann surface and the Tangent Space. To validate the robustness and validity of our approach, we test the algorithm on different datasets, subjects and tasks using the MOABB framework.

Right hand vs Left hand

We consider 3 different datasets -BNCI2014001, BNCI2014004 and Zhou2016 -to classify a right vs left hand task classification using within-session and cross-session evaluation procedures.

The various approaches of this study are compared in Table 3. A detailed picture of the statistical significance of these results is provided in Fig. 6 andA2. Out of the three datasets, Zhou2016 is generally giving good results with all methods. For this dataset, the best result is obtained with ShallowConvNet but is followed very closely by our best approach -ACM+TGSP+SVM (Grid) -with a non-significant decrease of 1% in the AUC score. With the other datasets, ACM+TGSP+SVM (Grid) brings a notable improvement over all other non-ACM methods with an average improvement of AUC by 4-16%. The improvements are very similar with both the WS and CS validation modes. Except in one case, all ACM methods improve on their non-ACM versions. The exceptional case is ACM+TGSP+SVM (MDOP) vs TGSP+SVM for the BNCI2014001 dataset with WS evaluation. The small regression (non significant) can be attributed to a failure of the MDOP algorithm to find optimal hyper-parameter values as the ACM+TGSP+SVM (Grid) clearly performs well in the same benchmark.

Right hand vs Feet

This classification task was performed on 4 different datasets -BNCI2014001, BNCI2014002, BNCI2015001 and Zhou2016 -to perform a right hand vs feet task classification using a withinsession and cross-session evaluation procedure. For sake of simplicity, we do not report the results for left hand vs feet task as they are very similar. There are no CS results for BNCI2014002 because this dataset is composed by only one session.

The various approaches of this study are compared in Table 4. A detailed picture of the statistical significance of these results is provided in Fig. A3 andA4. Generally speaking, the results are better than with the Left hand -Right hand case, showing that the Right hand -Feet task is generally the easiest of the two. Our best approach -ACM+TANG+SVM (Grid) -brings 0-12% of improvements over all other methods. The only method that gives some competitive results is again ShallowConvNet on the Zhou2016 dataset, which, as in the previous case, generally gives good results with all methods. Except in three cases, all ACM methods improve on their non-ACM counterparts. The exceptional cases again arise for ACM+TGSP+SVM (MDOP) vs TGSP+SVM. The MDOP algorithm seems to have more difficulty to find optimal hyper-parameter values with this specific method.

Right hand vs Left hand vs Feet

In this last case, we use 2 different datasets (BNCI2014001 and Zhou2016) to perform a three tasks classification -right hand vs left hand vs feet -using the within-session and cross-session evaluation procedures.

The various approaches considered in this study are compared in Table 5. A detailed picture of the statistical significance of these results is provided in Fig. A5 andA6. With the BNCI2014001 dataset, ACM+TGSP+SVM (Grid) brings a notable performance improvement with an average improvement of Accuracy of 6-12% over all non DL methods and of 0-37% with respect to DL methods. DeepConvNet seems to perform particularly bad on this benchmark. Removing it gives improvements of 0-17% over DL methods. ShallowConvNet is competitive for the CS evaluation method. For the Zhou2016 dataset, the best method is ShallowConvNet. The ACM+TGSP+SVM (Grid) method is a close second with a performance loss of 2-8%. This dataset appears to be noticeably more challenging for 3 tasks classification than for binary classification. As for the other tasks, ACM based methods improve on their non-ACM counterparts, with exceptions arising only with the ACM+TGSP+SVM (MDOP) algorithm.

Discussion

The numerous tests on different classification tasks, both binary and multi class, have shown some notable observations that we summarize here. In general, with either ACM or non ACM approaches, the Riemannian method based on tangent space classification is more effective than Riemann surface classification approaches (as shown in Tables 3, 4 and5). ACM-based methods consistently demonstrate enhancements over their standard covariance matrix counterparts, excelling in both Riemannian manifold and tangent space classification. This trend is more marked with the grid search hyper-parameters estimation than with that based on MDOP, which sometimes seems to have trouble to find a good set of hyper-parameters. This is counterbalanced by the fact that MDOP is generally faster than grid search, while often still providing improved results.

ACM improves performance in both WS and CS evaluation

The utilization of the ACM methodology enhances performance in both the WS and CS cases, demonstrating its stability and applicability in the more complex CS situations, which involve more variability. This is important for clinical or real-world scenarios, for which CS is a more realistic condition. This finding shows that the ACM methodology brings for some new and

ACM+TGSP+SVM (Grid) outperforms all DL pipelines. ShallowCon-vNet is sometimes competitive

Globally, deep learning methods perform less favorably than the non DL ACM ones. One exception is ShallowConvNet that sometimes is better than the best ACM method ACM+TGSP+SVM (Grid). However, the improvements it provides are never statistically significant. This success can be attributed to the specific tailored feature created using the ACM compared to the automatic feature extraction achieved by DL models.

hyper-parameter stability

The performance improvements are naturally less important when the baseline methods are performing better initially, but this definitely shows that the augmented covariance approach brings some valuable information into the classification problem for these datasets. The selection of hyper-parameters that produces the best performance is through the use of grid search, however, the modality based on non-linear dynamics concepts (MDOP) is also shown to produce in general better performance than standard algorithms (with a few exceptions). However, the selection of the hyper-parameters through grid search is time consuming in the training phase. Generally speaking, hyper-parameters values are not stable across subjects and sessions. Fig. 7 (a) shows one example of the spread of the order and lag hyper-parameters over subjects and sessions for BNCI2014001. Comparatively, it can be seen that the MDOP unified approach gives more condensed point sets. Yet, this method based on nonlinear systems theory places an important theoretical foundation for finding efficiently proper hyper-parameters which generally achieve performance improvements for both WS and CS evaluations, even if far from optimal. Fig. 7 (b) shows the set of parameters that provide an accuracy within 99% of the best one. On this example, we see that there are numerous almost equivalent sets of hyperparameters. This shows that focusing on the actual optimal value is maybe too strict. Fig. 7 (c) is a voting map which counts for each couple of hyper-parameters when they are found in this 99% region over subjects × sessions. The distribution of optimal hyper-parameters is not unimodal demonstrating that finding fixed values for these hyper-parameters that work in every situation is difficult. Other classification tasks and datasets exhibit the same behavior.

Computational time

In this section, we report the computational times of the proposed algorithms with respect to their scores. To allow a fair comparison, we have run all the algorithms on the same hardware, a Dell C6420 dual-Xeon Cascade Lake SP Gold 6240 @ 2.60GHz. The results for the BNCI2014001 dataset for the right hand vs left hand classification with WS evaluation procedure are shown in Fig. 8. The MDOP algorithms (ACM+MDM (MDOP), ACM+TGSP+SVM (MDOP)) exhibit a similar speed as ACM+MDM (Grid) while delivering slightly lower performance. On the other hand, the ACM+TANG+SVM (Grid) algorithm outshines in terms of performance but comes at the cost of using substantial computational resources. This increased computational load primarily arises from the necessity to explore a larger parameter space during grid search optimization, particularly due to the SVM component. Overall, the best trade off between algorithmic performance and computational efficiency is provided by ACM+MDM (Grid). Alternatively, opting for hyper-parameter selection through MDOP can be advantageous, especially when the classification algorithm itself requires extensive computational resources.

Finally, note that the augmented methods increase computational time because of the increased dimensionality of the augmented covariance compared to the standard covariance [START_REF] You | Geometric learning of functional brain network on the correlation manifold[END_REF], as shown in Fig. 9 for the BNCI2014001 dataset for the right hand vs feet classification with WS evaluation procedure. In this graph we show the computational performance of ACM+MDM without the grid search hyper-parameter selection but setting specifically the order. Yet, even in the higher orders, the total classification time of a new sample with the augmented approaches is compatible with real time computations.

Limitations

While our approach has been shown to achieve improvement in performance in many BCI applications in the Motor Imagery task, it has at least one limitation : as other methods, it relies on the assumption of the validity of the sample covariance estimator. When the span of the used time windows are way greater than the dimension of covariance matrix, this assumption remains valid. But this dimension increases quickly with higher AR orders. In that case, the matrix is still a symmetric positive matrix but is not definite and Riemannian metrics can no longer be used. To solve this issue, it is possible to adopt shrinkage or other methods to estimate covariance [START_REF] Ledoit | A well-conditioned estimator for large-dimensional covariance matrices[END_REF][START_REF] Bartz | Covariance shrinkage for autocorrelated data[END_REF][START_REF] Sabbagh | Manifold-regression to predict from MEG/EEG brain signals without source modeling[END_REF].

Conclusion

In this research, we explored the use of the augmented covariance matrix for the classification of Motor Imagery task in EEG-BCI. This methodology amounts to integrate temporal information into the normal spatial covariance, and improves quite dramatically classification results over the state-of-the-art when used with Riemann distance based classification algorithms. There remain some issues with the definiteness of the augmented covariance matrix, which stems from the augmentation of its dimensionality when the order increases. Another issue is the choice of the hyper-parameters which is costly when done exhaustively (grid search). Yet, we show that methods inspired by the Takens' theorem can provide some reasonable parameter values at the cost of a decrease of the improvements (but with improvements still in most cases).

The current study focused on offline classification. Our next direction will be to extend the approach for online classification. To do this, it will be necessary to test the algorithms on narrower epoch sliding windows. Another important development of the augmented approach would be to check its efficiency with smaller training dataset i.e. with either a low number of electrodes or with fewer training samples. This would provide an important usability improvement for BCIs, potentially reducing the setup and/or the calibration times. We have already seen that performance over the state-of-the-art is improved in the case of BNCI2014004, which possesses only 3 electrodes, but this needs to be generalized to other datasets and to a smaller set of training samples.

Another work direction is to analyze more in depth the distribution of hyper-parameters with the objective of identifying a methodology to make the classifier parameter free (e.g. by using ensemble classifiers methodology).

Acknowledgment

This work has been partially financed by a EUR DS4H/Neuromod fellowship. The authors are grateful to the OPAL infrastructure from Université Côte d'Azur for providing resources and support.

Data and Code Availability

The codes used to produce the results of this study are publicly available in this Github repository: https://github.com/carraraig/Augmented_Covariance_BCI. In the repository, we A Supplementary Material

A.1 Riemann geometry

Over the past 10 years, Riemannian approaches have shown interesting performances for SPD matrix classification [START_REF] Barachant | Riemannian geometry applied to BCI classification[END_REF]. The key idea behind this method is to map the covariance matrix into an appropriate geometrical space. The space of SPD matrix is a differentiable manifold with a natural Riemann structure [START_REF] Förstner | A metric for covariance matrices[END_REF]. The spaces of n × n real square and real symmetric matrices are denoted by M (n) and S(n), respectively: S(n) = S ∈ M (n), S T = S . The space of SPD matrices is defined as P (n) = {P ∈ S(n), P > 0}. In this way, we can represent a matrices P ∈ S(n) as a point on a Riemann manifold of dimension n(n + 1)/2. Since the space of SPD matrix is a manifold of non zero curvature, we cannot expect that concepts of Euclidean geometry preserve the manifold structure. It is possible to define a distance between two SPD matrix P 1 and P 2 in P (n) as the length of the geodesic that connects P 1 and P 2 on the Riemann manifold. In this paper, we use the affine-invariant metric [START_REF] Moakher | A Differential Geometric Approach to the Geometric Mean of Symmetric Positive-Definite Matrices[END_REF] δ R (P 1 , P 2 ) = log P 1 . The choice of the metric is not unique but we choose this metric because it possesses relevant properties for BCI applications: invariance under reordering, invariance under congruent transformation and invariance under inversion [START_REF] Congedo | Riemannian geometry for EEG-based brain-computer interfaces; a primer and a review[END_REF]. To define the mean of m SPD matrices (P 1 , ..., P m ), we use the concept of Fréchet mean [START_REF] Moakher | A Differential Geometric Approach to the Geometric Mean of Symmetric Positive-Definite Matrices[END_REF]

P = argmin P∈P (n) m i=1 δ 2 R (P, P i ) (8)
P is called the geometric mean in the Riemannian sense or the center of mass of (P 1 , ..., P m ). Note that, for m > 2, no closed-form expression is know, so we need to employ an iterative algorithm [START_REF] Pennec | A Riemannian Framework for Tensor Computing[END_REF]. Due to the properties of Eq. ( 7), the geometric mean inherits the properties of congruence invariance [START_REF] Congedo | Riemannian geometry for EEG-based brain-computer interfaces; a primer and a review[END_REF]. For every point P i ∈ P (n), in a complete Riemann space, we can define a point S i ∈ S(n) in the tangent space T P to the Riemann surface at point P ∈ P (n). The tangent space has zero curvature by definition and allows us to use standard ML algorithms and other classical tools. Mappings from the Riemann surface to the tangent space and from the tangent space to the Riemann surface, are provided respectively by the Exp and Log maps [START_REF] Moakher | A Differential Geometric Approach to the Geometric Mean of Symmetric Positive-Definite Matrices[END_REF] (see Fig. A1):

Exp P (S i ) = P 1/2 Exp P -1/2 S i P -1/2 P 1/2 Log P (P i ) = P 1/2 Log P -1/2 P i P -1/2 P 1/2 (9)

In the context of a Riemann geometry based classification algorithm, we estimate a spatial covariance matrix for each epoch with the standard approach of the spatial sample covariance estimator

Cov(X) = 1 T ime -1 T ime i=1 X i X T i , ( 10
)
where T is the transpose operator. This is not the only option for estimating the covariance. We could have used other estimators, but this one has the advantage of being unbiased when the number of electrodes is much smaller than the temporal instants contained in an epoch.

As the covariance matrix belongs to the SPD space, it is possible to classify the different mental states using the Riemannian framework. Note that the classification can be perform either

Figure 1 :

 1 Figure 1: Example of an EEG dataset composed of 4 electrodes. For a motor imagery task for the classification between right and left hand. The experimental paradigm consists in asking the subject to imagine the movement of a body part (left or right hand in this case), at specific times instants marked with green lines. For each green line, an epoch of data (blue boxes) is collected and associated with its left/right hand label. A subset of these epochs forms the training dataset.

Figure 2 :

 2 Figure 2: (a) Consider the EEG signal with 4 channels and an epoch (blue box), and create a new dataset by concatenating its values at time t with those at time t + τ . (b) Result of the previous procedure at every point of the dataset. This is a higher dimensional dataset for a fixed delay τ and with dimension d × p (here d = 4 and p = 2). The augmented covariance matrix of the original signal Γ aug is exactly the standard covariance matrix of this concatenated dataset.

Figure 4 :

 4 Figure4: Grid search evaluation procedures. In order to avoid data leakage, a Nested Cross Validation procedure is adopted. Plot (a) shows the within-session evaluation procedure: the outer loop uses a stratified 5-fold CV, while the inner loop, in charge of tuning the model's hyperparameters, uses a stratified 3-fold CV. Plot (b) shows the cross-session evaluation procedure, which closely resembles the within-session one, with the key distinction being that the outer loop adopts a leave-one-session-out CV strategy for testing.

Figure 5 :

 5 Figure 5: The figure schematically illustrates the ACM methodology. The presented example uses the BNCI2014004 dataset that has only 3 electrodes (in red on the top left plot). The measurement process of the original dynamic system is this a 3-dimensional time series. The process then begins with the extraction of epoch signal representing left and right hand tasks.The selection of hyper-parameters can be made via grid search using the nested approach or employing MDOP algorithm. Specifically, the MDOP methodology selects different sets of delays for each embedding and for each epoch and then obtains a unique value for delay and embedding dimension using an average across delays and epochs. Once these parameters are selected, we use the phase space reconstruction process to obtain a dynamic system equivalent to the original one (in this figure, we see an embedding corresponding to D E = 2 and τ = 10). From this high-dimensional space, the covariance matrix is extracted and classified using algorithms on the Riemannian manifold (MDM) or on the Tangent space (SVM) through projection.

Figure 6 :

 6 Figure6: Results for Right vs Left hand classification, using WS evaluation. Plot (a) provides a combined meta analysis (over all datasets) of the different pipelines. It shows the significance that the algorithm on the y-axis is better than the one on the x-axis. The color represents the significance level of the difference of accuracy, in terms of t-values. We only show significant interactions (p < 0.05). Plot (b) is a box plot of the score with the error of the different pipelines for every dataset considered. Plots (c), (d), (e), (f), (g) and (h) show the meta analysis of ACM+TGSP+SVM (Grid) against respectively Cov+EN, FgMDM, CSP+LDA, DeepConvNet, EEGNet and ShallowConvNet. We show the standardized mean differences of p-values computed as one-tailed Wilcoxon signed-rank test for the hypothesis given as title of the plot. The gray bar denotes the 95% interval. * stands for p < 0.05, ** for p < 0.01, and *** for p < 0.001.

Figure 7 :Figure 8 :

 78 Figure 7: Plot (a) shows the best hyper-parameters for the ACM+MDM algorithm for each session and subject in the left hand right hand task, using the BNCI2014001 dataset. The red box defines the grid search domain. The heat-map in plot (b) shows the area that contains the hyper-parameter that provide a classification score of 99% of the maximum score obtained. These results are computed for Subject 1 in BNCI2014001 Session E in a right hand vs left hand task. The heat-map of plot (c) shows the maximal occurrence of best area of hyper-parameters across subjects and sessions of the BNCI2014001 dataset. It is obtained by counting the number of non empty boxes of plots similar to plot (b) over all sessions and subjects.

Figure 9 :

 9 Figure9: Total computational time (training+testing) for the MOABB within-session evaluation for subject 1 of BNCI2014001 with the task right hand vs feet using ACM+MDM without the grid search hyper-parameter selection but setting specifically the order. Computational times at order 1 correspond to standard covariance. Similar behaviors occur in the other evaluation procedures, subjects, tasks and datasets.

  where ∥ • ∥ F is the Frobenius norm, Log() is the logarithm of a matrix and λ i are the eigenvalues of P

Table 1 :

 1 Table 1 contains all the details about those datasets. Dataset considered during this study

	Dataset	subjects	channels sampling rate	sessions	tasks	trials/class Epoch (s)
	BNCI2014001 [22]	9	22	250 Hz	2	4	144	[2, 6]
	BNCI2014002 [23]	14	15	512 Hz	1	2	80	[3, 8]
	BNCI2014004 [24]	9	3	250 Hz	5	2	360	[3, 7.5]
	BNCI2015001 [25]	12	13	512 Hz	3	2	200	[0, 5]
	Zhou2016 [26]	3	14	250 Hz	3	3	160	[0, 5]

Table 2 :

 2 Considered pipelines with classifications on the Riemann surface or the Tangent Space.

	Surface	Standard algorithm	Augmented algorithm
			Grid Search	Unified Approach
	Riemann Surface	MDM	ACM+MDM (Grid)	ACM+MDM (MDOP)
	Tangent Space	TGSP+SVM	ACM+TGSP+SVM (Grid) ACM+TGSP+SVM (MDOP)

Table 3 :

 3 Performance (AUC) of Right hand vs Left hand classification. The color code refers to the type of method (see section 2.6).

	Pipeline	Eval BNCI2014001 BNCI2014004	Zhou2016
	ACM+MDM (Grid)	WS	0.89 ± 0.13	0.82 ± 0.16	0.91 ± 0.07
	ACM+MDM (MDOP)	WS	0.86 ± 0.14	0.81 ± 0.16	0.91 ± 0.08
	ACM+TGSP+SVM (Grid)	WS	0.92 ± 0.11	0.82 ± 0.16	0.94 ± 0.06
	ACM+TGSP+SVM (MDOP) WS	0.87 ± 0.13	0.82 ± 0.16	0.93 ± 0.06
	CSP+LDA	WS	0.84 ± 0.15	0.79 ± 0.16	0.91 ± 0.10
	FgMDM	WS	0.88 ± 0.11	0.79 ± 0.15	0.91 ± 0.08
	MDM	WS	0.84 ± 0.14	0.78 ± 0.16	0.89 ± 0.08
	Cov+EN	WS	0.88 ± 0.12	0.80 ± 0.16	0.93 ± 0.08
	TGSP+SVM	WS	0.88 ± 0.12	0.79 ± 0.15	0.92 ± 0.08
	EEGNet	WS	0.76 ± 0.21	0.69 ± 0.20	0.94 ± 0.02
	DeepConvNet	WS	0.82 ± 0.17	0.72 ± 0.19	0.93 ± 0.07
	ShallowConvNet	WS	0.85 ± 0.15	0.72 ± 0.18	0.95 ± 0.06
	ACM+MDM (Grid)	CS	0.89 ± 0.13	0.83 ± 0.15	0.94 ± 0.04
	ACM+MDM (MDOP)	CS	0.86 ± 0.13	0.81 ± 0.15	0.93 ± 0.04
	ACM+TGSP+SVM (Grid)	CS	0.91 ± 0.10	0.83 ± 0.15	0.93 ± 0.06
	ACM+TGSP+SVM (MDOP) CS	0.87 ± 0.12	0.81 ± 0.15	0.94 ± 0.05
	CSP+LDA	CS	0.84 ± 0.15	0.81 ± 0.14	0.92 ± 0.08
	FgMDM	CS	0.85 ± 0.14	0.80 ± 0.14	0.90 ± 0.09
	MDM	CS	0.83 ± 0.15	0.79 ± 0.15	0.92 ± 0.05
	Cov+EN	CS	0.86 ± 0.13	0.81 ± 0.14	0.93 ± 0.07
	TGSP+SVM	CS	0.86 ± 0.14	0.80 ± 0.14	0.91 ± 0.07
	EEGNet	CS	0.75 ± 0.21	0.75 ± 0.17	0.95 ± 0.05
	DeepConvNet	CS	0.82 ± 0.17	0.75 ± 0.17	0.94 ± 0.07
	ShallowConvNet	CS	0.87 ± 0.12	0.73 ± 0.14	0.94 ± 0.06

Table 4 :

 4 Performance (AUC) of Right hand vs Feet classification. The color code refers to the type of method (see section 2.6).

	Pipeline	Eval BNCI2014001 BNCI2015001	Zhou2016	BNCI2014002
	ACM+MDM (Grid)	WS	0.96 ± 0.05	0.90 ± 0.11	0.94 ± 0.04	0.82 ± 0.15
	ACM+MDM (MDOP)	WS	0.93 ± 0.08	0.89 ± 0.12	0.94 ± 0.05	0.78 ± 0.17
	ACM+TGSP+SVM (Grid)	WS	0.97 ± 0.03	0.93 ± 0.08 0.98 ± 0.02 0.89 ± 0.11
	ACM+TGSP+SVM (MDOP) WS	0.94 ± 0.08	0.92 ± 0.09	0.97 ± 0.03	0.83 ± 0.16
	CSP+LDA	WS	0.92 ± 0.10	0.89 ± 0.11	0.94 ± 0.04	0.85 ± 0.13
	FgMDM	WS	0.94 ± 0.08	0.90 ± 0.10	0.96 ± 0.03	0.85 ± 0.12
	MDM	WS	0.91 ± 0.10	0.86 ± 0.13	0.92 ± 0.05	0.77 ± 0.15
	Cov+EN	WS	0.95 ± 0.07	0.91 ± 0.09	0.97 ± 0.03	0.86 ± 0.12
	TGSP+SVM	WS	0.95 ± 0.07	0.91 ± 0.09	0.96 ± 0.03	0.86 ± 0.12
	EEGNet	WS	0.89 ± 0.15	0.90 ± 0.12	0.95 ± 0.03	0.84 ± 0.17
	DeepConvNet	WS	0.88 ± 0.13	0.88 ± 0.14	0.96 ± 0.04	0.88 ± 0.11
	ShallowConvNet	WS	0.92 ± 0.10	0.91 ± 0.11	0.98 ± 0.02	0.88 ± 0.12
	ACM+MDM (Grid)	CS	0.94 ± 0.06	0.91 ± 0.09	0.95 ± 0.04
	ACM+MDM (MDOP)	CS	0.92 ± 0.09	0.90 ± 0.10	0.95 ± 0.05
	ACM+TGSP+SVM (Grid)	CS	0.96 ± 0.05	0.93 ± 0.07 0.99 ± 0.01
	ACM+TGSP+SVM (MDOP) CS	0.93 ± 0.07	0.92 ± 0.09	0.97 ± 0.04
	CSP+LDA	CS	0.89 ± 0.14	0.88 ± 0.11	0.94 ± 0.05
	FgMDM	CS	0.93 ± 0.08	0.90 ± 0.10	0.95 ± 0.06
	MDM	CS	0.90 ± 0.11	0.87 ± 0.12	0.93 ± 0.06
	Cov+EN	CS	0.94 ± 0.07	0.90 ± 0.10	0.96 ± 0.04
	TGSP+SVM	CS	0.94 ± 0.07	0.88 ± 0.12	0.96 ± 0.05
	EEGNet	CS	0.87 ± 0.16	0.90 ± 0.12	0.98 ± 0.02
	DeepConvNet	CS	0.86 ± 0.16	0.90 ± 0.12	0.97 ± 0.03
	ShallowConvNet	CS	0.90 ± 0.12	0.92 ± 0.10	0.99 ± 0.01

Table 5 :

 5 Performance (Accuracy) of Right hand vs Left hand vs Feet classification. The color code refers to the type of method (see section 2.6).

	Pipeline	Eval BNCI2014001	Zhou2016
	ACM+MDM (Grid)	WS	0.78 ± 0.12	0.79 ± 0.04
	ACM+MDM (MDOP)	WS	0.76 ± 0.14	0.79 ± 0.06
	ACM+TGSP+SVM (Grid)	WS	0.84 ± 0.11	0.84 ± 0.04
	ACM+TGSP+SVM (MDOP) WS	0.78 ± 0.15	0.82 ± 0.04
	CSP+LDA	WS	0.74 ± 0.16	0.81 ± 0.06
	FgMDM	WS	0.77 ± 0.14	0.81 ± 0.06
	MDM	WS	0.72 ± 0.15	0.75 ± 0.07
	Cov+EN	WS	0.79 ± 0.14	0.83 ± 0.06
	TGSP+SVM	WS	0.78 ± 0.14	0.83 ± 0.05
	EEGNet	WS	0.67 ± 0.20	0.84 ± 0.04
	DeepConvNet	WS	0.47 ± 0.08	0.56 ± 0.06
	ShallowConvNet	WS	0.76 ± 0.16	0.86 ± 0.04
	ACM+MDM (Grid)	CS	0.73 ± 0.12	0.72 ± 0.08
	ACM+MDM (MDOP)	CS	0.70 ± 0.13	0.72 ± 0.09
	ACM+TGSP+SVM (Grid)	CS	0.76 ± 0.11	0.78 ± 0.06
	ACM+TGSP+SVM (MDOP) CS	0.71 ± 0.14	0.76 ± 0.10
	CSP+LDA	CS	0.66 ± 0.16	0.73 ± 0.13
	FgMDM	CS	0.69 ± 0.13	0.72 ± 0.10
	MDM	CS	0.66 ± 0.13	0.70 ± 0.09
	Cov+EN	CS	0.70 ± 0.12	0.74 ± 0.10
	TGSP+SVM	CS	0.67 ± 0.15	0.73 ± 0.10
	EEGNet	CS	0.66 ± 0.21	0.84 ± 0.06
	DeepConvNet	CS	0.41 ± 0.10	0.68 ± 0.12
	ShallowConvNet	CS	0.76 ± 0.16 0.86 ± 0.05
	4.1 Tangent Space classification outperforms the manifold method, ACM
	enhances performance in both scenarios.

      Γ 0 Γ -1 Γ -2 • • • Γ 1 Γ 0 Γ -1 • • • Γ 2 Γ 1 Γ 0 • • • . . .

      Γ 0 Γ -1 Γ -2 • • • Γ 1 Γ 0 Γ -1 • • • Γ 2 Γ 1 Γ 0 • • • . . . . . . . . . . . .

ACM+MDM(Grid)

have made available the code along with the results table that includes additional scoring metrics not explicitly reported in the paper. Moreover we release the function AugmentedDataset in the MOABB library in order to facilitate the use of the ACM. Figure A1: A Riemann manifold and the maps from/to the tangent space at point P.. on the Riemann surface, using algorithms specifically defined on the Riemann manifold as Minimum Distance to the Mean (MDM), or by moving points to tangent space and using standard classification algorithms such as Support Vector Machine (SVM) or Logistic Regression [START_REF] Barachant | A Brain-Switch using Riemannian Geometry[END_REF].