# BT-ACM+TS+SVM: Enhancing Motor Imagery BCI classification with Block-Toeplitz Augmented Covariance Matrices

Igor Carrara, Theodore Papadopoulos

<figure><img src="image_1.png" /></figure>

## To cite this version:

Igor Carrara, Theodore Papadopoulos. BT-ACM+TS+SVM: Enhancing Motor Imagery BCI classification with Block-Toeplitz Augmented Covariance Matrices. 2024. hal-04797229

HAL Id: hal-04797229

https://inria.hal.science/hal-04797229v1

Preprint submitted on 22 Nov 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# BT-ACM+TS+SVM: ENHANCING MOTOR IMAGERY BCI CLASSIFICATION WITH BLOCK-TOEPLITZ AUGMENTED COVARIANCE MATRICES

Igor Carrara¹ Théodore Papadopoulos¹

¹ Centre Inria Université Côte d'Azur, Equipe Cronos, France
igor.carrara@inria.fr and theodore.papadopoulos@inria.fr

## Abstract

*Objective:* Electroencephalography signals are multidimensional datasets. We introduce an enhancement to the augmented covariance method (ACM), exploiting more thoroughly its mathematical properties, to improve motor imagery classification. *Methods:* Standard ACM emerges as a combination of phase space reconstruction of dynamical systems and of Riemannian geometry. Indeed, it is based on the construction of a Symmetric Positive Definite matrix to improve classification. But this matrix also has a Block-Toeplitz structure that was previously ignored. This property allows a shift in perspective, advocating for an analysis of the ACM within its authentic domain in order to fully use its complex structure and inherent properties. This work treats such matrices in the real manifold to which they belong: the set of Block-Toeplitz SPD matrices. After some manipulation, this set can be seen as the product of an SPD manifold and a Siegel Disk Space. *Results:* The Block-Toeplitz Augmented Covariance Matrix (BT-ACM) approach performs better than any state-of-the-art methods, or at worst obtains statistically compatible performances. We will test our approach on several datasets, several subjects and several tasks using the MOABB framework with a within-session evaluation procedure. *Conclusion:* The proposed approach simplifies the computation by treating each component into its proper geometrical space, bringing improvements not only in terms of performance but also in terms of computational time and carbon emission with respect to the standard ACM methodology, offering a more sustainable and effective solution. *Significance:* This extension represents a meaningful advancement in the field, offering novel insights and methodologies that enrich the existing body of knowledge surrounding Riemannian geometry's role in data classification.

**Keywords** Augmented Covariance, Block-Toeplitz, Brain-computer interface, Riemannian geometry, Takens's Theorem.

## 1 Introduction

Brain-computer interfaces (BCI) can be defined as devices that enable direct communication between humans and machines through the interpretation of brain signals. The primary goal of such technology is to restore communication capabilities and even rehabilitate motor functions, such as walking, in individuals with disabilities, including those suffering from spinal cord injuries [1].

One possible mode of signal acquisition in the field of BCI is through the use of Electroencephalography (EEG). This technology captures the electrical signals produced by neuronal activity within the brain through the placement of electrodes on the scalp i.e., using a noninvasive modality. EEG possesses a high temporal resolution and is relatively inexpensive compared to other neuroimaging modalities. During this research, we will then focus on a particular
paradigm, namely Motor Imagery (MI), where a subject alters their brain activity by imagining themselves performing a movement without physically executing the action.

With a paper published in 2010 [2], the Riemann distance-based classification algorithm established itself as the state-of-the-art for the EEG-BCI community. The core innovation of this method is to treat the spatial covariance, extracted from the EEG signal, not through the conventional Euclidean spaces but using the inherent structure as Symmetric Positive Definite (SPD) matrices, which is actually a differentiable manifold that can be equipped with a Riemannian metric [3].

Recently, an extension of this work has been proposed through the introduction of an Augmented Covariance Matrix (ACM) [4]. This method combines the concept of differential geometry with the phase space reconstruction procedure. By augmenting the spatial covariance matrix with additional dimensions derived from the temporal dynamics of the EEG signal, this method provides a more comprehensive representation of the brain's electrical activity. The ACM matrix generated by this procedure turns out to be an SPD matrix, ensuring that it can be effectively analyzed within the Riemannian manifold framework. However, the ACM methodology also introduces a limitation related to the potential for rank deficiency problems, that can arise when the phase space reconstruction process introduces excessive redundancy. Despite this challenge, the ACM approach offers a promising avenue for enhancing BCI systems' accuracy and reliability.

The current focus is to extend this methodology through the use of all the properties of the ACM matrix. This matrix possess a Block-Toeplitz structure, enriching the attributes of the matrix beyond its SPD properties. Recently, some algorithms and procedures have been proposed to deal with Block-Toeplitz SPD matrix [5], with applications in diverse fields such as audio processing and radar signal analysis [6]. The Block-Toeplitz structure of the ACM matrix, characterized by its repetitive blocks descending diagonally, lends itself to efficient computational techniques and algorithms that can exploit this pattern for enhanced signal processing and classification tasks.

The idea of this research is thus to endow the smooth manifold of Block-Toeplitz SPD matrices with a Riemannian metric, thus allowing the ACM matrix to be treated within its true manifold membership. This approach provides a new more specific metric to use for BCI classification algorithms. The strength of the approach lies in its ability to deconstruct the manifold into its constituent elements: a Symmetric Positive Definite (SPD) manifold and several Siegel Disk Spaces. This process is possible thanks to an appropriate conversion of the blocks of the ACM matrix into the Verblusky coefficients [7]. By discerningly analyzing each component within its respective geometrical domain, this method significantly alleviates the computational demand traditionally associated with the ACM methodology.

To the best of the author's knowledge, such an approach is a complete novelty regarding the applications to MI BCI-EEG classification. Furthermore, previous algorithms in the radar field primarily operated directly on the manifold. In contrast, this research proposes the strategy of mapping to the Tangent Space, an approach that has been identified for its superior discriminatory capabilities [8].

In order to validate the methods, we will test our approach against thirteen pipelines of the state-of-the-art (both Machine Learning (ML) and Deep Learning (DL)), on several datasets, several subjects and several tasks using the MOABB framework [9], with a within-session evaluation procedure. The decision to use MOABB places a strong emphasis on ensuring that our method is both reproducible and stable across various datasets.

The article is structured as follows: in section 2, we describe the theoretical framework of the methodology presented, the considered dataset and the evaluation procedure used. In section 3, we list the results obtained using within-session for the task right hand versus left hand and right hand versus feet. Finally, in section 4, we describe the findings of this research and the limitation. Section 5 summarizes the results of our study.
# 2 Material and Methods

## 2.1 Augmented Covariance Matrix for MI-BCI classification

The EEG signal can be represented as a multivariate time series, denoted as $\mathbf{X} \in \mathbb{R}^{d \times T}$, where $T$ represents the total number of sampled data points, and $d$ indicates the number of electrodes used in the EEG recording.

Since this paper focuses on MI task, we split the EEG signals into smaller sections known as epochs, each representing a snapshot of brain activity during various tasks or mental states. The objective is thus to develop a method that can accurately identify the specific task or mental state associated with these EEG segments.

Over the past decade, the field of EEG-BCI has been revolutionized by the Riemannian approaches applied for the classification of the SPD matrix classification [2]. At the heart of this approach lies the idea of mapping the covariance matrix onto its proper geometrical space: the smooth manifold of SPD matrices. On this manifold is possible to endow several Riemannian metrics, finally converting the SPD matrix set to a Riemannian manifold [3].

The spaces of SPD matrices are composed by square real symmetric matrices that are positive definite, formally the space of $\mathbb{SPD}_d$ matrices is defined as

$$ \mathbb{SPD}_d = \{ \mathbf{M} \in \mathbb{R}^{d \times d} | x^T \mathbf{M} x > 0 \ \forall x \in \mathbb{R}^d \setminus \{0\} \} \quad (1) $$

In this way, we can represent a matrices $\mathbf{M} \in \mathbb{SPD}_d$ as a point on a Riemann manifold of dimension $d(d+1)/2$.

The ACM [4] that extends the standard Riemannian methodology by combining it with the phase space reconstruction (PSR) approach that is grounded in the Takens theorem [10]. ACM combines spatial and temporal information from the signal, allowing an enrichment of the information extracted from the signal and improving classification performance.

The ACM methodology can be derived using the connection with the Takens delay-embedding theorem [10] in the context of nonlinear dynamical systems. In particular, we use the formulation of Takens theorem for stochastic systems [11]. The idea of using Takens theorem is based on the idea that time series obtained from experimental observations that capture only a fraction of the complex dynamics of the underlying system, can nonetheless be used to reconstruct the system's full dynamical behavior. This is achieved using a uniform embedding procedure: consider a time series $s(n)$ created through a measurement process, the PSR technique generates a point $\mathbf{s}_E(n)$ of a $D_E$-dimensional space constructed as

$$ \mathbf{s}_E(n) = [s(n), s(n - \tau), \dots, s(n - (D_E - 1)\tau)]^T \quad (2) $$

where $\tau$ is a positive integer called the embedding delay and $D_E$ is the embedding dimension. $\mathbf{s}_E(n)$ is an uniform embedding of the original phase space. This space, if constructed correctly, effectively captures the geometric structure and dynamics of the original system.

The ACM matrix is obtained by expanding the original EEG signal using the PSR approach with an embedding dimension $p$ to get a new $dp \times (T - p\tau)$ time series, parameterized by the fixed delay $\tau$. The Augmented Covariance Matrix $\Gamma_{aug}$ is defined as the covariance matrix of this new augmented dataset:

$$ \Gamma_{Aug} = \begin{bmatrix} \Gamma_0 & \Gamma_{-1} & \Gamma_{-2} & \cdots \\ \Gamma_1 & \Gamma_0 & \Gamma_{-1} & \cdots \\ \Gamma_2 & \Gamma_1 & \Gamma_0 & \cdots \\ \vdots & \vdots & \vdots & \ddots \\ \Gamma_{D_E-1} & \Gamma_{D_E-2} & \Gamma_{D_E-3} & \cdots \end{bmatrix}, \quad (3) $$

The $\Gamma_{Aug}$ exhibits positive definiteness when noise is present in the signal [6]. Interestingly, if constructed correctly, this matrix has an additional characteristic that has not been fully
explored until now: the matrix possesses a Block-Toeplitz structure, with blocks of dimension *d* by *d*. More formally, Γ<sub>Aug</sub> belongs to the space $BTSPD_{d×p}$ of Block-Toeplitz and SPD matrices, i.e., SPD matrices of size *dp* × *dp* with constant blocks of size *d* × *d* along all diagonals. This opens up new possibilities for enhancing the ACM formulation by mapping the ACM matrix to the most suitable geometric space that fully captures both its Block-Toeplitz and SPD natures.

To make the Block Toeplitz structure of the matrix evident, it was necessary to slightly modify the procedure of *AugmentedDataset*, i.e. the algorithm used to generate the signal expansion $X_{d×p}(t)$ for the ACM algorithm. Originally, the algorithm, upon expanding the signal, would select only the overlapping segments of the signal, inadvertently introducing edge effects that obscured the natural formulation of Γ<sub>Aug</sub> as BT-SPD matrices. The revised strategy, *AugmentedDatasetBT*, involves expanding the signal in the same initial manner but then incorporating null values in areas lacking overlapping signal segments. This approach effectively prevents any distortion effects that could compromise the matrix's BT-SPD structure. The procedure is outlined in Figure 1.

<figure><img src="image_3.png" /><figcaption>Figure 1: The presented example uses the BNCI2014004 dataset that has only 3 electrodes. In this figure, we consider an embedding corresponding to $D_E = 3$ and $\tau = 10$. (a) *Augmented-Dataset* is created, the final covariance estimator, representing the $\Gamma_{aug}$, computed over the overlapping part of the signal (basically discarding the blue component). In this way, each row contains a slightly different signal. (b) *AugmentedDatasetBT* is created, the final covariance estimator, representing the $\Gamma_{aug}$, is computed over all the signal and the uncommon part (blue zone) are filled using zero value. This procedure specifically guarantees that the final covariance has a Block Toeplitz structure. This procedure allows the information contained in each row to be the same up to the shift.</figcaption></figure>

## 2.2 Mapping ACM to Verblusky coefficients

While $\Gamma_0$ is an SPD matrix, the other fundamental block ($\Gamma_1, ..., \Gamma_{p-1}$) must also adhere to specific properties since the complete matrix $\Gamma_{Aug}$ belongs to $BTSPD_{d×p}$. Nevertheless, such properties are not obvious at first look. However, by implementing a minor adjustment to the coefficients, it is possible to ensure their belonging within the domain of the Siegel disk [7, 12], defined as

$$ \Omega_i \in \mathbb{SD}_d \quad \mathbb{SD}_d = \{\mathbf{M} \in \mathbb{C}^{d \times d}; \mathbf{I} - \mathbf{MM}^H > 0\} \quad (4) $$
where $^H$ is the hermitian operator. We use the partial ordering of the set of complex matrices: we note $A > B$ when the difference $A - B$ is a positive definite matrix. The coefficients that have undergone such modification are also known as Verblunsky coefficients.

We now present the transformation [7] that allows us to map the fundamental block matrix $\Gamma_{Aug}$ to the Verblunsky coefficients $\Omega_i$. Basically, it is possible to define an isomorphism

$$ \begin{gathered} \mathbb{B}_{d \times p} \rightarrow \mathbb{SPD}_d \times \mathbb{SD}_d^{p-1} \\ \Gamma_{Aug} \mapsto (\Gamma_0, \Omega_1, \dots, \Omega_{p-1}). \end{gathered} \qquad (5) $$

This transformation is described through the Algorithm 1. Consider the initial matrix $\Gamma_{Aug}$ decomposed in its constituent blocks $(\Gamma_0, \dots, \Gamma_{p-1})$. The initialization of the recursive transformation is set to $\mathbf{P}_0 = \mathbf{\Gamma}_0$. The subsequent coefficients are composed recursively using the Algorithm 1. Note that in the algorithm $\tilde{\Gamma}_{l-1}$ identifies the sub-matrix of $\Gamma_{Aug}$ confined to the first $l-1$ row and column.

**Algorithm 1 ACM conversion to Verblunsky coefficients**

Input ACM matrix $\Gamma_{Aug} \in \mathbb{BTSPD}_{d \times p}$

Initialization

$$ P_0 = \Gamma_0 $$

for $l = 0 : p-1$ do

$$ \Omega_{l+1} = \mathbf{L}_l^{-1/2}(\Gamma_{l+1} - \mathbf{M}_l)\mathbf{K}_l^{-1/2} $$

using

$$ \mathbf{L}_l = \mathbf{P}_0 - (\Gamma_1, \dots, \Gamma_l)\tilde{\Gamma}_{l-1}^{-1}(\Gamma_1, \dots, \Gamma_l)^H $$

$$ \mathbf{K}_l = \mathbf{P}_0 - (\Gamma_1^H, \dots, \Gamma_l^H)\tilde{\Gamma}_{l-1}^{-1}(\Gamma_1^H, \dots, \Gamma_l^H)^H $$

$$ \mathbf{M}_l = (\Gamma_1, \dots, \Gamma_l)\tilde{\Gamma}_{l-1}^{-1}(\Gamma_1^H, \dots, \Gamma_l^H)^H $$

end for

$$ \textbf{Return} (\Gamma_0, \Omega_1, \dots, \Omega_{p-1}) \in \mathbb{SPD}_d \times \mathbb{SD}_d^{p-1} $$

The homeomorphism introduced by the Verblunsky transformation clarifies the space to which the off-diagonal blocks belong, revealing that they are constrained to lie within the Siegel Disk. This transformation thus plays a crucial role in making the essential characteristics of these blocks clear.

## 2.3 Information Theory for Block Toeplitz SPD matrices

We now proceed to formalize the Riemannian metric for matrices $\Gamma_{Aug} \in \mathbb{B}_{d \times p}$. For an in-depth exploration of this formulation, readers are encouraged to refer to [6, 5]. The set of Block-Toeplitz SPD $\mathbb{BTSPD}_{d \times p}$ is considered as a Kähler manifold [5], on which it is possible to define a Kähler potential $\Phi$, computed as the entropy of the signal [13], as follows

$$ \Phi(\Gamma_{Aug}) = -\log(\det(\Gamma_{Aug})) - \log(\pi e) \qquad (6) $$

This function is also known as the Kähler potential [14, 15]. After applying some decomposition properties of the determinant of $\Gamma_{Aug}$ we can rewrite the potential as

$$ \begin{aligned} \Phi(\Gamma_{Aug}) = & -n \log(\det \Gamma_0) - \sum_{l=1}^{n-1} (n-l) \log(\det(\mathbf{I} - \Omega_1 \Omega_1^H)) \\ & - \log(\pi e) \end{aligned} \qquad (7) $$
where $\Gamma_{Aug}$ is identified with coefficient $(\Gamma_0, ..., \Gamma_{p-1})$. Finally, it is possible to define the metric of the considered matrices by computing the Hessian matrix of the Kähler potential.

$$ds^2 = p \operatorname{trace} (\mathbf{P}_0^{-1} d\mathbf{P}_0 \mathbf{P}_0^{-1} d\mathbf{P}_0) + \sum_{l=1}^{p-1} (p-l) \operatorname{trace} \left( (\mathbf{I} - \Omega_l \Omega_l^H)^{-1} d\Omega_l (\mathbf{I} - \Omega_l^H \Omega_l)^{-1} d\Omega_l^H \right) \quad (8)$$

The first term of Equation (8) is identifiable as the metric for the SPD<sub>d</sub> space. The other terms represent the metrics for the Siegel disk spaces SD<sub>d</sub> repeated *p* − 1 times, i.e. the space of Block-Toeplitz SPD matrices is equipped with a product Riemannian metric over SPD<sub>d</sub> × SD<sub>d</sub><sup>*p*−1</sup>.

Interestingly, since the final space is a product of Riemannian manifold it is possible to treat each component independently, using the unique geometric and analytical tools pertinent to its own manifold. This consideration will be crucial in constructing the proposed methodology. Let us now recall some of the properties of the metric spaces under discussion, a structured overview of the relevant details can be viewed in Table 1.

<table><caption>Table 1: The table provides a comprehensive overview of the fundamental concepts of differential geometry for the considered spaces SPD<sub>d</sub> and SD<sub>d</sub>. In the table we use <sup>H</sup> is the conjugate transpose operation, <sup>-1</sup> as the inverse operator.</caption><thead><tr><th></th><th>SPD<sub>d</sub></th><th>SD<sub>d</sub></th></tr></thead><tbody><tr><td>Definition</td><td>SPD<sub>d</sub> = {M ∈ ℝ<sup>d×d</sup>; x<sup>T</sup>Mx > 0, ∀x ∈ ℝ<sup>d</sup>\{0}</td><td>SD<sub>d</sub> = {M ∈ ℝ<sup>d×d</sup>; I - MM<sup>H</sup> > 0}</td></tr><tr><td>Metric</td><td>ds<sup>2</sup><sub>SPD<sub>d</sub></sub> = trace (P<sup>-1</sup>dPP<sup>-1</sup>dP)</td><td>ds<sup>2</sup><sub>SD<sub>d</sub></sub> = trace ((I - Ω<sub>l</sub>Ω<sub>l</sub><sup>H</sup>)<sup>-1</sup> dΩ<sub>l</sub> (I - Ω<sub>l</sub><sup>H</sup>Ω<sub>l</sub>)<sup>-1</sup> dΩ<sub>l</sub><sup>H</sup>)</td></tr><tr><td>Logarithmic Map</td><td>log<sub>A</sub>(B) = A<sup>1/2</sup> log(A<sup>-1/2</sup>BA<sup>-1/2</sup>)A<sup>1/2</sup></td><td>log<sub>Ω</sub>(ψ) = (I - Ω<sub>l</sub>Ω<sub>l</sub><sup>H</sup>)<sup>1/2</sup>V'(I - Ω<sub>l</sub><sup>H</sup>Ω)<sup>1/2</sup></td></tr><tr><td>Exponential Map</td><td>exp<sub>A</sub>(V) = A<sup>1/2</sup> exp(A<sup>-1/2</sup>V A<sup>-1/2</sup>)A<sup>1/2</sup></td><td>exp<sub>Ω</sub>(V) = (I - Ω<sub>l</sub>Ω<sub>l</sub><sup>H</sup>)<sup>-1/2</sup>(ψ' + Ω)(I + Ω<sub>l</sub><sup>H</sup>ψ')<sup>-1</sup>(I - Ω<sub>l</sub><sup>H</sup>Ω)<sup>1/2</sup></td></tr><tr><td>Curvature</td><td>-1 ≤ K(σ) ≤ 0</td><td>-4 ≤ K(σ) ≤ 0</td></tr></tbody></table>

## 2.4 Proposed Algorithm: BT-ACM+TS+SVM

The proposed algorithm, Block Toeplitz Augmented Covariance Matrix (BT-ACM) with Tangent Space projection (TS) and SVM classifier, combines all the mentioned elements together. A visual representation of the methodology employed in this chapter is available in Figure 2.

① **Phase Space Reconstruction:** The core concept involves expanding the signal through the application of nonlinear theory, using the approach outlined for the Augmented-DatasetBT. Such expansion is controlled by the hyper-parameters $D_E$ and $\tau$, carefully chosen using a grid-search procedure. Note that for $D_E = \tau = 1$, we obtain, as a special case, the state-of-the-art [2] approach.

② **BT-ACM:** On the expanded signal the covariance is computed, using a regularization through the Oracle Approximating Shrinkage Estimator (OAS) [16] in order to ensure that the resulting matrices are always positive definite. The obtained matrix is $\Gamma_{Aug} \in \text{BTSPD}_{d \times p}$.
<figure><img src="image_4.png" /><figcaption>Figure 2: Illustration of the BL-ACM + TS + SVM 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 thus a 3-dimensional time series. The method then begins with the extraction of the epochs representing left and right-hand tasks. The selection of hyper-parameters can be made via grid search using the nested approach. Once these parameters are selected, we use the phase space reconstruction process to obtain a dynamic system equivalent to the original one hidden in the brain (In this figure, we see an embedding corresponding to $D_E = 3$ and $\tau = 10$). From this high-dimensional space, the BL-ACM matrix is extracted, and the main blocks are converted into Verblunsky coefficients. After this step, each component is mapped through the tangent space using the appropriate transformation depending on the Riemannian manifold and subsequently vectorized. The final step is the application of a classification algorithm based on SVM.</figcaption></figure>

③ **ACM to Verblunsky coefficients:** the next step is then to extract the $p$ fundamental blocks of the matrix, each of dimension $d \times d$. Following this extraction, we apply the transformation detailed in Algorithm 1, which effectively converts the matrix into Verblunsky coefficients.

$$ \Gamma_{Aug} \mapsto (\Gamma_0, \Omega_1, \dots, \Omega_{p-1}) \qquad (9) $$

④ **Tangent Space Projection:** The final space is a product of manifold endowed with a product metric. It is therefore possible to estimate independently each component. So consider having two matrices $\Gamma_{Aug1}$ and $\Gamma_{Aug2}$ both belonging to $\text{BTSPD}_{d \times p}$, we can compute their mean by simply computing the mean of each component in their respectively manifold

$$ \begin{gathered} \Gamma_{Aug1} \mapsto (\Gamma_{0,1}, \Omega_{1,1}, \dots, \Omega_{p-1,1}) \\ \Gamma_{Aug2} \mapsto (\Gamma_{0,2}, \Omega_{1,2}, \dots, \Omega_{p-1,2}) \\ \downarrow \quad \downarrow \quad \downarrow \\ \mathcal{G} \leftarrow (\mathcal{G}_{\Gamma_0}, \mathcal{G}_{\Omega_1}, \dots, \mathcal{G}_{\Omega_{p-1}}) \end{gathered} \qquad (10) $$

For the proposed algorithm, we focus on a Tangent Space methodology since it has been shown to obtain superior performance when compared to Riemannian surface classification [8]. This approach is based on projecting each component onto the tangent space, using the appropriate transformation that corresponds to the manifold of each component, using exactly the same methodology described in Equation 10.
⑤ **Vectorization:** Following the projection onto the Tangent Space, we vectorize the matrices. A common approach for the part concerning $SPD_d$ is to use the methodology of minimal representation [17] since, in this case, we are dealing with symmetric matrices, the number of unique parameters to be considered reduces to $d(d+1)/2$. In this study, we consciously chose to avoid that procedure in order to maintain uniformity across each block, irrespective of whether it is associated with $SPD_d$ or with $SD_d$, facilitating a more straightforward integration of blocks from both spaces within our analysis. Actually, for the part concerning $SD_d$, such minimal representation is not applicable since the number of independent parameters are $d^2$ for each matrix, so we use standard vectorization.

⑥ **Classification:** After the vectorization procedure, we combine the vectors from the different components of the space $SPD_d \times SD_d^{p-1}$ into a single vector, which is then classified with an SVM algorithm.

## 2.5 Dataset

In order to evaluate the effectiveness of our proposed methodology, we use publicly accessible datasets provided by the MOABB [9] framework. This choice was driven by a commitment to conduct open and reproducible research. We selected nine motor imagery datasets, consisting of over 250 subjects. Detailed information regarding these datasets, including their specific characteristics and the number of subjects and sessions, is summarized in Table 2.

<table><caption>Table 2: Dataset considered during this study.</caption><thead><tr><th>Dataset</th><th>subjects</th><th>channels</th><th>sampling rate</th><th>sessions</th><th>tasks</th><th>trials/class</th><th>Epoch (s)</th></tr></thead><tbody><tr><td>BNCI2014001 [18]</td><td>9</td><td>22</td><td>250 Hz</td><td>2</td><td>4</td><td>144</td><td>[2, 6]</td></tr><tr><td>BNCI2014002 [19]</td><td>14</td><td>15</td><td>512 Hz</td><td>1</td><td>2</td><td>80</td><td>[3, 8]</td></tr><tr><td>BNCI2014004 [20]</td><td>9</td><td>3</td><td>250 Hz</td><td>5</td><td>2</td><td>360</td><td>[3, 7.5]</td></tr><tr><td>BNCI2015001 [21]</td><td>12</td><td>13</td><td>512 Hz</td><td>3</td><td>2</td><td>200</td><td>[0, 5]</td></tr><tr><td>Cho2017 [22]</td><td>52</td><td>64</td><td>512 Hz</td><td>1</td><td>2</td><td>100</td><td>[0, 3]</td></tr><tr><td>Lee2019-MI [23]</td><td>54</td><td>62</td><td>1000 Hz</td><td>2</td><td>2</td><td>100</td><td>[0, 4]</td></tr><tr><td>PhysionetMI [24]</td><td>109</td><td>64</td><td>160 Hz</td><td>1</td><td>4</td><td>23</td><td>[0, 3]</td></tr><tr><td>Weibo2014 [25]</td><td>10</td><td>60</td><td>200 Hz</td><td>1</td><td>7</td><td>80</td><td>[3, 7]</td></tr><tr><td>Zhou2016 [26]</td><td>4</td><td>14</td><td>250 Hz</td><td>3</td><td>3</td><td>160</td><td>[0, 5]</td></tr></tbody></table>

On the considered datasets, we use established preprocessing protocols tailored for the MI task, incorporating a standard band-pass filter in the 8 to 32 Hz frequency range [27]. For the pipelines based on the Filter Bank approach, we divided the EEG signal into six non-overlapping frequency bands: 8–12 Hz, 12–16 Hz, 16–20 Hz, 20–24 Hz, 24–28 Hz, and 28–32 Hz. The duration of each epoch within our study is intentionally aligned with the task conditions’ length, which is subject to variation across the datasets employed. In the datasets considered, all classes are balanced.

## 2.6 Evaluation Procedure

This work is based on a Within-Session (WS) evaluation procedure to assess the performance of our methodologies. This procedure focuses on a singular session for each participant. The MOABB framework employs a Nested Cross-Validation methodology [28], structured with an outer loop of 5-Fold Cross-validation and an inner one composed by a 3-fold Cross Validation, as described in Figure 3. The metrics used to evaluate the performance will always be ROC-AUC since we are dealing only with a binary classification problem. To analyze the ranking of the different pipelines we used the statistical framework provided by MOABB.
<figure><img src="image_5.png" /><figcaption>Figure 3: Nested cross-validation procedure. The outer loop uses a stratified 5-fold CV for splitting the training and testing data, while the inner loop applied on the training dataset of each fold uses a stratified 3-fold CV to tune the model's hyper-parameters.</figcaption></figure>

## 2.7 Studied Pipelines

In order to evaluate the proposed methodology we test it in comparison with several pipelines from the state-of-the-art of the BCI literature. In order to provide a complete analysis of the performance of our pipeline we use 13 pipelines from both Machine Learning (ML) and Deep Learning (DL). All the details of the pipeline considered are listed in Table 3. 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 re-sampling 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, refer to Table A2.

<table><caption>Table 3: Pipelines considered in this study are organized into two distinct sections within the table: the first part is dedicated to the traditional classical ML pipelines while the second part focuses on the DL pipelines for MI.</caption><thead><tr><th>Pipeline</th><th>Feature Extraction</th><th>Classifier</th></tr></thead><tbody><tr><td>CSP + LDA [29]</td><td>Common Spatial Patterns (CSP) with OAS covariance estimator</td><td>Optimized Shrinkage LDA</td></tr><tr><td>FBCSP+SVM [30]</td><td>Filter Bank Common Spatial Patterns (FBCSP) using OAS</td><td>Optimized Support Vector Machine (SVM)</td></tr><tr><td>MDM [2]</td><td>Spatial Covariance using OAS</td><td>Mean Distance to Mean (MDM)</td></tr><tr><td>FgMDM [2]</td><td>Spatial Covariance using OAS</td><td>Minimum Distance to Mean with geodesic filtering (FgMDM)</td></tr><tr><td>TS + EN [31]</td><td>Spatial Covariance using OAS mapped to TS</td><td>Optimized Elastic Network (EL)</td></tr><tr><td>TS + SVM [2]</td><td>Spatial Covariance using OAS mapped to TS</td><td>Optimized SVM</td></tr><tr><td>ACM + TS + SVM [4]</td><td>ACM with Sample Covariance Estimator mapped to TS</td><td>Optimized SVM</td></tr><tr><td>ShallowConvNet [32]</td><td>Standardized and resample EEG signal at 250Hz</td><td>Convolutional Neural Network (CNN)</td></tr><tr><td>DeepConvNet [32]</td><td>Standardized and resample EEG signal at 250Hz</td><td>CNN</td></tr><tr><td>EEGNet 8 2 [33]</td><td>Standardized and resample EEG signal at 128Hz</td><td>CNN with architecture EEGNet</td></tr><tr><td>EEGTCNet [34]</td><td>Standardized and resample EEG signal at 250Hz</td><td>CNN with architecture EEGTCNet</td></tr><tr><td>EEGITNet [35]</td><td>Standardized and resample EEG signal at 128Hz</td><td>CNN with architecture EEGITNet</td></tr><tr><td>EEGNeX 8 32 [36]</td><td>Standardized and resample EEG signal at 128Hz</td><td>CNN with architecture EEGNeX</td></tr></tbody></table>

When possible, the hyper-parameters of the classification models were optimized using a Grid Search procedure. Table A1 summarizes the domains over which these grid searches are performed.

# 3 Results

In this section, we discuss the results of our proposed method and compare its performance with the current state-of-the-art in the field.

## 3.1 Right Hand vs Left Hand

Regarding the classification task of Right vs Left Hand, we use 7 datasets: BNCI2014001, BNCI2014004, Cho2017, Lee2019-MI, PhysionetMI, Weibo2014 and Zhou2016.
<table><caption>Table 4: Performance (AUC) of Right hand vs Left hand classification.</caption><thead><tr><th>Models</th><th>BNCI2014001</th><th>BNCI2014004</th><th>Cho2017</th><th>Lee2019-MI</th><th>PhysionetMI</th><th>Weibo2014</th><th>Zhou2016</th></tr></thead><tbody><tr><td>BT-ACM+TS+SVM (Our)</td><td>0.89 ± 0.11</td><td><strong>0.83 ± 0.14</strong></td><td><strong>0.76 ± 0.14</strong></td><td>0.84 ± 0.13</td><td>0.67 ± 0.19</td><td>0.84 ± 0.15</td><td>0.95 ± 0.06</td></tr><tr><td>CSP+LDA</td><td>0.82 ± 0.17</td><td>0.80 ± 0.15</td><td>0.71 ± 0.15</td><td>0.77 ± 0.17</td><td>0.66 ± 0.17</td><td>0.81 ± 0.13</td><td>0.93 ± 0.07</td></tr><tr><td>FBCSP+SVM</td><td>0.84 ± 0.16</td><td>0.80 ± 0.16</td><td>0.68 ± 0.16</td><td>0.75 ± 0.17</td><td>0.59 ± 0.14</td><td>0.77 ± 0.19</td><td>0.94 ± 0.05</td></tr><tr><td>MDM</td><td>0.82 ± 0.15</td><td>0.78 ± 0.16</td><td>0.63 ± 0.14</td><td>0.70 ± 0.14</td><td>0.55 ± 0.17</td><td>0.59 ± 0.16</td><td>0.91 ± 0.07</td></tr><tr><td>FgMDM</td><td>0.87 ± 0.12</td><td>0.79 ± 0.15</td><td>0.73 ± 0.13</td><td>0.81 ± 0.14</td><td><strong>0.69 ± 0.19</strong></td><td>0.78 ± 0.15</td><td>0.93 ± 0.07</td></tr><tr><td>TS+EL</td><td>0.86 ± 0.13</td><td>0.80 ± 0.15</td><td><strong>0.76 ± 0.14</strong></td><td><strong>0.85 ± 0.13</strong></td><td>0.68 ± 0.20</td><td><strong>0.85 ± 0.12</strong></td><td>0.94 ± 0.06</td></tr><tr><td>TS+SVM</td><td>0.87 ± 0.14</td><td>0.79 ± 0.15</td><td>0.75 ± 0.14</td><td>0.84 ± 0.14</td><td>0.68 ± 0.20</td><td>0.84 ± 0.14</td><td>0.93 ± 0.06</td></tr><tr><td>ACM+TS+SVM</td><td><strong>0.92 ± 0.10</strong></td><td><strong>0.83 ± 0.15</strong></td><td>0.74 ± 0.15</td><td>0.83 ± 0.14</td><td>0.64 ± 0.21</td><td><strong>0.85 ± 0.13</strong></td><td>0.95 ± 0.05</td></tr><tr><td>ShallowConvNet</td><td>0.86 ± 0.14</td><td>0.72 ± 0.18</td><td>0.74 ± 0.15</td><td>0.76 ± 0.15</td><td>0.65 ± 0.16</td><td>0.79 ± 0.13</td><td><strong>0.96 ± 0.06</strong></td></tr><tr><td>DeepConvNet</td><td>0.82 ± 0.16</td><td>0.72 ± 0.19</td><td>0.72 ± 0.13</td><td>0.71 ± 0.16</td><td>0.60 ± 0.17</td><td>0.74 ± 0.16</td><td>0.94 ± 0.06</td></tr><tr><td>EEGNet</td><td>0.77 ± 0.19</td><td>0.70 ± 0.20</td><td>0.67 ± 0.16</td><td>0.66 ± 0.16</td><td>0.60 ± 0.16</td><td>0.67 ± 0.22</td><td>0.95 ± 0.03</td></tr><tr><td>EEGTCNet</td><td>0.68 ± 0.21</td><td>0.70 ± 0.20</td><td>0.58 ± 0.13</td><td>0.57 ± 0.13</td><td>0.56 ± 0.13</td><td>0.63 ± 0.18</td><td>0.82 ± 0.09</td></tr><tr><td>EEGITNet</td><td>0.75 ± 0.16</td><td>0.65 ± 0.15</td><td>0.57 ± 0.12</td><td>0.59 ± 0.12</td><td>0.53 ± 0.11</td><td>0.59 ± 0.14</td><td>0.69 ± 0.15</td></tr><tr><td>EEGNeX 8 32</td><td>0.66 ± 0.13</td><td>0.67 ± 0.17</td><td>0.53 ± 0.11</td><td>0.55 ± 0.10</td><td>0.51 ± 0.11</td><td>0.58 ± 0.16</td><td>0.62 ± 0.15</td></tr></tbody></table>

The various approaches of this study are compared in Table 4. A detailed picture of the statistical significance of these results is provided in Figure 4.

Overall, our method scores the best in 2 datasets - BNCI2014001 and Cho2017 - and obtains similar performance with respect to ACM+TS+SVM and TS+EL respectively. Across other evaluated datasets, our approach consistently delivers results that closely rival those of other leading algorithms, with a marginal performance deviation of no more than 1% in the AUC score with the only exception of the BNCI2014001 dataset. For instance, while the Zhou2016 dataset reveals ShallowConvNet as the top performer, our method, BT-ACM+TS+SVM, follows it by a thin margin, demonstrating robust effectiveness.

Moreover, a comprehensive analysis across multiple datasets underscores a statistically significant performance enhancement achieved by BT-ACM+TS+SVM compared to all considered algorithms, with the only exception of TS+EL, where the outcomes are remarkably similar (Figure 4 (a) and (d)). This evidence collectively affirms the superiority of our method, not only in achieving high-performance benchmarks but also in maintaining competitive results across datasets.

## 3.2 Right Hand vs Feet

Regarding the classification task of Right Hand vs Feet, we use 7 datasets: AlexMI, BNCI2014001, BNCI2014002, BNCI2015001, PhysionetMI, Weibo2014 and Zhou2016.

<table><caption>Table 5: Performance (AUC) of Right hand vs Feet classification.</caption><thead><tr><th>Models</th><th>AlexMI</th><th>BNCI2014001</th><th>BNCI2014002</th><th>BNCI2015001</th><th>PhysionetMI</th><th>Weibo2014</th><th>Zhou2016</th></tr></thead><tbody><tr><td>BT-ACM+TS+SVM (Our)</td><td>0.84 ± 0.17</td><td>0.96 ± 0.05</td><td>0.87 ± 0.12</td><td>0.92 ± 0.09</td><td><strong>0.94 ± 0.07</strong></td><td>0.92 ± 0.05</td><td><strong>0.97 ± 0.02</strong></td></tr><tr><td>CSP+LDA</td><td>0.77 ± 0.18</td><td>0.92 ± 0.10</td><td>0.81 ± 0.15</td><td>0.89 ± 0.11</td><td>0.86 ± 0.14</td><td>0.89 ± 0.06</td><td>0.95 ± 0.03</td></tr><tr><td>FBCSP+SVM</td><td>0.84 ± 0.16</td><td>0.80 ± 0.16</td><td>0.68 ± 0.16</td><td>0.75 ± 0.17</td><td>0.59 ± 0.14</td><td>0.77 ± 0.19</td><td>0.94 ± 0.05</td></tr><tr><td>MDM</td><td>0.74 ± 0.21</td><td>0.89 ± 0.10</td><td>0.78 ± 0.14</td><td>0.86 ± 0.13</td><td>0.82 ± 0.12</td><td>0.65 ± 0.10</td><td>0.92 ± 0.04</td></tr><tr><td>FgMDM</td><td>0.80 ± 0.18</td><td>0.94 ± 0.08</td><td>0.85 ± 0.11</td><td>0.90 ± 0.10</td><td>0.90 ± 0.11</td><td>0.89 ± 0.05</td><td>0.96 ± 0.03</td></tr><tr><td>TS+EL</td><td>0.81 ± 0.21</td><td>0.95 ± 0.07</td><td>0.86 ± 0.11</td><td>0.91 ± 0.09</td><td>0.90 ± 0.07</td><td>0.92 ± 0.04</td><td><strong>0.97 ± 0.03</strong></td></tr><tr><td>TS+SVM</td><td>0.79 ± 0.20</td><td>0.94 ± 0.08</td><td>0.86 ± 0.12</td><td>0.91 ± 0.09</td><td><strong>0.94 ± 0.07</strong></td><td>0.92 ± 0.05</td><td>0.96 ± 0.03</td></tr><tr><td>ACM+TS+SVM</td><td><strong>0.87 ± 0.11</strong></td><td><strong>0.97 ± 0.03</strong></td><td>0.87 ± 0.11</td><td><strong>0.93 ± 0.08</strong></td><td>0.93 ± 0.09</td><td><strong>0.93 ± 0.04</strong></td><td><strong>0.97 ± 0.03</strong></td></tr><tr><td>ShallowConvNet</td><td>0.64 ± 0.18</td><td>0.93 ± 0.08</td><td><strong>0.88 ± 0.12</strong></td><td>0.91 ± 0.11</td><td>0.75 ± 0.15</td><td>0.89 ± 0.06</td><td><strong>0.97 ± 0.02</strong></td></tr><tr><td>DeepConvNet</td><td>0.62 ± 0.19</td><td>0.88 ± 0.12</td><td><strong>0.88 ± 0.11</strong></td><td>0.88 ± 0.13</td><td>0.72 ± 0.16</td><td>0.79 ± 0.13</td><td>0.96 ± 0.04</td></tr><tr><td>EEGNet</td><td>0.64 ± 0.16</td><td>0.89 ± 0.15</td><td>0.84 ± 0.16</td><td>0.90 ± 0.12</td><td>0.74 ± 0.16</td><td>0.78 ± 0.15</td><td>0.95 ± 0.03</td></tr><tr><td>EEGTCNet</td><td>0.61 ± 0.22</td><td>0.75 ± 0.19</td><td>0.74 ± 0.19</td><td>0.77 ± 0.19</td><td>0.57 ± 0.13</td><td>0.62 ± 0.12</td><td>0.86 ± 0.16</td></tr><tr><td>EEGITNet</td><td>0.48 ± 0.10</td><td>0.76 ± 0.13</td><td>0.71 ± 0.18</td><td>0.72 ± 0.17</td><td>0.55 ± 0.12</td><td>0.63 ± 0.12</td><td>0.80 ± 0.17</td></tr><tr><td>EEGNeX 8 32</td><td>0.52 ± 0.15</td><td>0.66 ± 0.14</td><td>0.70 ± 0.20</td><td>0.72 ± 0.20</td><td>0.52 ± 0.12</td><td>0.60 ± 0.12</td><td>0.65 ± 0.17</td></tr></tbody></table>

The various approaches of this study are compared in Table 5. A detailed picture of the statistical significance of these results is provided in Figure A1.

Overall, our method scores the best in 2 datasets - PhysionetMI and Zhou2016 - and obtains similar performance when compared to other algorithms from the state-of-the-art. Across the
other evaluated datasets, our approach consistently delivers results that closely rival those of other leading algorithms, with a marginal performance deviation of no more than 1% in the AUC score with the only exception of the AlexMI dataset.

Moreover, a comprehensive analysis across multiple datasets underscores a statistically significant performance enhancement achieved by BT-ACM+TS+SVM compared to all considered algorithms, with the sole exception of TS+EL, TS+SVM and ACM+TS+SVM, where the statistical analysis shows that they have no significant difference (Figure A1 (a), (c), (d) and (e)). This evidence collectively affirms the superiority of our method, setting it among the top-performing algorithms in the MI BCI-EEG domain.

# 4 Discussion

## 4.1 BT-ACM+TS+SVM approach is at worst comparable with the state-of-the-art

Our analysis across multiple datasets demonstrates that the BT-ACM+TS+SVM algorithm not only competes with but frequently surpasses the performance of current state-of-the-art methodologies. Specifically, for the task of Right Hand vs. Left Hand classification, our proposed algorithm consistently outperforms existing benchmarks, establishing a new standard in classification performance. The only exception to this trend is when compared to the TS+EL algorithm, where our results are statistically indistinguishable.

Regarding the classification Right Hand vs. Feet (a simpler classification than the previous one according to [8]), the BT-ACM+TS+SVM algorithm continues to exhibit exceptional performance. In this case, the algorithm obtains statistically indistinguishable results from TS+EL, TS+SVM and ACM+TS+SVM.

This consistent high-level performance across diverse datasets underlines the BT-ACM+TS+SVM algorithm's capacity to maintain stability and reliability in its predictive outcomes.

## 4.2 BT-ACM+TS+SVM approach is at worst comparable with ACM+TS+SVM

As we analyzed above, the BT-ACM+TS+SVM algorithm shows superior performance compared to the ACM+TS+SVM algorithm in the case of the right hand and left hand versus classification (Figure A1 (c)), while it shows statistically comparable results in the case of the right hand vs left hand task (Figure A1 (c)).

Furthermore, the BT-ACM+TS+SVM algorithm proves its adaptability and robustness by effectively addressing the challenge of rank deficiency often encountered in ACM matrices. This issue becomes particularly pronounced when dealing with datasets that incorporate a high number of electrodes, combined with a high order parameter.

We also noticed improved stability over changes of the SVM parameter, which thus does not need to be optimized as we did in the TS+SVM and ACM+TS+SVM cases (see Table A1).

### 4.2.1 Improvement of computational efficiency

In this section, we look at the computational performance of the proposed approach against the ACM methodology, with a specific focus on two key aspects: the execution time and its environmental impact. For conducting this analysis, we have run all the algorithms on the same hardware, a Dell C6420 dual-Xeon Cascade Lake SP Gold 6240 @ 2.60GHz.

In order to conduct a fair comparison, we consider ACM+TS+SVM and BT-ACM+TS+SVM with the same number of parameters to optimize, i.e. we optimize the order and the lag of the augmentation procedure in the range [1–10] without any optimization of the SVM parameter.

Computational performance was tested on the dataset BNCI2014001 composed of 9 subjects over 2 sessions. Figure 5 (a) shows the differences in computational time. We see that the
difference in computational time is statistically significant as we pass from a mean time of (94.58 ± 2.62)*s* for the ACM methodology versus (79.71 ± 0.80)*s* for the BT-ACM approach. The results have the same level of statistical significance for carbon emissions, evaluated using CodeCarbon package [37].

This behavior is even more marked when using the Lee2019-MI dataset that has a higher number of electrodes. Figure 5 (b) shows the differences in computational time. The difference in computational time is statistically significant as we pass from a mean time of (1042.00 ± 89.80)*s* for the ACM methodology versus (397.78 ± 5.75)*s* for the BT-ACM approach. The results have the same level of statistical significance for carbon emissions.

It is also noteworthy that the ACM methodology exhibits a significant variability in computational times across different sessions and subjects, indicating a fluctuation in performance consistency. In contrast, the BT-ACM approach demonstrates enhanced stability, showcasing a more uniform and predictable computational time.

These results underscore not only the technical superiority of our proposed methodology in terms of speed and efficiency but also its contribution towards reducing the environmental footprint.

### 4.2.2 Subspace Analysis

To provide a more detailed evaluation of the enhancements attributed to the new algorithm, we specifically examine its performance in comparison to traditional Riemannian mapping and the application of ACM. This analysis is conducted using the BNCI2014001 dataset, with a particular focus on subject 6 during session T. This selection is motivated by the observation that, in this scenario, the BT-ACM+TS+SVM algorithm achieves the highest ROC-AUC score among the algorithms under consideration.

To highlight the change in feature space representation due to BT-ACM, we use the *t*-SNE algorithm [38] to provide a two-dimensional representation. For a reliable illustration of the data topology, we experimented with different values of the perplexity parameter [39]. We selected 35 as perplexity but other values give a similar mapping, thus not influencing the conclusion.

The results are shown in Figure 6, where it is evident that the BT-ACM+TS+SVM algorithm significantly alters the feature space, resulting in a pronounced circular structure.

## 4.3 Hyper-parameter stability

The model presented in this study depends on 2 hyper-parameters: the embedding dimension and the lag. To gain a deeper understanding of their impact on model performance and to assess the stability of these parameters, we conducted a comprehensive analysis of their spread over several datasets.

The results are represented by Figure 7, which delineates the distribution of optimal hyperparameters across the entirety of the grid search space. In particular for the BNCI2014001 dataset (Figure 7 (a)), we observe a clustering of optimal hyperparameters towards lower orders, whereas the lag parameter exhibits a broader distribution.

However, this conclusion is not applicable to the BNCI2014004 dataset (Figure 7 (b)). There is a noticeable shift towards higher order and lag values. This divergence in optimal hyper-parameter configurations between the two datasets may be attributable to the differences in their characteristics, notably the number of EEG channels (22 for the former and 3 for the latter). However, by recomputing results for the BNCI2014001 dataset selecting the same electrodes as in BNCI2014004 (Figure 7 (c)), we get results aligned with Figure 7 (a) thus showing that the number of electrodes is not the only parameter affecting the distribution of hyper-parameters.
# 5 Conclusion

Throughout this research, we focused on the uses of the Block-Toeplitz Augmented Covariance matrix for MI classification. This methodology extends the current ACM by comprehensively using the Block-Toeplitz properties of the BT-ACM matrix. This approach transforms the classification challenge, enabling a distinct analysis within the separate domain of SPD and Siegel Disk matrices.

This procedure achieves performance generally superior to the state-of-the-art, or at worst comparable. However, the improvement over the standard ACM is not only in terms of performance but also in terms of significant reductions in computational costs and carbon emissions.

An interesting future exploration emerges from the fact that we are not directly using the property of the blocks contained in the BT-ACM matrix since the metric does not belong to any known manifold. In order to treat the problem, we were forced to use the Verblusky coefficient transform which introduces possible errors and complications. This means that it might be interesting to develop the mathematical framework for the direct treatment of such BTSPD matrices.

A discernible limitation of our proposed algorithm is its performance in comparison to the TS+EL model. While the results achieved by our algorithm and TS+EL are statistically similar, indicating comparable performance and efficacy in the application, there exists a notable disparity in computational efficiency. Specifically, TS+EL demonstrates a superior computational performance, characterized by faster processing times and reduced computational resource consumption. Addressing this limitation will be a focal point of our future work, aiming to refine algorithm architecture to achieve a more balanced trade-off between performance and computational cost. However, we attempted to combine the BT-ACM with EL, resulting in lower performance compared to the proposed approach.

# 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/BL-ACM. In the repository, we have made available the code along with the results table. Moreover, all the code will be submitted to *PyRiemann* upon the publication.

# References

[1] Henri Lorach et al. "Walking naturally after spinal cord injury using a brain-spine interface". In: *Nature* (2023), pp. 1-8.

[2] Alexandre Barachant et al. "Riemannian geometry applied to BCI classification". In: *International conference on latent variable analysis and signal separation*. Springer. 2010, pp. 629-636.

[3] Wolfgang Förstner and Boudewijn Moonen. "A metric for covariance matrices". In: *Geodesy-the Challenge of the 3rd Millennium*. Springer, 2003, pp. 299-309.

[4] Igor Carrara and Théodore Papadopoulos. "Classification of BCI-EEG based on augmented covariance matrix". In: *arXiv preprint arXiv:2302.04508* (2023).
[5] Ben Jeuris and Raf Vandebril. "The Kahler mean of block-Toeplitz matrices with Toeplitz structured blocks". In: *SIAM Journal on matrix analysis and applications* 37.3 (2016), pp. 1151-1175.

[6] Yann Cabanes. "Multidimensional complex stationary centered Gaussian autoregressive time series machine learning in Poincaré and Siegel disks: application for audio and radar clutter classification". PhD thesis. Université de Bordeaux, 2022.

[7] Holger Dette and Jens Wagener. "Matrix measures on the unit circle, moment spaces, orthogonal polynomials and the Geronimus relations". In: *Linear algebra and its applications* 432.7 (2010), pp. 1609-1626.

[8] Sylvain Chevallier et al. "The largest EEG-based BCI reproducibility study for open science: the MOABB benchmark". In: *arXiv preprint arXiv:2404.15319* (2024).

[9] Bruno Aristimunha et al. "Mother of all BCI Benchmarks". In: *Zenodo* (2023).

[10] Floris Takens. "Detecting strange attractors in turbulence". In: *Dynamical systems and turbulence, Warwick 1980*. Springer, 1981, pp. 366-381.

[11] Jaroslav Stark et al. "Takens embedding theorems for forced and stochastic systems". In: *Nonlinear Analysis: Theory, Methods & Applications* 30.8 (1997), pp. 5303-5314.

[12] Bernd Fritzsche and Bernd Kirstein. "An extension problem for non-negative Hermitian block Toeplitz matrices. II". In: *Mathematische Nachrichten* 131.1 (1987), pp. 287-297.

[13] Frederic Barbaresco. "Interactions between symmetric cone and information geometries: Bruhat-tits and siegel spaces models for high resolution autoregressive doppler imagery". In: *LIX Fall Colloquium on Emerging Trends in Visual Computing*. Springer. 2008, pp. 124-163.

[14] Werner Ballmann. *Lectures on Kähler manifolds*. Vol. 2. European mathematical society, 2006.

[15] Le Yang. "Medians of probability measures in Riemannian manifolds and applications to radar target detection". PhD thesis. Université de Poitiers, 2011.

[16] Yilun Chen et al. "Shrinkage algorithms for MMSE covariance estimation". In: *IEEE transactions on signal processing* 58.10 (2010), pp. 5016-5029.

[17] Oncel Tuzel, Fatih Porikli, and Peter Meer. "Pedestrian detection via classification on riemannian manifolds". In: *IEEE transactions on pattern analysis and machine intelligence* 30.10 (2008), pp. 1713-1727.

[18] Michael Tangermann et al. "Review of the BCI competition IV". In: *Frontiers in neuroscience* (2012), p. 55.

[19] David Steyrl et al. "Random forests in non-invasive sensorimotor rhythm brain-computer interfaces: a practical and convenient non-linear classifier". In: *Biomedical Engineering/Biomedizinische Technik* 61.1 (2016), pp. 77-86.

[20] R. Leeb et al. "Brain-computer communication: motivation, aim, and impact of exploring a virtual apartment". In: *Neural Systems and Rehabilitation Engineering, IEEE Transactions on* 15.4 (2007), pp. 473-482.

[21] Josef Faller et al. "Autocalibration and recurrent adaptation: Towards a plug and play online ERD-BCI". In: *IEEE Transactions on Neural Systems and Rehabilitation Engineering* 20.3 (2012), pp. 313-319.

[22] Hohyun Cho et al. "EEG datasets for motor imagery brain-computer interface". In: *Giga-Science* 6.7 (2017), gix034.

[23] Min-Ho Lee et al. "EEG dataset and OpenBMI toolbox for three BCI paradigms: An investigation into BCI illiteracy". In: *GigaScience* 8.5 (2019), giz002.
[24] Gerwin Schalk et al. "BCI2000: a general-purpose brain-computer interface (BCI) system". In: *IEEE Transactions on biomedical engineering* 51.6 (2004), pp. 1034-1043.

[25] Weibo Yi et al. "Evaluation of EEG oscillatory patterns and cognitive process during simple and compound limb motor imagery". In: *PloS one* 9.12 (2014), e114853.

[26] Bangyan Zhou et al. "A fully automated trial selection method for optimization of motor imagery based brain-computer interface". In: *PloS one* 11.9 (2016), e0162657.

[27] Chang S Nam, Anton Nijholt, and Fabien Lotte. *Brain-computer interfaces handbook: technological and theoretical advances*. CRC Press, 2018.

[28] Gavin C Cawley and Nicola LC Talbot. "On over-fitting in model selection and subsequent selection bias in performance evaluation". In: *The Journal of Machine Learning Research* 11 (2010), pp. 2079-2107.

[29] Fabien Lotte and Cuntai Guan. "Learning from other Subjects Helps Reducing Brain-Computer Interface Calibration Time". In: *International Conference on Audio Speech and Signal Processing (ICASSP)*. Dallas, United States, Mar. 2010, pp. 614-617. URL: http://hal.inria.fr/inria-00441670.

[30] Kai Keng Ang et al. "Filter bank common spatial pattern (FBCSP) in brain-computer interface". In: *2008 IEEE international joint conference on neural networks (IEEE world congress on computational intelligence)*. IEEE. 2008, pp. 2390-2397.

[31] Marie-Constance Corsi et al. "Functional connectivity ensemble method to enhance BCI performance (FUCONE)". In: *IEEE Transactions on Biomedical Engineering* (2022).

[32] Robin Tibor Schirrmeister et al. "Deep learning with convolutional neural networks for EEG decoding and visualization". In: *Human brain mapping* 38.11 (2017), pp. 5391-5420.

[33] Vernon J Lawhern et al. "EEGNet: a compact convolutional neural network for EEG-based brain-computer interfaces". In: *Journal of neural engineering* 15.5 (2018), p. 056013.

[34] Thorir Mar Ingolfsson et al. "EEG-TCNet: An accurate temporal convolutional network for embedded motor-imagery brain-machine interfaces". In: *2020 IEEE International Conference on Systems, Man, and Cybernetics (SMC)*. IEEE. 2020, pp. 2958-2965.

[35] Abbas Salami, Javier Andreu-Perez, and Helge Gillmeister. "EEG-ITNet: An Explainable Inception Temporal Convolutional Network for Motor Imagery Classification". In: *IEEE Access* 10 (2022), pp. 36672-36685.

[36] Xia Chen et al. "Toward reliable signals decoding for electroencephalogram: A benchmark study to EEGNeX". In: *arXiv preprint arXiv:2207.12369* (2022).

[37] Courty Benoit et al. *CodeCarbon v2.2.7*. Zenodo. July 2023. URL: https://doi.org/10.5281/zenodo.8181237.

[38] Laurens Van der Maaten and Geoffrey Hinton. "Visualizing data using t-SNE." In: *Journal of machine learning research* 9.11 (2008).

[39] Martin Wattenberg, Fernanda Viégas, and Ian Johnson. "How to use t-SNE effectively". In: *Distill* 1.10 (2016), e2.
<figure><img src="image_6.png" /><figcaption>Figure 4: 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 BT-ACM+TS+SVM against respectively ACM+TS+SVM, TS+EN, TS+SVM, 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.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5: Performance (ROC AUC) versus Computational Time for 1-Fold. Each point represents the result of a subject in a session. The size of the blob indicates CO₂ (g CO₂ equivalent) consumptions. (a) Analysis for the dataset BNCI2014001. (b) Analysis for the dataset Lee2019-MI. The plot shows the statistical significance level where * stands for $p < 0.05$, ** for $p < 0.01$, and *** for $p < 0.001$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Visualization of EEG features using two-dimensional *t*-SNE projection on the binary classification task for right versus left hand. The results represent the *t*-SNE after the tangent space mapping using different algorithms using subject 6 of dataset BNCI2014001 over the session T. The results are obtained for a perplexity value of 35.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7: (a) Heatmap representing the number of times that a certain combination of hyper-parameter of lag and order was selected as best between folds, session and subject in the BNCI2014001 dataset. (b) Heatmap representing the number of times that a certain combination of hyper-parameter of lag and order was selected as best between folds, session and subject in the BNCI2014004 dataset. (c) Heatmap representing the number of times that a certain combination of hyper-parameter of lag and order was selected as best between folds, session and subject in the BNCI2014001 dataset with the selection of the 3 electrodes common with BNCI2014004.</figcaption></figure>
A Supplementary Material

<table><caption>Table A1: Parameters for ML pipelines. * For dataset Cho2017, Lee2019-MI, PhysionetMI and Weibo the parameter were optimized over the range [1, 2, 3, 4, 5] due to computational constraint.</caption><thead><tr><th>Pipeline</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>"FBCSP+SVM"</td><td>C</td><td>numpy.logspace(-2, 2, 10)</td></tr><tr><td>"TS+EN"</td><td>l1 ratio</td><td>[0.2, 0.30, 0.45, 0.65, 0.75]</td></tr><tr><td>"TS+SVM"</td><td>C<br>Kernel</td><td>[0.5, 1, 1.5]<br>"linear", "rbf"</td></tr><tr><td>"ACM+TS+SVM"</td><td>C<br>Kernel<br>Order *<br>Lag *</td><td>[0.5, 1, 1.5]<br>"linear", "rbf"<br>[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]<br>[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]</td></tr></tbody></table>
<table><caption>Table A2: Table listing hyper parameters used in the DL pipeline, these hyper parameters are consistent across all DL pipelines. These hyper parameters represent widely accepted and established settings commonly used in the Deep Learning community.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Epoch</td><td>300</td></tr><tr><td>Batch Size</td><td>64</td></tr><tr><td>Validation Split</td><td>0.2</td></tr><tr><td>Loss</td><td>Sparse Categorical Crossentropy</td></tr><tr><td>Optimizer</td><td>Adam<br>Learning Rate = 0.001</td></tr><tr><td>Callbacks ES</td><td>Early Stopping<br>Patience = 75<br>Monitor = Validation Loss</td></tr><tr><td>Callbacks LR</td><td>ReduceLROnPlateau<br>Patience = 75<br>Monitor = Validation Loss<br>Factor = 0.5</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>Figure A1: Results for Right vs Feet 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 BT-ACM+TS+SVM against respectively ACM+TGSP+SVM, TS+EN, TS+SVM, 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.</figcaption></figure>