MFMSRNet: An Interpretable Multi-frequency and Multi-scale Riemannian Network for Motor Imagery Decoding

1 Introduction

Brain computer interface (BCI) based on electroencephalogram (EEG) has achieved notable progress in motor imagery (MI) decoding, driven largely by advances in deep learning architectures [Ajmeria et al., 2022; Altaheri et al., 2023]. Convolutional and attention-based neural networks enable automatic feature extraction from raw EEG signals and have demonstrated competitive classification performance across benchmark datasets [Zhang and Li, 2025]. However, the black-box nature of these models poses a significant challenge for practical deployment, as they provide limited transparency into the learned features and underlying decision processes [Rudin and Cynthia, 2019]. This lack of interpretability not only hinders neuroscientific insight but also limits reliability in safety-critical applications, where understanding the mapping between brain dynamics and model outputs is essential for clinical trust [Xu and Yang, 2025].

Beyond interpretability, a fundamental limitation arises from the geometric mismatch between EEG representations and the learning space [Yger et al., 2017]. Multi-channel EEG signals exhibit rich inter-channel dependencies, often represented through covariance or functional connectivity (FC) matrices [Varoquaux et al., 2010; Dai et al., 2019]. Standard neural architectures treat multi-channel EEG as vectors or tensors in Euclidean space, implicitly assuming linearity and disregarding the true geometry of connectivity representations [He and Wu, 2020]. Nevertheless, covariance descriptors, correlation patterns, and FC matrices reside on non-Euclidean manifolds, specifically the manifold of symmetric positive definite (SPD) matrices, where conventional Euclidean operations lack geometric fidelity. Prior studies have shown that Riemannian geometry can be leveraged to better preserve the intrinsic structure of EEG features [Tibermacine et al., 2024; Paillard et al., 2025]. For example, Riemannian geometry-based spatial filtering can enhance discriminative power while maintaining robustness to noise compared to Euclidean counterparts, highlighting the importance of respecting geometric structure in EEG classification tasks [Pan et al., 2025].

Despite the advantages of Riemannian geometry in EEG representation learning, existing manifold-based methods for MI decoding have several key limitations [Bronstein et al., 2017; Gerken et al., 2023]. First, most methods focus on covariance or single-band connectivity, ignoring the multi-frequency nature of MI-related neural dynamics [Wen et al., 2024]. Different frequency bands contribute unequally to MI discrimination, especially the mu (8-13 Hz) and beta (13-30 Hz) bands, where event-related desynchronization (ERD) is most prominent and informative. Current multi-frequency fusion strategies often use simple concatenation or Euclidean attention, which can distort the intrinsic geometry of SPD matrices [Suh and Hyung, 2021; Dissanayake et al., 2022]. Second, prior Riemannian networks typically model functional connectivity at a single (whole-brain) spatial scale, failing to capture the hemispheric lateralization and localized sensorimotor interactions critical for MI tasks [Tanamachi et al., 2023]. Although recent geometric deep learning models introduce manifold-aware mappings or attention mechanisms, they rarely integrate multi-scale spatial structure within a unified framework [Shi et al., 2025; Fei et al., 2025]. Finally, while manifold representations are often regarded as more interpretable, most existing works provide limited neurophysiological validation beyond performance improvements, leaving the relationship between learned geometric features and established brain network properties largely unexplored [Dan et al., 2022; Atz et al., 2021].

To address these challenges, this work introduces Multifrequency and Multi-scale Riemannian Network (MFMSRNet), an end-to-end geometry-aware framework that integrates kernelized FC, Riemannian attention-based fusion, and multiscale geometric representations for MI decoding. By constructing FC matrices on the SPD manifold and fusing multifrequency information through attentive geometric operations, the proposed method not only achieves competitive performance across intra-subject and inter-subject benchmarks but also yields representations that are more interpretable in terms of underlying neural interactions. Our contributions can be summarized as follows:

• We propose MFMSRNet, an end-to-end Riemannian deep learning framework for MI EEG decoding that jointly models FC, multi-frequency fusion, and multiscale spatial structure on the manifold.
• We introduce a kernelized phase-synchrony FC representation that captures non-linear inter-channel dependencies while ensuring SPD, enabling principled Riemannian analysis of EEG connectivity.
• We develop an attention-based Riemannian fusion mechanism that adaptively integrates multi-frequency connectivity representations directly on the manifold, preserving the intrinsic geometric structure to support multiscale spatial learning.
• Beyond improved classification performance, we conduct systematic interpretability analyses, including graph-theoretic small-world topology evaluation and frequency- and region-aware attention profiling, revealing neurophysiologically meaningful patterns underlying the learned representations.

2 RelatedWork

MI EEG decoding has been extensively studied using deep learning methods due to their strong capacity for automatic feature extraction from raw or minimally processed signals. Convolutional neural networks (CNNs) have been widely adopted to model temporal and spatial EEG patterns, with compact architectures such as EEGNet demonstrating competitive performance across multiple BCI benchmarks while maintaining low computational complexity [Lawhern et al., 2018]. Subsequent works extend this paradigm by introducing multi-scale convolutions, frequency-aware branches, and attention mechanisms to better capture spectral variability and inter-channel dependencies [Mane et al., 2020;

Song et al., 2022]. Transformer-based architectures further enhance long-range dependency modeling and provide visualization capabilities through attention maps [Zhang and Li, 2025; Fu et al., 2025; Zhang et al., 2023]. Despite these advances, most deep learning approaches represent multi-channel EEG signals in Euclidean space, implicitly assuming linear relationships between channels [Zhang and Etemad, 2023]. Such representations are insufficient for modeling structured interchannel interactions and often obscure the underlying FC patterns critical to MI decoding. These limitations motivate the development of connectivity-aware and geometry-preserving frameworks, as adopted in MFMSRNet, which explicitly models EEG interactions through FC on non-Euclidean manifolds.

To address the geometric mismatch between EEG connectivity representations and Euclidean learning spaces, Riemannian geometry-based methods have been proposed to operate directly on SPD matrices. Early studies demonstrated that covariance-based features modeled on SPD manifolds yield improved robustness and discriminability for MI classification [Yger et al., 2017; Kalaganis et al., 2022]. More recent approaches integrate Riemannian geometry into deep architectures via bilinear mappings, eigenvalue rectification, and Log-Euclidean projections, enabling end-to-end learning of manifold-valued representations [Suh and Hyung, 2021; Dissanayake et al., 2022]. Extensions of this paradigm include tensor-based Riemannian networks and SPD-valued graph models that incorporate time-frequency information to capture localized neural dynamics [Ju and Guan, 2022; Ju and Guan, 2024]. However, most existing Riemannian methods focus on covariance descriptors or single-band representations, and multi-frequency information is often fused through handcrafted strategies or Euclidean operations, potentially compromising the intrinsic manifold structure. In contrast, the proposed MFMSRNet performs attention-guided fusion of multi-frequency FC directly on the Riemannian manifold, preserving geometric consistency while enabling adaptive spectral integration.

FC analysis provides a principled way to characterize interregional neural interactions and has been increasingly combined with attention mechanisms to enhance MI decoding. Attention-based CNNs and feature selection networks emphasize informative frequency bands or spatial regions, leading to improved classification performance [Altaheri et al., 2022; Qin et al., 2024]. More recent works introduce self-attention mechanisms on correlation or covariance manifolds, enabling adaptive weighting of connectivity patterns under Riemannian constraints [Wang et al., 2024; Hu et al., 2025]. Nevertheless, most existing attention-based FC methods focus on either frequency-level or region-level selection in isolation, and few consider connectivity patterns across multiple spatial scales within a unified framework. Moreover, although manifold-based representations are often claimed to be more interpretable, systematic validation against established neurophysiological network properties remains limited. MFMSRNet addresses these limitations by jointly modeling global, hemispheric, and local connectivity structures on the SPD manifold, while providing quantitative interpretability analyses based on graph-theoretic measures and attention-based frequency and region profiling.

3 ProposedMethod

This work proposes a unified end-to-end framework for decoding MI EEG based on multi-frequency FC matrices construction, Riemannian geometry-aware feature learning, attentionguided frequency fusion, and multi-scale connectivity modeling. The entire pipeline (Figure 1) preserves the intrinsic geometry of FC matrices-modeled as SPD matrices-while enabling adaptive and discriminative representation across frequency bands.

3.1 Kernelized Phase-Synchrony FCMatrices

Given multi-channel EEG signals X ∈ RN×C×T , where N denotes the number of trials, C the number of channels, and T the number of time points, band-pass filtering is first applied to extract frequency-specific components.

To characterize functional interactions between EEG channels, we adopt phase-based connectivity modeling. Phase synchronization has been shown to be particularly suitable for MI EEG analysis, as it is less sensitive to amplitude variability and volume conduction effects, while effectively capturing oscillatory coupling underlying MI processes [Liu et al., 2016]. Accordingly, for each frequency band Xn, analytic signals are obtained via the Hilbert transform, and the instantaneous phase ϕc(t) is extracted for each channel.

The phase-locking value (PLV) Pij between channels i and j is computed as:

where i represents the imaginary unit, Pij ∈ [0, 1], with larger values indicating stronger phase synchronization.

Although PLV provides a robust measure of phase coupling, it yields a bounded similarity that does not directly conform to the geometric structure required by Riemannian learning. To address this limitation, we propose a Kernelized PLV (KPLV) representation, which maps PLV-based relationships into the space of SPD matrices.

Specifically, PLV is first converted into a distance measure:

A Gaussian kernel is then applied to obtain a similarity matrix:

where σ controls the decay rate of similarity, which is empirically determined by grid search on the verification set. This kernel operation introduces nonlinearity and ensures that the resulting FC matrix is positive definite, thereby enabling principled Riemannian geometry-based analysis. Compared with raw PLV or linear correlation, KPLV captures higher-order dependencies while explicitly preserving the manifold structure of FC.

To enhance numerical stability, we apply a regularization term:

where λ is a small positive constant. The matrix is then tracenormalized:

and symmetrized to eliminate numerical asymmetry introduced by floating-point computation:

Each Kfinal corresponds to the kernelized FC of a specific frequency band. Let Fn = K(n) final denote the SPD matrix for the n-th band. The set {Fn}Bn=1 then serves as input to the subsequent Riemannian feature learning and attention-based fusion modules, where B is the number of frequency bands.

3.2 Attention-BasedMulti-frequency FCMatrices Fusion

Multi-frequency fusion is essential for MI decoding because different frequency ranges contribute unequally to motor execution and imagery, particularly in the mu and beta bands. However, directly fusing multi-branch SPD representations is non-trivial, as naive concatenation or Euclidean averaging may violate the underlying manifold geometry and fail to capture frequency-specific relevance. To address this, we design an SPD-constrained attention-based fusion mechanism operating in the Log-Euclidean domain, as illustrated in Figure 2.

To fuse SPD matrices, each encoded matrix is projected to the tangent space via the matrix logarithm:

This mapping converts the non-Euclidean SPD manifold into a Euclidean vector space where linear operations are valid. To obtain a compact representation for attention learning, we extract a band-level embedding from each Ln by concatenating its diagonal elements and row-wise means, denoted as zn.

An attention scoring network g(·), implemented as a lightweight MLP (or alternatively a multi-head self-attention module), is then applied to compute the importance score of each frequency band:

The attention weights are obtained via a softmax normalization:

which satisfies αn ≥ 0 and ∑Bn=1 αn = 1.

This yields interpretable band relevance, highlighting which frequency bands contribute most to MI classification. We then fuse multi-frequency FC in the Log-Euclidean domain:

Finally, the fused SPD matrix is reconstructed via the matrix exponential:

This guarantees that the fused output remains a valid SPD matrix, enabling further processing in Riemannian layers.

3.3 Bilinear and RectificationMatrixMapping Module

The FC matrices extracted by the proposed KPLV module are SPD matrices, which lie on a connected Riemannian manifold denoted by S C ++ [Tuzel et al., 2007]. Directly applying conventional Euclidean operations to such matrices may distort their intrinsic geometry. To enable discriminative learning while preserving manifold structure, we introduce a Bilinear Softplus Rectification (BiSRe) module for SPD-valued feature transformation. The fused FC matrix is treated as an SPD-valued neural network input and processed through two stacked BiSRe layers to perform nonlinear geometric transformations directly on the SPD manifold.

As illustrated in Figure 3, each BiSRe module consists of a Bilinear Mapping (BiMap) layer followed by a SoftPlus Eigenvalue Rectification (ReEig) layer. Let F (k-1) ∈ S C ++ denote the input SPD matrix to the k-th BiSRe layer. The BiMap layer applies a congruence transformation defined as:

where Wk ∈ RC×C is a learnable full-rank matrix initialized close to the identity. This operation maps the input SPD matrix to another SPD matrix without vectorization, thereby preserving the intrinsic functional connectivity structure among EEG channels. The positive definiteness of the output is guaranteed by the property of congruence transformations, i.e., WkFW ⊤ k ≻ 0 whenever F ≻ 0.

To ensure that the mapped matrix remains on the SPD manifold and to introduce nonlinearity into the Riemannian network, an eigenvalue rectification layer is applied after the BiMap transformation.

Specifically, the eigen-decomposition of F (k) b is computed as:

where U and Λn denote the eigenvectors and eigenvalues, respectively. In classical Riemannian neural networks, the ReEig layer introduces nonlinearity via hard-thresholding of small eigenvalues with a predefined constant ε . Although this guarantees strict positive definiteness of the output matrix, it results in a non-smooth eigenvalue transformation, potentially causing unstable gradients when eigenvalues are near the threshold.

To address this, we replace the standard hard thresholding with a Softplus-based eigenvalue rectification. The rectification function is defined as:

where ε > 0 is a small constant to prevent numerical degeneration. The output of the ReEig layer is reconstructed as:

This operation guarantees that all eigenvalues remain strictly positive, thereby preserving the SPD property of the matrix while introducing nonlinearity analogous to the ReLU function in conventional convolutional neural networks. Two BiSRe layers per frequency band enable the encoder to learn non-linear SPD-to-SPD transformations, capturing MI task-relevant connectivity structure within each band.

3.4 Multi-Scale Riemannian Network

After refining Ffused with stacked BiSRe modules to obtain the discriminant SPD representation on the Riemannian manifold, it is further processed by a designed multi-scale Riemannian network to capture the functional interactions related to MI at different spatial resolutions. This design is motivated by neurophysiological evidence that MI modulates brain connectivity at global, hemispheric, and local sensorimotor scales. To facilitate scale-aware modeling, EEG channels are reordered as [Cleft, Cright], with channels within each hemisphere grouped by anatomical regions, enabling structured extraction of multi-scale connectivity patterns on the SPD manifold.

The network consists of three complementary branches.

Whole-brain branch. The full-channel FCmatrix is directly processed by two cascaded BiSRe layers to encode global brain network organization.

Hemispheric branch. To capture lateralized MI patterns, the FC matrix is partitioned into intra- and inter-hemispheric principal submatrices corresponding to left-left, left-right, right-left, and right-right connectivity. According to Theorem 1, each submatrix, which remains SPD as a principal submatrix of Ffused, is embedded using shared BiSRe layers, and the resulting features are concatenated to model hemispheric interactions relevant to left/right MI.

Theorem1. LetM ∈ S C ++. For any index set I ⊂ {1, . . . , C} with |I| = w, the principal submatrix M[I] is in S w ++. (The proof is provided in Appendix A.)

Local branch. Fine-grained regional connectivity is extracted using a sliding-window mechanism on the FC matrix, producing a set of local SPD submatrices that capture micro-scale sensorimotor interactions. These local representations are mapped via shared BiSRe layers and aggregated to characterize localized connectivity patterns.

After scale-specific Riemannian encoding, the resulting SPD matrices are projected onto the tangent space using the Log-Euclidean mapping and vectorized by extracting upper triangular elements. Feature vectors from the whole-brain, hemispheric, and local branches are concatenated and fed into a fully connected classifier for MI decoding.

By jointly aggregating connectivity representations across multiple spatial scales while preserving manifold structure, the proposed network captures complementary global, lateralized, and local brain network characteristics, leading to more discriminative and interpretable MI representations.

4 Experiments

4.1 Datasets and Evaluation Protocols

Experiments are conducted on two complementary MI EEG datasets: the standardized BCI Competition IV 2a (22×22 FC matrix) and a clinical lower-limb MI (LL-MI) dataset from stroke patients (40×40 FC matrix). The two complement each other and are used to evaluate methods in standard benchmarks and real clinical scenarios. Performance is evaluated using accuracy, F1 score, and Cohen's kappa. All tests run on an RTX 5090 D GPU (32GB) and an i7 CPU @ 3.4GHz. The average end-to-end processing time per trial is approximately 0.22s, with core MI inference around 0.1s.

BCI-IV 2a. The dataset [Tangermann et al., 2012] comprises EEG recordings from 9 healthy subjects performing four-class MI tasks: left hand (LH), right hand (RH), both feet (BF), and tongue (Tg). Owing to its standardized experimental design and well-established evaluation protocols, this dataset serves as a canonical benchmark for MI decoding. Both intra-subject and inter-subject evaluation protocols are adopted to assess the robustness and generalization ability of the proposed method in a multi-class setting.

LL-MI. The dataset [Liu et al., 2025] contains EEG recordings from 27 post-stroke patients performing lower-limb MI tasks, representing a challenging and clinically relevant scenario. In this work, we focus on the gait-phase-encoded sequential electrical stimulation paradigm, which includes gait MI and idle-state tasks. Evaluation is conducted across subjects within this paradigm to examine the robustness and interpretability of the proposed framework in rehabilitationoriented BCI applications.

4.2 Comparison with Prior Art

The results in Table 1 show that our MFMSRNet consistently surpasses representative baselines in both intra- and inter-subject evaluations on BCI-IV 2a and LL-MI datasets. In both settings, Euclidean-based models (e.g., EEGNet, EEG-Inception) perform less effectively. This limitation suggests that shallow spatial-temporal convolutions alone are inadequate for capturing discriminative patterns in multiclass MI. Transformer-hybrid architectures such as EEGConformer achieve moderate improvements but lag behind geometry-aware and attention-enhanced methods, particularly in the more challenging inter-subject evaluation, where crossindividual variability dominates.

Tensor-CSPNet and Graph-CSPNet, which incorporate geometric insights into SPD manifold structures, deliver stronger performance than purely Euclidean networks, yet they remain constrained by their focus on second-order statistics without adaptive frequency fusion. Attention-driven models such as MFANet and CorAtt-MIX further narrow this gap by leveraging multi-feature attention or correlation manifold self-attention, demonstrating the value of adaptive weighting mechanisms. However, these approaches still do not fully exploit both FC geometry and adaptive multi-frequency integration.

Notably, our MFMSRNet achieves the highest accuracy, F1 scores, and Cohen's kappa values across all benchmarks, substantially surpassing CorAtt-MIX and other manifold baselines. On the LL-MI dataset, MFMSRNet attains 94.69%

accuracy, 94.81% F1 score, and 0.8938 kappa, indicating not only strong discriminative capability but also robust generalization in a binary classification paradigm. These consistent improvements underscore that the proposed MFMSRNet yields a more expressive and resilient representation of EEG dynamics compared to existing state-of-the-art methods.

4.3 Explainability Analysis

The proposed framework produces fused FC representations with inherent interpretability by explicitly modeling EEG interactions on the Riemannian manifold of SPD matrices, enabling direct neurophysiological analysis [Wang et al., 2025].

Graph-theoretic Interpretability Analysis

To assess whether the fused FC representations preserve meaningful brain organization, we conduct a graph-theoretic analysis based on network topology metrics (see Appendix B for metric calculations), with a focus on small-world properties [Gallego-Molina et al., 2022]. As shown in Figure 4, all fused FC networks exhibit significant small-world properties (> 1), indicating a balance between local segregation and global integration [Bassett and Bullmore, 2017].

For the BCI-IV 2a dataset (Figure 4 (a)), hand MI shows significantly higher small-world properties than Tg and BF imagery (p < 0.05), linked to greater local clustering without compromising path length - consistent with the denser cortical representation of hand movements. No significant differences were observed between Tg and BF imagery or between LH and RH.

For the LL-MI dataset (Figure 4 (b)), the Idle state consistently shows a significantly higher small-world property than Gait MI (p < 0.01). This reflects higher clustering and shorter path lengths at rest, while Gait MI reduces clustering and slightly increases path length, indicating task-related network reconfiguration [Neyland et al., 2021; Su et al., 2025].

Frequency- and Region-aware Interpretability Analysis

We further investigate the interpretability of the proposed attention-based Riemannian fusion from spectral and regional perspectives. In the LL-MI dataset (Figure 5 (a)), Idle is dominated by mu and low-beta bands, while Gait MI assigns higher attention to high-beta and gamma bands, reflecting task-relevant spectral adaptation.

Region-wise FC analysis (Figure 5 (b)), after trace normalization and significance filtering, reveals distinct spatial patterns. Gait MI shows stronger connectivity within the central sensorimotor region, reflecting enhanced sensorimotor integration during gait imagery. In contrast, Idle shows more distributed connectivity across bilateral motor and frontal-parietal regions, associated with internally oriented processing without explicit motor demands.

Similar trends appear in the BCI-IV 2a dataset (Figure 5 (c, d)). All MI classes are dominated by mu and beta rhythms, with BF and Tg imagery showing higher contributions from higher-frequency bands. Attention distributions are symmetric between LH and RH MI. Region-wise FC reveals contralateral motor cortex dominance during hand imagery, while foot and tongue imagery involve stronger parietal and parieto-occipital engagement.

Unlike prior methods that provide mainly post-hoc or feature-level interpretability, MFMSRNet offers intrinsic connectivity-level explanations grounded in Riemannian geometry, enabling frequency- and region-specific neurophysiological interpretation.

5 Ablation Studies

To verify the effectiveness of multi-frequency FC fusion, ablation studies are conducted comparing the performance of individual frequency bands against the fused FC. As shown in Figure 6 (sparsity level fixed at 0.05 for both datasets), while individual bands exhibit certain advantages in capturing specific inter-group differences, the fused FC comprehensively integrates unique discriminative features from all bands, leading to significantly improved overall classification performance.

The decoding network further improves recognition performance by fusing FC across whole-brain (Sw), hemispheric

(Sh), and local region (Sl) scales. Ablation experiments confirm the effectiveness of this multi-region selection mechanism, showing that models using each scale alone or in partial combinations are outperformed by the full three-scale fusion (Table 2). Specifically, on the BCI-IV 2a dataset, the combined Sw, Sh, and Sl fusion yields the best results across all metrics, while achieving 94.69% accuracy on the LL-MI dataset. These findings demonstrate that the proposed coarse-to-fine feature fusion strategy substantially enhances model performance.

6 Conclusion

This work presents MFMSRNet, an interpretable end-to-end MI decoding framework that integrates kernelized FC construction, attention-based Riemannian fusion, and multi-scale geometric feature learning. By modeling EEG FC as SPD matrices and performing adaptive fusion directly on the Riemannian manifold, the proposed method effectively captures non-linear dependencies and preserves the intrinsic geometric structure of brain networks. Extensive experiments on two public MI datasets demonstrate that the proposed framework consistently outperforms state-of-the-art baselines in both intra- and inter-subject settings. Moreover, comprehensive explainability analyses, including small-world topology, frequency-specific attention patterns, and region-level connectivity strength, reveal neurophysiologically meaningful representations aligned with established motor control mechanisms.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (No.62372338) and the Key Research and Development Projects in Hubei Province (No.2023BCB133).