Routing on the Stiefel Manifold:
When Does Adaptive Subspace Selection Help for Cross-Domain EEG Decoding?
May 29, 2026
Abstract
Cross-domain EEG decoding remains challenging despite advances in Riemannian deep learning: covariance matrices from different subjects occupy systematically distinct regions of the SPD manifold, yet existing domain adaptation methods either require target-domain calibration data or learn subject-specific components that cannot generalise across domains. We propose to address this through dynamic Stiefel routing: a pool of \(K\) expert projection filters on the Stiefel manifold, each specialised for a different region of the SPD manifold, with each input covariance routed to the most appropriate filter via cross-attention — adapting the subspace projection per sample. A central finding of this work is that this approach, implemented naively, provably collapses to ensemble averaging: when routing weights are uniform, the adaptive filter reduces exactly to an equal-contribution combination of experts, indistinguishable from a single fixed filter. We identify the conditions under which this degeneracy can be broken. Three structural properties are necessary and sufficient: a symmetric anchor \(W_\text{base} \in \mathrm{St}(n,k)\) that removes proximity bias among experts; a frozen domain-discriminative query encoder — realised as a domain-separation projection of the log-covariance tangent space — that decouples routing from task optimisation; and a decoupled key alignment loss that trains expert keys toward stable domain attractors. Together they produce, to the best of our knowledge, the first genuinely committed and domain-structured routing on SPD manifolds, with consistent gains across three datasets spanning contrasting ratios of dimensionality versus number of domains: balanced accuracy improves from \(0.773{\to}0.823\), \(0.757{\to}0.809\), and \(0.801{\to}0.839\), with the domain-separation projection and alignment strategy determined automatically by a single data-driven rule and no dataset-specific hyperparameter search. A full-rank rectangular extension is indicated, and a concrete implementation path is described.
1 Introduction↩︎
1.0.0.1 The problem.
The introduction of Riemannian geometry to brain-computer interfaces marked a turning point in EEG-based motor imagery decoding. By representing neural signals as covariance matrices on the symmetric positive definite manifold \(\mathrm{SPD}(n)\) and exploiting its rich metric structure, Riemannian methods consistently outperformed conventional feature engineering and established a new state of the art [1], [2]. Deep networks soon followed: SPDNet [3] and its successors extended the Riemannian paradigm to learned representations, stacking BiMap, ReEig, and LogEig layers that respect the manifold geometry throughout the forward pass and brought further performance gains.
Despite this progress, cross-domain decoding remains stubbornly plateaued. The main culprit identified with high inter-subject variability: electrode placement, scalp conductivity, neural synchrony patterns, and oscillatory frequency profiles all differ between subjects, causing their covariance matrices to occupy systematically distinct regions of \(\mathrm{SPD}(n)\). A model trained on many subjects must simultaneously serve these wildly different input distributions at test time — but the standard BiMap filter \(Y = W^\top X W\) with \(W \in \mathrm{St}(n,k)\) projects each covariance matrix \(X \in \mathrm{SPD}(n)\) through the same fixed subspace, no matter its domain of origin. This is a fundamental design mismatch: the network is geometrically aware of the manifold structure within a subject, yet completely blind to the between-subject structure across the manifold. No matter how deep the architecture becomes, a single shared projection cannot account for the subject-specific region of \(\mathrm{SPD}(n)\) each input occupies — and performance saturates as a result. We address this by learning to dynamically route each input covariance to a sample-specific projection at inference time, without retraining or calibration.
1.0.0.2 An idea and its failure.
An idea to address this is to learn \(K\) expert filters \(\{W_j\}_{j=1}^K \subset \mathrm{St}(n,k)\) and route each input sample to the most appropriate expert via a weighted Riemannian barycenter. We implemented this routing on the Stiefel manifold and immediately encountered a striking failure mode: the routing weights collapse to near-uniform (\(\alpha_i \approx \mathbf{1}/K\)), and the method reduces to ensemble averaging of a diverse expert pool. The performance gain over the fixed-filter baseline is real — but it comes entirely from having \(K\) diverse filters whose mean is better than any single one, not from any per-sample adaptation. This failure is not an accident of initialisation or optimisation. We argue it is a consequence of a fundamental gradient degeneracy: the cross-entropy signal reaching the routing weights \(\alpha\) is proportional to the deviation of each expert’s tangent vector from the current weighted mean, which is vanishingly small when experts are similar — and experts remain similar precisely because they receive no differentiated gradient signal to specialise. Routing collapse and expert homogeneity are self-reinforcing.
1.0.0.3 Diagnosing and breaking the degeneracy.
We introduce a key diagnostic that separates these two regimes: the K=1 proxy gap \(\Delta_\text{K=1}\), defined as the difference in balanced accuracy between the adaptive model and the single-expert ensemble baseline evaluated on the same test set. When routing is uniform, \(\Delta_\text{K=1} = 0\) exactly. Positive \(\Delta_\text{K=1}\) is therefore a necessary condition for genuine routing — it cannot be mimicked by ensemble quality alone.
Armed with this diagnostic, we identify three structural properties that taken together break the degeneracy, each motivated by an analogy with Learning to Prompt (L2P) [4] for continual learning:
A symmetric anchor \(W_\text{base}\) that decouples the reference point from the expert pool, removing the \(W_1\)-proximity bias that artificially dampens routing gradients toward one expert.
A frozen domain-discriminative query encoder — a between-domain separation projection (DSP) of the log-covariance tangent space — that provides routing with domain-identity information independently of the task loss, analogous to L2P’s frozen Visual Transformer (ViT) [5] backbone.
A decoupled key alignment loss that trains expert keys exclusively toward stable query attractors, preventing the moving-target problem that arises when keys compete with classification gradients.
Each component is independently motivated, ablated in isolation, and together they produce, to the best of our knowledge, the first domain-structured routing on SPD manifolds.
1.0.0.4 Contributions.
We prove that naive adaptive routing on \(\mathrm{St}(n,k)\) collapses to a degenerate solution (uniform weights, single-expert ensemble filter) and characterise the responsible gradient mechanism.
We propose Domain-Adaptive Stiefel Pool (DASP) layer, a drop-in replacement for the standard BiMap layer that combines the three structural fixes above, and introduce the K=1 proxy gap as a principled diagnostic for genuine routing value.
We derive a unified scaling methodology for the DSP and alignment loss, governed by the tangent-to-domain ratio \(\rho\), that requires no dataset-specific tuning and generalises across three datasets with \(\rho \in \{3.3, 28, 203\}\). Expert count \(K\) follows a sub-linear scaling with \(D\) that is characterised empirically across our datasets.
We extend the expert parametrisation to full-rank rectangular matrices and describe two concrete implementation paths.
2 Routing Collapse: The Degenerate Solution↩︎
2.0.0.1 Setup.
Let \(\mathcal{D} = \{(X_i, y_i, d_i)\}_{i=1}^N\) where \(X_i \in \mathrm{SPD}(n)\) is a pre-conditioned EEG covariance matrix, \(y_i \in \{1,\ldots,C\}\) is the class label, and \(d_i \in \{1,\ldots,D\}\) is the domain index. Our Domain-Adaptive Stiefel Pool layer maintains \(K\) expert filters \(\{W_j\}_{j=1}^K \subset \mathrm{St}(n,k)\) and computes attention weights \[\label{eq:alpha} \alpha_i = \mathrm{softmax}(q_i E^\top / \sqrt{m}) \in \Delta^{K-1}\tag{1}\] where \(E \in \mathbb{R}^{K \times m}\) is a learnable matrix of \(K\) expert keys, one per expert filter, trained to align with the query vectors of their respective domains, and \(q_i \in \mathbb{R}^m\) is a query vector dependent on the input \(X_i\) and domain \(d_i\). The \(j\)-th component \(\alpha_{i,j}\) of \(\alpha_i\) gives the weight assigned to expert \(j\) for sample \(i\). The layer then applies the per-sample filter via a Riemannian interpolation on the tangent space \(\mathcal{T}_{W_\text{base}}\mathrm{St}(n,k)\) \[\label{eq:ensemble95expert} W_i = \mathcal{R}_{W_\text{base}}\!\left(\sum_j \alpha_{i,j}\, \mathcal{P}_{W_\text{base}}(W_j)\right)\tag{2}\] where \(\mathcal{R}_{W_\text{base}}(\Delta)\) is the QR retraction of \(\Delta \in \mathcal{T}_{W_\text{base}}\mathrm{St}(n,k)\) back onto the Stiefel manifold, and \(\mathcal{P}_{W_\text{base}}(W)\) is the orthogonal projection of \(W \in \mathrm{St}(n,k)\) onto the tangent space, both at the basepoint \(W_\text{base} \in \mathrm{St}(n,k)\).
2.0.0.2 The degenerate solution.
When \(\alpha_i = \mathbf{1}/K\) for all samples \(i\), the per-sample filter reduces to: \[W_i = \mathcal{R}_{W_\text{base}}\!\Bigl(\tfrac{1}{K}\textstyle\sum_j \mathcal{P}_{W_\text{base}}(W_j)\Bigr) \label{eq:degenerate}\tag{3}\] i.e.every sample is projected through the same fixed filter — the one-step approximation of the Fréchet mean of the expert pool. This is exactly the K=1 proxy: a single fixed filter that carries no sample-specific information.
If \(\alpha_i = \mathbf{1}/K\) for all \(i\), then \(W_i = \mathcal{R}_{W_\text{base}}(\bar{\Delta})\) is constant across samples, with \(\bar{\Delta} = \frac{1}{K}\sum_j \mathcal{P}_{W_\text{base}}(W_j)\), and consequently \(\mathrm{bAcc}_\text{adaptive} = \mathrm{bAcc}_\text{K=1}\), giving \(\Delta_\text{K=1} = 0\).
Since the K=1 proxy uses the same single-expert ensemble filter, the two models produce identical predictions on every input. The K=1 proxy gap is therefore a necessary condition for genuine routing: \(\Delta_\text{K=1} = 0\) implies \(\alpha_i \approx \mathbf{1}/K\) regardless of absolute accuracy.
2.0.0.3 Collapse is self-reinforcing.
The routing weights \(\alpha_i\) are updated by gradients that flow back from the loss through the BiMap output \(Y_i = W_i^\top X_i W_i\). Differentiating \(Y_i\) with respect to the barycenter weights reveals that the sensitivity of the loss to changing \(\alpha_{i,j}\) is proportional to how much expert \(j\) deviates from the current mixture in tangent space — that is, how different \(\mathcal{P}_{W_\text{base}}(W_j)\) is from the interpolated expert \(\bar{\Delta}_i = \sum_j \alpha_{i,j} \mathcal{P}_{W_\text{base}}(W_j)\). When all experts are similar, every expert’s tangent vector is close to the mean, these deviations are small, and the routing weights receive negligible gradient regardless of how informative the upstream signal is. Crucially, experts remain similar precisely because near-uniform routing gives them identical gradient signals and no pressure to specialise — which in turn keeps their tangent vectors close together, keeping routing gradients weak. Routing collapse and expert homogeneity are mutually reinforcing: neither can resolve without the other moving first.
2.0.0.4 Empirical evidence.
Table 2 (row “Adaptive (naive)”) shows this collapse in practice: on Weibo2014 database [6], \(\Delta_\text{K=1} = -0.00\) despite a \(+4.5\%\) gain over the fixed-filter baseline. The gain is real but comes entirely from ensemble diversity — the equal-weight barycenter of \(K\) diverse experts is a better filter than any single one. The routing mechanism contributes nothing.
3 Domain-Adaptive Stiefel Pool (DASP)↩︎
We now describe the three structural fixes that break the degeneracy identified in §2. Each fix has an independent motivation and is validated in isolation in §4.3.
3.1 Fix 1: Symmetric Anchor \(W_\text{base}\)↩︎
We introduce a dedicated learnable parameter \(W_\text{base} \in \mathrm{St}(n,k)\), initialised independently of the expert pool and updated freely by the optimiser while remaining on \(\mathrm{St}(n,k)\) via an orthogonal parametrisation. This symmetrises expert competition: no expert has a geometric proximity advantage from the choice of reference point.
3.2 Fix 2: Frozen DSP Query Encoder↩︎
A query multi-layer perceptron (MLP) trained end-to-end with cross-entropy learns class-discriminative features — precisely the wrong inductive bias for routing, which requires domain-discriminative features. This is the Riemannian analogue of the core insight behind L2P [4]: its frozen ViT backbone provides task-agnostic, domain-informative queries that decouple key learning from task optimisation.
The between-domain signal in tangent space lives in a subspace of rank at most \(D-1\), where \(D\) is the number of training domains. When the tangent space is large relative to \(D\) (high \(\rho = n(n+1)/(2D)\)), this signal is buried in noise and the MLP cannot recover it.
We compute the tangent representation at the identity \(I_n \in \mathrm{SPD}(n)\) as \(x_i = \mathrm{upper}(\log(X_i)) \in \mathbb{R}^{n(n+1)/2}\), where \(\log(\cdot)\) denotes the matrix logarithm and \(\mathrm{upper}(\cdot)\) extracts the upper-triangular part with \(\sqrt{2}\)-scaled off-diagonal entries [7]. Then, we project \(x_i\) through a frozen matrix: \[P_\text{DSP} \in \mathbb{R}^{n(n+1)/2 \times r},\quad \text{columns} = \text{top-}r\text{ DSP filter eigenvectors of }S_B \label{eq:DSP}\tag{4}\] where \(S_B = \sum_{d=1}^D n_d(\mu_d - \mu)(\mu_d-\mu)^\top\) is the between-domain scatter matrix, computed once from training tangent-space representations and frozen. \(P_\text{DSP}\) carries no gradient; it acts as a fixed domain-discriminative feature extractor analogous to L2P’s frozen ViT. We use \(r=40\) to span both the \(D{-}1\) between-domain directions and additional within-domain structure (§4.3). We denote the (optionally projected) tangent vector as \[\tilde{x}_i = \begin{cases} x_i P_\text{DSP} & \text{if } \rho > \rho^* \\ x_i & \text{if } \rho \leq \rho^* \end{cases} \label{eq:xtilde}\tag{5}\] so that \(\tilde{x}_i \in \mathbb{R}^r\) in the high-\(\rho\) regime and \(\tilde{x}_i \in \mathbb{R}^{n(n+1)/2}\) in the low-\(\rho\) regime.
3.3 Fix 3: Decoupled Key Alignment Loss↩︎
Even with a stable, domain-informative query, expert keys trained via cross-entropy receive a moving target: the task loss pushes keys toward whatever representation helps classification, which may not align with domain identity. We instead train keys
exclusively via a surrogate alignment loss: \[\mathcal{L}_\text{align} = \lambda_\text{align}\, \mathbb{E}_i\!\Bigl[1 - \cos\bigl( \mathrm{sg}(q_i),\; E_{j^*_i}\bigr)\Bigr],\quad j^*_i =
\mathop{\mathrm{arg\,max}}_j \alpha_{i,j} \label{eq:align}\tag{6}\] where \(\mathrm{sg}(\cdot)\) denotes stop-gradient on the query (analogous to L2P’s frozen query function), and only the highest-weighted key
\(E_{j^*_i}\) is updated per sample. Pulling all keys toward every query simultaneously would cause key collapse toward the global mean, defeating the purpose of the pool; the winner-take-all update instead creates
competitive pressure for keys to specialise toward their respective domains. When decouple_keys=True, \(E\) is additionally detached from the cross-entropy gradient path entirely, so keys are optimised
exclusively by \(\mathcal{L}_\text{align}\) toward stable per-domain attractors. The total loss is \(\mathcal{L} = \mathcal{L}_\text{CE} + \mathcal{L}_\text{align}\).
3.4 The Full Layer↩︎
The complete DASP layer is shown in Figure 1. It takes as input a batch of SPD matrices \(\{X_i\} \subset \mathrm{SPD}(n)\) and domain indices \(\{d_i\}\), and returns \(Y_i = W_i^\top X_i W_i \in \mathrm{SPD}(k)\), where the per-sample filter \(W_i\) is the attention-weighted Riemannian interpolation of Eq. 2 , with weights \(\alpha_i\) from Eq. 1 . The three components below specify how \(q_i\), \(W_\text{base}\), and the manifold operations are instantiated.
Query construction. The query \(q_i\) (Eq. 1 ) is formed by combining the projected tangent vector \(\tilde{x}_i\) (Eq. 5 ) with a learnable domain embedding \[q_i = \mathrm{MLP}\bigl([\tilde{x}_i;\; \mathrm{Emb}(d_i)]\bigr) \in \mathbb{R}^m \label{eq:query}\tag{7}\] where the MLP has one hidden layer of dimension \(2m\) with GELU activation.
Dedicated anchor. The reference point \(W_\text{base}\) in Eqs. 2 –3 is instantiated as the dedicated learnable parameter \(W_\text{base} \in \mathrm{St}(n,k)\) described in §3.1, replacing the dynamic one-step barycenter of naive routing.
Manifold operations. Since the true Riemannian log-map on \(\mathrm{St}(n,k)\) has no closed form, we use the tangent projection as a first-order approximation \[\mathcal{P}_{W_\text{base}}(W) = W - W_\text{base}\,\mathrm{sym}(W_\text{base}^\top W) \;\in\; \mathcal{T}_{W_\text{base}}\mathrm{St}(n,k) \label{eq:stproj}\tag{8}\] where \(\mathrm{sym}(A) = \tfrac{1}{2}(A + A^\top)\) denotes the symmetric part of \(A\).
The QR retraction maps tangent vectors back onto the manifold: letting \(W_\text{base} + \Delta = QR\) be the thin QR decomposition \[\mathcal{R}_{W_\text{base}}(\Delta) = Q \;\in\; \mathrm{St}(n,k) \label{eq:stexp}\tag{9}\] Both are valid approximations of the true log and exp maps when \(W_\text{base}\) and the target are not too far apart [8], [9]. The BiMap output \(Y_i = W_i^\top X_i W_i \in \mathrm{SPD}(k)\) is then fed to the standard ReEig \(\to\) LogEig \(\to\) Linear tail [3].
3.5 Routing Diagnostics↩︎
Three diagnostics track whether routing is genuine or degenerate:
K=1 proxy gap \(\Delta_\text{K=1} = \mathrm{bAcc}_\text{adaptive} - \mathrm{bAcc}_\text{K=1}\), where \(\mathrm{bAcc}_\text{K=1}\) is the accuracy of \(W_\text{base}^\dagger\) (the single equal-contribution filter) evaluated on the test set. \(\Delta_\text{K=1} = 0\) iff routing is uniform; positive values in a majority of cross-validation folds constitute our success criterion.
Routing entropy \(\bar{H} = \mathbb{E}_i\!\left[\frac{H(\alpha_i)}{\log K}\right] \in [0,1]\): near 1 is uniform (degenerate); near 0 is one-hot (committed).
Domain alignment ratio \(r_d = \mathrm{Var}_d(\bar{\alpha}_d)\,/\,\mathrm{Var}_i(\alpha_i)\), where \(\bar{\alpha}_d = \frac{1}{|I_d|}\sum_{i:\,d_i=d}\alpha_i \in \Delta^{K-1}\) is the mean routing vector for domain \(d\): values above 0.10 indicate domain-structured routing.
3.6 Unified Scaling Methodology↩︎
The degree to which domain identity is recoverable from the tangent space depends on the tangent-to-domain ratio \[\rho = \frac{n(n+1)/2}{D}
\label{eq:ratio}\tag{10}\] When \(\rho \gg 1\), domain identity is buried in high-dimensional noise and the DSP is essential. When \(\rho \approx 1\), the tangent space is
compact relative to the number of domains, and the query MLP finds domain structure unaided. A threshold \(\rho^* = 50\) cleanly separates these regimes across all three datasets, determining whether the DSP projection and
alignment loss are activated. In the high-\(\rho\) regime, expert keys must be decoupled from the cross-entropy gradient (decouple_keys=True) so that the alignment loss can train them exclusively toward stable
domain attractors, as described in §3.3. In the low-\(\rho\) regime, the alignment loss is disabled and decouple_keys=False ensures that keys still receive a training signal
from the classification objective. Expert count \(K\) follows a sub-linear scaling with \(D\): empirically, \(K \approx D-1\) works well for small \(D\) (\(D \leq 10\)) and \(K \approx D/2\) for larger \(D\) (\(D \geq 20\)), reflecting
diminishing returns in routing capacity as the number of domains grows. Exact values per dataset are given in 8. These rules are summarised in Table 1.
| Parameter | Rule | Condition |
|---|---|---|
| \(K\) | sub-linear in \(D\); see 8 | always |
| \(d_\text{emb}\) | \(20\) | always |
| DSP | enabled, \(r=40\) | \(\rho > 50\) |
| DSP | disabled | \(\rho \leq 50\) |
| \(\lambda_\text{align}\) | \(0.05\), decouple_keys=True | \(\rho > 50\) |
4 Experiments↩︎
4.1 Setup↩︎
We evaluate on three motor imagery datasets with contrasting geometry, spanning the low- and high-\(\rho\) regimes. All experiments use the Adam optimiser [10] at learning rate \(0.01\), batch size \(32\), and 5-fold cross-validation with 70/15/15% trial-level train/val/test splits. The primary metric is balanced accuracy (\(\mathrm{bAcc}\)). We consider binary classification between right-hand and feet motor imagery across all datasets. Raw EEG signals are bandpass filtered between 8 and 32 Hz and subjected to an adaptive amplitude thresholding artifact rejection. Covariance matrices are then pre-conditioned per domain by a sequence of three operations: Tikhonov regularization (diagonal loading) with parameter \(10^{-8}\), whitening using the Euclidean mean of the training set and multiplied by the ratio trace/squared Frobenius norm.
Stiefel-constrained parameters are maintained via the orthogonal parametrisation of [11]. SPDNet layers (BiMap, ReEig, LogEig), and
operations like the matrix logarithm are provided by the spd_learn library [12]. The baseline is a standard
single-stage SPDNet [3] with a fixed BiMap layer applying the same subspace dimension \(k\) as our
method.
All experiments were run on a single CPU workstation (Intel Core i7, 32 GB RAM). Each 5-fold cross-validation run takes approximately 20–30 minutes for Weibo2014 (\(n{=}60\)) and 5–10 minutes for the BNCI datasets (\(n \leq 22\)). The full ablation study (Table 3, 6 configurations \(\times\) 5 folds) required approximately 3 hours of compute. No GPU was used.
Weibo2014 [6]: \(n{=}60\), \(D{=}9\), \(K=8\), \(\rho{=}203\). High-\(\rho\) regime; DSP essential.
BNCI2015001 [13]: \(n{=}13\), \(D{=}28\), \(K=16\), \(\rho{=}3.3\). Low-\(\rho\) regime; domain structure naturally accessible.
BNCI2014001 [14]: \(n{=}22\), \(D{=}9\), \(K=8\), \(\rho{=}28\). Low-\(\rho\) regime; sits near the \(\rho^*{=}50\) boundary, testing the scaling rule at an intermediate channel count.
4.1.0.1 Subspace dimension \(k\).
The output dimension \(k\) of the BiMap layer controls the degree of dimensionality reduction applied to the covariance matrices. For Weibo2014 (\(n{=}60\)) we use \(k{=}30\), halving the input dimension following standard practice in SPDNet-based decoding [3]. For the BNCI datasets (\(n{=}13\), \(n{=}22\)) the input dimension is already compact, and aggressive reduction would discard between-domain discriminative structure needed by the routing mechanism. We therefore use more conservative reductions: \(k{=}12\) for BNCI2015001 and \(k{=}20\) for BNCI2014001, preserving a larger fraction of the tangent space geometry. Exact values are reported in 8.
4.2 Main Results↩︎
| Dataset | Model | bAcc | Std | \(\Delta_\text{base}\) | \(\Delta_\text{K=1}\) | Folds | \(\bar{H}\) | \(r_d\) |
|---|---|---|---|---|---|---|---|---|
| Weibo2014 | Baseline BiMap | 0.773 | 0.020 | |||||
| Adaptive (naive) | 0.818 | 0.031 | +0.045 | \(-\)0.00 | 0/5 | 0.836 | 0.281 | |
| DASP (Ours)\(^\dagger\) | 0.823 | 0.011 | +0.050 | +0.029 | 4/5 | 0.674 | 0.710 | |
| BNCI2015001 | Baseline BiMap | 0.757 | 0.017 | |||||
| Adaptive (naive) | 0.796 | 0.012 | +0.040 | +0.050 | 5/5 | 0.494 | 0.888 | |
| DASP (Ours)\(^\dagger\) | 0.809 | 0.009 | +0.053 | +0.247 | 5/5 | 0.482 | 0.539 | |
| BNCI2014001 | Baseline BiMap | 0.801 | 0.023 | |||||
| DASP (Ours)\(^\dagger\) | 0.839 | 0.009 | +0.038 | +0.170 | 5/5 | 0.597 | 0.828 |
Table 2 tells a clear story. The naive adaptive model improves over the fixed-filter baseline on all datasets, but on Weibo2014 the K=1 proxy gap is zero: every percentage point of gain comes from ensemble diversity, none from routing. Our full method breaks this degeneracy: \(\Delta_\text{K=1}\) is positive in 4 out of 5 folds on Weibo2014 (\(+0.029\)) and in 5 out of 5 folds on both BNCI datasets (\(+0.247\), \(+0.170\)). Routing entropy confirms the transition from degenerate (\(\bar{H}=0.836\)) to committed (\(\bar{H}=0.674\)) routing on Weibo2014, while domain alignment (\(r_d = 0.710\)) confirms the routing is structured by subject identity, not noise.
A secondary pattern is notable: \(\Delta_\text{K=1}\) is substantially larger in the low-\(\rho\) regime (\(+0.247\), \(+0.170\)) than in the high-\(\rho\) regime (\(+0.029\)). In the low-\(\rho\) regime the compact tangent space enables expert specialisation, making the single-expert ensemble filter a poor proxy and per-sample routing highly valuable.
LOSO generalisation results and directions for future work are presented in Appendices 10 and 11.
4.3 Ablation: What Breaks the Degeneracy?↩︎
| ID | Configuration | bAcc | \(\Delta_\text{base}\) | \(\Delta_\text{K=1}\) | \(\bar{H}\) | \(r_d\) | |
|---|---|---|---|---|---|---|---|
| Ref | Naive adaptive | 0.818 | +0.045 | \(-\)0.000 | 0.836 | 0.281 | |
| — Each fix in isolation — | lation —* | | ||||||
| B1 | \(+\) load-balance loss | 0.818 | +0.045 | +0.001 | 0.973 | 0.288 | |
| B2 | \(+\) entropy loss | 0.819 | +0.046 | +0.006 | 0.768 | 0.278 | |
| B4 | \(+\) alignment loss only | 0.818 | +0.045 | \(-\)0.000 | 0.853 | 0.386 | |
| B12\(^\dagger\) | \(+\) \(W_\text{base}\) only | 0.805 | +0.032 | +0.027 | 0.884 | 0.335 | |
| DSP\(^\dagger\) | \(+\) \(W_\text{base}\) + DSP | 0.810 | +0.037 | +0.016 | 0.776 | 0.690 | |
| — Full method — |
|
||||||
| DASP (Ours)\(^\dagger\) | \(W_\text{base}\) + DSP + align + decouple | 0.823 | +0.050 | +0.029 | 0.674 | 0.710 | |
We follow a strict one-factor-at-a-time protocol on Weibo2014. The results expose a clear causal chain (Table 3).
Load-balancing loss (B1) is a negative control. The Switch Transformer load-balancing loss [15] is minimised at uniform routing, pushing entropy toward 0.973. \(\Delta_\text{K=1}\) remains zero. This confirms our claim: high entropy causes degeneracy; we cannot resolve it by pushing entropy further up. The only path forward is to provide routing with information about what to differentiate — which requires structural, not just loss-level, changes.
\(W_\text{base}\) is the most impactful single change. Replacing the \(W_1\)-anchored reference point with a dedicated symmetric anchor immediately produces \(\Delta_\text{K=1} = +0.027\) (3 out of 5 folds), while domain alignment climbs from 0.281 to 0.335. The mechanism is geometric: symmetrising expert competition allows all \(K\) experts to receive equal gradient signal, enabling specialisation.
Frozen DSP amplifies the domain signal. Adding the DSP on top of \(W_\text{base}\) causes domain alignment to jump from 0.335 to 0.690, confirming that the projection provides the stable, domain-discriminative query signal needed for key convergence. Without it, the query MLP learns class-discriminative features that happen to be domain-informative only incidentally.
Alignment loss + decoupled keys seals the convergence. Alignment alone (B4, without \(W_\text{base}\)) raises domain alignment but not \(\Delta_\text{K=1}\), because the query provides a moving target. Combined with \(W_\text{base}\) and DSP, alignment produces the full method’s entropy of 0.674 and domain alignment of 0.710 — both conditions for domain-structured committed routing satisfied simultaneously.
4.4 Validation of the Scaling Methodology↩︎
BNCI2014001 (\(\rho = 28 < \rho^*\)) behaves identically to BNCI2015001 (\(\rho = 3.3\)) despite a \(7\times\) higher ratio: routing is domain-structured without DSP (\(r_d = 0.828\)), and adding DSP hurts rather than helps, confirming that the scaling rule generalises to intermediate channel counts. The contrast with Weibo2014 (\(\rho = 203\)) is sharp: without DSP, the query MLP cannot extract domain identity from the 1830-dimensional tangent space and routing collapses (\(\Delta_\text{K=1} \approx 0\)).
5 Full-Rank Rectangular Extension↩︎
The Stiefel constraint \(W_j \in \mathrm{St}(n,k)\) restricts each expert to the orthonormal group, preserving subspace orientation but precluding per-direction amplitude scaling. We describe here a theoretical generalisation to full-rank rectangular experts that we leave for empirical validation to future work (see §11).
The idea is to learn full-rank rectangular experts \(W_j \in \mathbb{R}^{n \times k}\), so that the Domain-Adaptive Pool layer \(Y = W_j^\top X W_j\) additionally encodes re-colouring transformation of the projected covariance. The space of full-rank rectangular matrices can be identified with the product of the Stiefel manifold \(\mathrm{St}(n,k)\) and the cone of symmetric positive definite matrices \(\mathrm{SPD}(k)\), making this a strict generalisation of the Stiefel prompt pool.
5.0.0.1 Implementation path 1 — Pseudo-polar decomposition.
Each expert is factored as \(W_j = Q_j P_j^{1/2}\), where \(Q_j \in \mathrm{St}(n,k)\) is an orthonormal factor maintained on the Stiefel manifold as before, and \(P_j \in \mathrm{SPD}(k)\) is a learnable symmetric positive definite matrix. The BiMap becomes \(Y = P_j^{1/2} Q_j^\top X Q_j P_j^{1/2}\), separating the rotation (\(Q_j\)) from the scaling (\(P_j\)) components. \(P_j\) can be parametrised via its Cholesky factor to guarantee positive definiteness throughout training [9], [16].
5.0.0.2 Implementation path 2 — Barrier regularisation.
Alternatively, each expert \(W_j \in \mathbb{R}^{n \times k}\) is learned as an unconstrained parameter, with a log-determinant barrier term added to the loss to enforce full column rank: \[\mathcal{L}_\text{barrier} = -\lambda_\text{barrier} \sum_{j=1}^K \log \det(W_j^\top W_j) \label{eq:barrier}\tag{11}\] As \(W_j\) approaches rank deficiency, \(\det(W_j^\top W_j) \to 0\) and \(\mathcal{L}_\text{barrier} \to +\infty\), creating an implicit manifold boundary that keeps experts full-rank without hard parametrisation. This approach requires no manifold-aware optimiser and is readily compatible with the existing training pipeline.
6 Related Work↩︎
6.0.0.1 Riemannian deep learning for EEG.
SPDNet [3] and extensions [17], [18] operate on SPD manifolds via BiMap, ReEig, and LogEig layers. Domain adaptation approaches include covariance recentering [19], Riemannian batch normalisation [17], and Riemannian Procrustes Analysis [20] — all requiring target-domain calibration data. Our method requires no test-time calibration and is, to the best of our knowledge, the first to achieve genuine input-adaptive routing on the SPD manifold.
6.0.0.2 Mixture of experts.
Sparse MoE [21] and Switch Transformer [15] use learned gating with load-balancing losses that encourage uniform utilisation — the opposite of our goal. We show explicitly (Table 3, B1) that applying such losses to our setting drives routing toward the degenerate solution.
6.0.0.3 Prompt-based continual learning.
L2P [4] and DualPrompt [22] select prompts via frozen pre-trained feature queries. Our work translates the key structural principles (frozen query, decoupled key training, alignment loss) from the Euclidean token setting to the Riemannian geometry of SPD manifolds.
7 Conclusion↩︎
We studied the problem of per-sample adaptive subspace selection in Riemannian deep learning for EEG decoding. Our central finding is that naive adaptive routing on the Stiefel manifold is provably degenerate: uniform weights reduce the layer to a single equal-experts-contribution filter, indistinguishable from the fixed-filter baseline in terms of per-sample adaptation. The K=1 proxy gap — the accuracy difference between the adaptive model and this degenerate solution — is the only diagnostic that separates genuine routing from ensemble averaging.
Three structural properties borrowed from Learning to Prompt break the degeneracy: a symmetric anchor, a frozen domain-discriminative query encoder, and a decoupled key alignment loss. Together they produce the first committed (\(\bar{H} = 0.67\)) and domain-structured (\(r_d = 0.71\)) routing on SPD manifolds, with consistent accuracy gains across three datasets of contrasting geometry. A unified data-driven scaling rule — governed by the tangent-to-domain ratio \(\rho\) — eliminates dataset-specific tuning.
8 Hyperparameter Summary↩︎
| Param | Weibo2014 | BNCI2015001 | BNCI2014001 | Description |
|---|---|---|---|---|
| \(n\) | 60 | 13 | 22 | Input SPD dimension |
| \(k\) | 30 | 12 | 20 | Output SPD dimension |
| \(K\) | 8 | 16 | 8 | Expert filters |
| \(m\) | 20 | 20 | 20 | Query/key/emb dim |
| \(r\) | 40 | 0 | 0 | DSP dim |
| \(\lambda_\text{align}\) | 0.05 | 0.0 | 0.0 | Alignment weight |
| decouple | True | False | False | Key gradient decoupling |
| lr | 0.01 | 0.01 | 0.01 | Adam [10] |
| Batch | 32 | 32 | 32 | |
| Epochs | 40 | 40 | 40 | Max, patience 10 |
9 Diagnostic Metric Definitions↩︎
9.0.0.1 Routing entropy.
\(\bar{H} = \frac{1}{N\log K}\sum_i H(\alpha_i)\), where \(H(\alpha_i) = -\sum_j \alpha_{i,j}\log\alpha_{i,j}\). Threshold \(\bar{H} < 0.70\) indicates sharp routing.
9.0.0.2 Domain alignment ratio.
\(r_d = \mathrm{Var}_d(\bar{\alpha}_d)/\mathrm{Var}_i(\alpha_i)\), where \(\bar{\alpha}_d = \frac{1}{|I_d|}\sum_{i:\,d_i=d}\alpha_i\) is the mean routing vector for domain \(d\), and \(|I_d|\) is the number of samples belonging to domain \(d\). Values above 0.10 indicate domain-structured routing.
9.0.0.3 Expert diversity.
Mean principal angle between column spaces of all expert pairs, computed via SVD of \(W_j^\top W_{j'}\) for all \(j \neq j'\). Values near \(\ang{45}\) indicate maximal diversity; near \(\ang{0}\) indicates collapse.
9.0.0.4 K=1 proxy.
The equal-expert contribution \(W_\text{base}^\dagger\) used as a single shared filter, evaluated on the test set. \(\Delta_\text{K=1} > 0.01\) in \(\geq 3/5\) folds is the success criterion for genuine routing beyond ensemble averaging.
10 LOSO Generalisation↩︎
10.0.0.1 Evaluation paradigm.
The main paper evaluates under a cross-domain protocol using a cross-validation procedure: all \(D\) domains (combination of subject-session) are seen during training, and test samples are drawn from the same domain pool. This protocol measures how well the routing mechanism adapts to known domains at test time, and is the standard evaluation in SPDNet-based motor imagery decoding [3].
A strictly harder paradigm is leave-one-subject-out (LOSO): each fold holds out one subject entirely from training, and the model must classify trials from a subject it has never seen. This is a genuine zero-shot transfer setting — the domain embedding \(\mathrm{Emb}(d_i)\) is undefined for the held-out subject, and the routing mechanism must operate without subject-specific supervision. LOSO therefore tests a fundamentally different capability: whether the routing generalises beyond the training domain distribution, rather than adapting within it.
10.0.0.2 Protocol.
We evaluate on Weibo2014 (\(D=9\), 8 LOSO folds, one subject held out per fold). Since the domain embedding cannot be used for the unseen subject, we compare four inference strategies:
Baseline BiMap: fixed filter, no routing.
Condition A — Input-only:
use_domain_emb=False; routing driven solely by the frozen between-domain encoder \(\tilde{x}_i\), with no domain embedding. This is the most principled zero-shot strategy since it uses only information available at test time.Condition B — Zero-mask: domain embedding zeroed for the held-out subject (\(\mathrm{Emb}(d_\text{new}) = \mathbf{0}\)).
Condition C — NN-inference: domain embedding initialised as the embedding of the nearest training subject, measured by L2 distance between subject-level tangent-space means.
| Condition | Description | bAcc | Std |
|---|---|---|---|
| Baseline BiMap | Fixed filter | 0.779 | 0.087 |
| Adaptive (full) | Domain emb for train subj. | 0.770 | 0.080 |
| A — Input-only | No domain emb | 0.776 | 0.074 |
| B — Zero-mask | Domain emb zeroed | 0.771 | 0.081 |
| C — NN-inference | Nearest training subject | 0.768 | 0.098 |
10.0.0.3 Discussion.
No adaptive variant consistently beats the baseline under LOSO, confirming that the method is designed for the cross-domain paradigm where subject identity is known at training time. The gap between the cross-domain result (\(0.823\)) and the best LOSO result (\(0.776\)) reflects the fundamental difficulty of zero-shot transfer: the routing mechanism learns domain attractors during training, and these attractors do not automatically generalise to unseen subjects.
The best strategy is input-only routing (Condition A, \(-0.003\) vs baseline, within variance), where the frozen between-domain encoder alone provides sufficient structure to make soft routing decisions without any domain supervision. This confirms that the encoder captures some subject-agnostic geometric structure in the tangent space, but not enough to match the performance achievable when the subject is known.
NN-based domain inference (Condition C) does not improve over zero-masking (\(0.768\) vs \(0.771\)), because L2 distance in tangent space is not a reliable proxy for routing similarity. A more geometrically principled alternative — computing the affine-invariant Riemannian distance \(d(X,Y) = \|\log(X^{-1/2}YX^{-1/2})\|_F\) between subject-level Fréchet means on \(\mathrm{SPD}(n)\) — is a natural direction for future work.
The tangent space projection at the identity \(I_n\) used in the query encoder is justified by the per-domain pre-conditioning applied to the covariance matrices (whitening and variance equalisation), which recenters each domain’s distribution around \(I_n\) before the network sees the data. Incorporating recentering directly within the model via a Riemannian batch normalisation layer [17], [23] is a natural extension that would remove the dependency on external pre-conditioning and remains an open direction.
11 Future Work↩︎
11.0.0.1 Top-\(p\) sparse routing.
Current routing computes a soft barycenter over all \(K\) experts for every sample, diluting gradient signals and making expert specialisation harder. Replacing the softmax over all \(K\) with a masked softmax over the top-\(p\) experts (straight-through estimator for the selection boundary) would create genuine winner-take-all pressure, analogous to L2P’s top-\(N\) prompt selection, and is expected to produce sharper routing without requiring additional losses.
11.0.0.2 Transductive LOSO adaptation.
The most promising near-term direction for improving LOSO performance is transductive inference: given a batch of unlabelled test samples from an unseen subject, optimise only the domain embedding \(\mathrm{Emb}(d_\text{new})\) for a small number of gradient steps on an unsupervised objective (e.g.routing entropy minimisation). All other parameters remain frozen. This requires no architectural changes and directly exploits the routing structure learned during training.
11.0.0.3 Riemannian NN-inference.
Replacing L2 tangent-space distance with the affine-invariant Riemannian distance \(d(X,Y) = \|\log(X^{-1/2}YX^{-1/2})\|_F\) between subject-level Fréchet means on \(\mathrm{SPD}(n)\) would give nearest-neighbour domain initialisation a more principled geometric basis, potentially improving Condition C in the LOSO setting.
11.0.0.4 Full-rank rectangular expert routing.
Section 5 introduces the full-rank rectangular generalisation of Stiefel expert routing and describes two concrete implementation paths; validating this extension empirically on non-pre-conditioned data, where between-subject amplitude differences are preserved, remains an open direction.
11.0.0.5 Per-domain Common Spatial Patterns (CSP) targets.
Maintaining a running exponential-moving average of per-domain per-class covariance estimates and computing domain-specific CSP targets would give routing supervision a stable per-domain attractor compatible with the alignment loss [24], addressing the instability observed with batch-level CSP supervision.