1 Introduction↩︎

Electroencephalography (EEG) measures multi-channel electric brain activity from the human scalp with millisecond precision [1]. Transient modulations in the rhythmic brain activity can reveal cognitive processes [2], affective states [3] and a person’s health status [4]. Unfortunately, these modulations exhibit low signal-to-noise ratio (SNR), domain shifts (i.e., changes in the data distribution) and have low specificity, rendering statistical learning a challenging task - particularly in the context of brain-computer interfaces (BCI) [5] where the goal is to predict a target from a short segment of multi-channel EEG data in real-time.

Under domain shifts, domain adaptation (DA), defined as learning a model from a source domain that performs well on a related target domain, offers principled statistical learning approaches with theoretical guarantees [6], [7]. DA in the BCI field mainly distinguishes inter-session and -subject transfer learning (TL) [8]. In inter-session TL, domain shifts, are expected across sessions mainly due to mental drifts (low specificity) as well as differences in the relative positioning of the electrodes and their impedances. Inter-subject TL is more difficult, as domain shifts are additionally driven by structural and functional differences in brain networks as well as variations in the performed task [9].

These domain shifts are traditionally circumvented by recording labeled calibration data and fitting domain-specific models [10], [11]. As recording calibration data is costly, models that are robust to scarce data with low SNR perform well in practice. Currently, tangent space mapping (TSM) models [12], [13] operating with symmetric, positive definite (SPD) covariance descriptors of preprocessed data are considered state-of-the-art (SoA) [10], [14], [15]. They are well suited for EEG data as they exhibit invariances to linear mixing of latent sources [16], and are consistent [13] and intrinsically interpretable [17] estimators for generative models that encode label information with a log-linear relationship in source power modulations.

Competitive, supervised calibration-free methods are one of the long-lasting grand challenges in EEG neurotechnology research [5], [10], [15], [18], [19]. Among the applied transfer learning techniques, including multi-task learning [20] and domain-invariant learning [21]–[23], unsupervised domain adaptation (UDA) [24] is considered as key to overcome this challenge [10], [19]. Contemporary methods cast the problem as a multi source and target UDA problem and address it with deep learning [25]–[28] or shallow, Riemannian geometry aware alignment methods [29]–[32]. Successful methods must cope with notoriously small and heterogeneous datasets (i.e., dozens of domains with a few dozens observations per domain and class). In a recent, relatively large scale inter-subject and -dataset TL competition with few labeled examples per target domain [19], deep learning approaches that aligned the first and second order statistics either in input [27], [33] or latent space [34] obtained the highest scores. Whereas, in a pure UDA competition [15] with a smaller dataset, Riemannian geometry aware approaches dominated. With the increasing popularity of geometric deep learning [35], [36] proposed an architecture based on SPD neural networks [37] to align SPD features in latent space and attained SoA scores. Despite the tremendous advances in recent years, the field still lacks methods that can consistently close the performance gap to state-of-the-art domain-specific methods.

To close this gap, we propose a machine learning framework around domain-specific batch normalization on the SPD manifold (Figure 1). The proposed framework is used to implement domain-specific TSM (Figure 1a), which requires tracking the domains’ Fréchet means in latent space as they are changing during training a typical TSM model in an end-to-end fashion (Figure 1b). After reviewing some preliminaries in section 2, we extend momentum batch normalization (MBN) [32] to SPDMBN that controls the Fréchet mean and variance of SPD data in section 3. In a theoretical analysis, we show under reasonable assumptions that SPDMBN can track and converge to the data’s true Fréchet mean, enabling, for the first time, end-to-end learning of feature extractors, TSM and tangent space classifiers. Building upon this insight, we combine SPDMBN with domain-specific batch normalization (DSBN) [38] to form SPDDSMBN (Figure 1a). A SPDDSMBN layer can transform domain-specific SPD inputs into domain-invariant SPD outputs (Figure 1c). Like DSBN, SPDDSMBN easily extends to multi-source, multi-target and online UDA scenarios. In section 4, we briefly review the generative model of EEG, before the proposed methods are combined in a simple, intrinsically interpretable network architecture, denoted TSMNet (Figure 1b). We obtain state-of-the-art performance in inter-session and -subject UDA on small and large scale EEG BCI datasets, and show in an ablation study that the performance increase is primarily driven by performing DSBN on the SPD manifold.

2 Preliminaries↩︎

2.0.0.1 Multi-source multi-target unsupervised domain adaptation

Let \(\mathcal{X}\) denote the space of input features, \(\mathcal{Y}\) a label space, and \(\mathcal{I}_d \subset \mathbb{N}\) an index set that contains unique domain identifiers. In the multi-source, multi-target unsupervised domain adaptation scenario considered here, we are given a set \(\mathcal{T}^{source} = \{\mathcal{T}_i | i \in \mathcal{I}_d^{source} \subset \mathcal{I}_d\}\) with \(|\mathcal{I}_d^{source}|=N\) domains. Each domain \(\mathcal{T}_i = \{(\mathbf{X}_{ij},y_{ij})\}_{j=1}^{M} \sim P_{XY}^i\) contains \(M\) observations of feature (\(\mathbf{X} \in \mathcal{X}\)) and label (\(y \in \mathcal{Y}\)) tuples sampled from a joint distribution \(P_{XY}^i\).¹ While the joint distributions can be different (but related) across domains, we assume that the class priors are the same (i.e., \(P_{Y}^i = P_{Y}\)). The goal is to learn a predictive function \(h : \mathcal{X}\times \mathcal{I}_d \rightarrow \mathcal{Y}\) that, once fitted to \(\mathcal{T}^{source}\), can generalize to unseen target domains \(\mathcal{T}^{target} = \{\mathcal{T}_l | l \in \mathcal{I}_d^{target} \subset \mathcal{I}_d, \mathcal{I}_d^{target} \cap \mathcal{I}_d^{source} = \emptyset\}\) merely based on unsupervised adaptation of \(h\) to each target domain \(\mathcal{T}_l\) once its label \(l\) and features \(\{\mathbf{X}_{lj}\}_{j=1}^{M} \sim P_{X}^l\) are revealed.

2.0.0.2 Riemannian geometry on \(\mathcal{S}^+_D\)

We start with recalling notions of geometry on the space of real \(D \times D\) symmetric positive definite (SPD) matrices \(\mathcal{S}^+_D = \{ \mathbf{Z} \in \mathbb{R}^{D \times D} : \mathbf{Z}^T = \mathbf{Z}, \mathbf{Z} \succ 0 \}\). The space \(\mathcal{S}^+_D\) forms a cone shaped Riemannian manifold in \(\mathbb{R}^{D\times D}\) [39]. A Riemannian manifold \(\mathcal{M}\) is a smooth manifold equipped with an inner product on the tangent space \(\mathcal{T}_\mathbf{Z}\mathcal{M}\) at each point \(\mathbf{Z} \in \mathcal{M}\). Tangent spaces have Euclidean structure with easy to compute distances \(\mathcal{T}_\mathbf{Z}\mathcal{M} \times \mathcal{T}_\mathbf{Z}\mathcal{M} \rightarrow \mathbb{R}^+\) which locally approximate Riemannian distances on \(\mathcal{M}\) induced by an inner product [40]. Logarithmic \(\mathrm{Log}_\mathbf{Z}:\mathcal{M} \rightarrow \mathcal{T}_\mathbf{Z}\mathcal{M}\) and exponential \(\mathrm{Exp}_\mathbf{Z}: \mathcal{T}_\mathbf{Z}\mathcal{M} \rightarrow \mathcal{M}\) mappings project points to and from tangent spaces.

Using the inner product \(\langle\mathbf{S}_1,\mathbf{S}_2\rangle_\mathbf{Z} = \mathrm{Tr}(\mathbf{Z}^{-1} \mathbf{S}_1 \mathbf{Z}^{-1} \mathbf{S}_2 )\) for points \(\mathbf{S}_1,\mathbf{S}_2\) in the tangent space \(\mathcal{T}_\mathbf{Z}\mathcal{S}^+_D\) (i.e., the space of real symmetric \(D \times D\) matrices) results in a globally defined affine invariant Riemannian metric on \(\mathcal{S}^+_D\) [39], [41], which can be computed in closed form: \[\delta(\mathbf{Z}_1,\mathbf{Z}_2) = ||\mathrm{log}(\mathbf{Z}_1^{-\frac{1}{2}}\mathbf{Z}_2\mathbf{Z}_1^{-\frac{1}{2}})||_F \label{eq:airm}\tag{1}\] where \(\mathbf{Z}_1\) and \(\mathbf{Z}_2\) are two SPD matrices, \(\mathrm{log}(\cdot)\) denotes the matrix logarithm², \(||\cdot||_F\) the Frobenius norm, and \(\mathrm{Tr}(\cdot)\) in the inner product the trace operator. Due to affine invariance, we have \(\delta_{}(\mathbf{A} \mathbf{Z}_1 \mathbf{A}^T, \mathbf{A}\mathbf{Z}_2 \mathbf{A}^T) = \delta_{}(\mathbf{Z}_1,\mathbf{Z}_2)\) for any invertible \(D \times D\) transformation matrix \(\mathbf{A}\). The exponential and logarithmic mapping are also globally defined in closed form as \[\mathrm{Log}_\mathbf{Z}(\mathbf{Z}_1) = \mathbf{Z}^{\frac{1}{2}} \mathrm{log}(\mathbf{Z}^{-\frac{1}{2}}\mathbf{Z}_1\mathbf{Z}^{-\frac{1}{2}}) \mathbf{Z}^{\frac{1}{2}} \label{eq:log}\tag{2}\] \[\mathrm{Exp}_\mathbf{Z}(\mathbf{S}_1) = \mathbf{Z}^{\frac{1}{2}} \mathrm{exp}(\mathbf{Z}^{-\frac{1}{2}}\mathbf{S}_1\mathbf{Z}^{-\frac{1}{2}}) \mathbf{Z}^{\frac{1}{2}} \label{eq:exp}\tag{3}\]

For a set of SPD points \(\mathcal{Z} = \{\mathbf{Z}_j \in \mathcal{S}^+_D\}_{j \le M}\), we will use the notion of Fréchet mean \(\mathbf{G}_{\mathcal{Z}} \in \mathcal{S}^+_D\) and Fréchet variance \(\nu^2_{\mathcal{Z}} \in \mathbb{R}^+\). The Fréchet mean is defined as the minimizer of the average squared distances \[\mathbf{G}_{\mathcal{Z}} = \mathrm{arg \min_{\mathbf{G}\in \mathcal{S}^+_D}} \frac{1}{M}\sum_{j=1}^{M} \delta_{}^2(\mathbf{G},\mathbf{Z}_j) \label{eq:geommean}\tag{4}\] For \(M=2\), there is a closed form solution expressed as \[\mathbf{G}_{\mathcal{Z}}(\gamma) = \mathbf{Z}_1 \#_\gamma \mathbf{Z}_2 = \mathbf{Z}_1^{\frac{1}{2}} \left(\mathbf{Z}_1^{-\frac{1}{2}}\mathbf{Z}_2\mathbf{Z}_1^{-\frac{1}{2}}\right)^\gamma \mathbf{Z}_1^{\frac{1}{2}} \label{eq:mean2points}\tag{5}\] with weight \(\gamma = 0.5\). Choosing \(\gamma \in [0,1]\) computes weighted means along the geodesic (i.e., the shortest curve) that connects both points. For \(M > 2\), (4 ) can be solved using the Karcher flow algorithm [42], which iterates between projecting the data to the tangent space (2 ) at the current estimate, arithmetic averaging, and projecting the result back (3 ) to obtain a new estimate. The Fréchet variance \(\nu_{\mathcal{Z}}^2\) is defined as the attained value at the minimizer \(\mathbf{G}_{\mathcal{Z}}\): \[\nu^2_{\mathcal{Z}} = \mathrm{Var}_{\mathcal{Z}}(\mathbf{G}_{\mathcal{Z}}) = \frac{1}{M}\sum_{j=1}^{M} \delta_{}^2(\mathbf{G}_{\mathcal{Z}},\mathbf{Z}_j) \label{eq:geomvar}\tag{6}\]

To shift a set of tangent space points to vary around a parametrized mean \(\mathbf{G}_\phi\), parallel transport on \(\mathcal{S}^+_D\) can be used [43]: \[\Gamma_{\mathbf{G}_{\mathcal{Z}}\rightarrow\mathbf{G}_\phi}(\mathbf{S}) = \mathbf{E}^T\mathbf{S}\mathbf{E}~,~~\mathbf{E} = (\mathbf{G}_{\mathcal{Z}}^{-1}\mathbf{G}_\phi)^{\frac{1}{2}} \label{eq:paralleltransport}\tag{7}\] While, parallel transport is generally defined for tangent space vectors \(\mathbf{S}\) [40], on \(\mathcal{S}^+_D\) the same operations also apply directly to points on the manifold (i.e., \(\mathbf{Z} \in \mathcal{Z}\)) [31], [44].

3 Domain-specific batch normalization on \(\mathcal{S}^+_D\)↩︎

In this section, we review relevant batch normalization (BN) [45] variants with a focus on \(\mathcal{S}^+_D\). We then present SPDMBN and show in a theoretical analysis that the running estimate converges to the true Fréchet mean under reasonable assumptions. At last, we combine the idea of domain-specific batch normalization (DSBN) [38] with SPDMBN to form a SPDDSMBN layer. Table 1 provides a brief overview of related and proposed methods.

3.0.0.1 Batch normalization

Batch normalization (BN) [45] is a widely used training technqiue in deep learning as BN layers speed up convergence and improve generalization via smoothing of the engery landscape [32], [47]. A standard BN layer applies slightly different transformations during training and testing to independent and identically distributed (iid) observations \(\mathbf{x}_j \in \mathbb{R}^d\) within the \(k\)-th minibatch \(\mathcal{B}_k\) of size \(M\) drawn from a dataset \(\mathcal{T}\). During training, the data are normalized using the batch mean \(\mathbf{b}_k\) and variance \(\mathbf{s}^2_k\), and then scaled and shifted to have a parametrized mean \(\mathbf{g}_\phi\) and variance \(\boldsymbol{\sigma}^2_\phi\). Internally, the layer updates running estimates of the dataset’s statistics (\(\mathbf{g}_k, \boldsymbol{\sigma}^2_k\)) during each training step \(k\); the updates are computed via exponential smoothing with momentum parameter \(\gamma\). During testing, the running estimates are used.

Using batch statistics to normalize data during training rather than running estimates introduces noise whose level depends on the batch size [32]; smaller batch sizes raise the noise level. The introduced noise regularizes the training process, which can help to escape poor local minima in initial learning but also lead to underfitting. Momentum BN (MBN) [32] allows small batch sizes while avoiding underfitting. Like batch renormalization [48], MBN uses running estimates during training and testing. The key difference is that MBN keeps two sets of running statistics; one for training and one for testing. The latter are updated conventionally, while the former are updated with momentum parameter \(\gamma_{train}(k)\) that decays over training steps \(k\). MBN can, therefore, quickly escape poor local minima during initial learning and avoid underfitting at later stages [32].

3.0.0.2 Batch normalization on \(\mathcal{S}^+_D\)

It is intractable to compute the Fréchet mean \(\mathbf{G}_{\mathcal{B}_k}\) for each minibatch \(\mathcal{B}_k = \{\mathbf{Z}_j \in \mathcal{S}^+_D\}_{j=1}^M\), as there is no efficient algorithm to solve (4 ). [44] proposed Riemannian Batch Normalization (RBN) as a tractable approximation. RBN approximateley solves (4 ) by aborting the iterative Karcher flow algorithm after one iteration. To transform \(\mathbf{Z}_j \in \mathcal{B}_k\) with estimated mean \(\mathbf{B}_k\) to vary around \(\mathbf{G}_\phi\), parallel transport (7 ) is used. The RBN input output transformation is then expressed as \[\mathrm{RBN}(\mathbf{Z}_j; \mathbf{G}_\phi, \gamma) = \Gamma_{\mathbf{B}_k\rightarrow\mathbf{G}_\phi}(\mathbf{Z}_j) = \mathbf{E}^T\mathbf{Z}_j\mathbf{E}~,~~\mathbf{E} = (\mathbf{B}_k^{-1}\mathbf{G}_\phi)^{\frac{1}{2}} ~, ~~~ \forall \mathbf{Z}_j\in\mathcal{B}_k \label{eq:rbn}\tag{8}\] Using (5 ), the running estimate of the dataset’s Fréchet mean can be updated in closed form \[\mathbf{G}_k = \mathbf{G}_{k-1} \#_{\gamma} \mathbf{B}_k \label{eq:spdrunningmean}\tag{9}\] In [46] we proposed an extension to RBN, denoted SPD batch renormalization (SPDBN) that controls both Fréchet mean and variance. Like batch renormalization [48], SPDBN uses running estimates \(\mathbf{G}_k\) and \(\nu^2_k\) during training and testing. To transform \(\mathbf{Z}_j \in \mathcal{B}_k\) to vary around \(\mathbf{G}_\phi\) with variance \(\nu^2_\phi\), each observation is first transported to vary around the identity matrix \(\mathbf{I}\), rescaled via computing matrix powers and finally transported to vary around \(\mathbf{G}_\phi\). The sequence of operations can be expressed as \[\mathrm{SPDBN}(\mathbf{Z}_j; \mathbf{G}_\phi, \nu^2_\phi, \varepsilon, \gamma) = \Gamma_{\mathbf{I}\rightarrow\mathbf{G}_\phi} \circ \Gamma_{\mathbf{G}_k\rightarrow\mathbf{I}}(\mathbf{Z}_j)^{\frac{\nu_\phi}{\nu_k + \varepsilon}}~, ~~~ \forall \mathbf{Z}_j\in\mathcal{B}_k \label{eq:spdbatchnorm}\tag{10}\] The standard backpropagation framework with extensions for structured matrices [49] and manifold-constrained gradients [40] can be used to propagate gradients through RBN and SPDBN layers and learn the parameters (\(\mathbf{G}_\phi, \nu_\phi\)).

3.0.0.3 Momentum batch normalization on \(\mathcal{S}^+_D\)

SPDBN [46] suffers from the same limitations as batch renormalization [48]. Consequently, we propose to extend MBN [32] to \(\mathcal{S}^+_D\). We list the pseudocode of our proposed extension, which we denote SPDMBN, in algorithm 2. SPDMBN uses approximations of batch-specific Fréchet means to update two sets of running estimates of the dataset’s Fréchet mean. As MBN [32], we decay \(\gamma_{train}(k)\) with a clamped exponential decay schedule \[\gamma_{train}(k) = 1 - \gamma_{min}^{\frac{1}{K-1} \mathrm{max}(K-k,0)} + \gamma_{min} \label{eq:momentumschedule}\tag{11}\] where \(K\) defines the training step at which \(\gamma_{min} \in [0,1]\) should be attained.

3.0.0.4 The running mean in SPDMBN converges to the Fréchet mean

Here, we consider models that apply a SPDMBN layer to latent representations generated by a feature extractor \(f_\theta : \mathcal{X} \rightarrow \mathcal{S}^+_D\) with learnable parameters \(\theta\).

We define a dataset that contains the latent representations generated with feature set \(\theta_k\) as \(\mathcal{Z}_{\theta_k} = \{ f_{\theta_k} (\mathbf{x}) | \mathbf{x} \in \mathcal{T} \}\), and a minibatch of \(M\) iid samples at training step \(k\) as \(\mathcal{B}_k\). We denote the Fréchet mean of \(\mathcal{Z}_{\theta_k}\) as \(\mathbf{G}_{\theta_k}\), the estimated Fréchet mean, defined in (9 ), as \(\mathbf{G}_k\), and the estimated batch mean as \(\mathbf{B}_k\). Since the batches are drawn randomly, we consider the batch and running means as random variables.

We assume that the variance \(\mathrm{Var}_{\theta_{k}}(\mathbf{B}_{k-1}) = \mathbb{E}_{\mathbf{B}_{k-1}} \{ \delta^2(\mathbf{B}_{k-1}, \mathbf{G}_{\theta_{k}}) \}\) of the previous batch mean \(\mathbf{B}_{k-1}\) with respect to the current Fréchet mean \(\mathbf{G}_{\theta_k}\) is bounded by the current variance \(\mathrm{Var}_{\theta_{k}}(\mathbf{B}_{k})\) and the norm of the difference in the parameters \[\mathrm{Var}_{\theta_k}(\mathbf{B}_{k-1}) \le (1 + || \theta_k - \theta_{k-1} ||) \mathrm{Var}_{\theta_k}(\mathbf{B}_{k}) \label{eq:assumption1}\tag{12}\] That is, across training steps \(k\) the parameter updates are required to change the first and second order moments of the distribution of \(\mathcal{Z}_{\theta_k}\) gradually so that the expected distance between \(\mathbf{G}_{\theta_k}\) and the previous batch mean \(\mathbf{B}_{k-1}\) is bounded. We conjecture that this is the case for feature extractors \(f_\theta\) that are smooth in the parameters and small learning rates, but leave the proof for future work.

The proof is provided in appendix 7.1 of the supplementary material and relies on the proof of the geometric law of large numbers [50].

Proposition 1 states that if (12 ) and (?? ) are met, the expected distance between the true Fréchet mean and the running mean is less or equal to the one of the batch mean. Consequently, the introduced noise level of SPDBN (equation 10 ) and SPDMBN (algorithm 2), which use \(\mathbf{G}_k\) to normalize batches during training, is smaller or equal to RBN (equation 8 ), which uses \(\mathbf{B}_k\).

Since \(\gamma\) controls the adaptation speed of \(\mathbf{G}_k\), proposition 1 also states that if \(\gamma\) converges to zero (=no adaptation), the parameter updates are required to converge to zero as well (=no learning). Hence, for a fixed \(\gamma \in (0,1)\), as in the case of SPDBN (equation 10 ), proposition 1 is fulfilled, if the learning rate for the parameters \(\theta\) is chosen sufficiently small. This can substantially slow down initial learning for standard choices of \(\gamma\) (e.g., 0.1 or 0.01). As a remedy, SPDMBN (algorithm 2) uses an adaptive momentum parameter, which allows larger parameter updates during initial training steps.

If we consider a late stage of learning, and in particular assume that after a certain number of iterations \(\kappa\) the parameters stay in a small ball with radius \(\rho\) around \(\theta^*\) (i.e., \(|| \theta_k - \theta^* || \le \rho ~~~\forall~k > \kappa\)) and the feature extractor is \(L\)-smooth in the parameters (i.e., \(\delta(f_\theta(\mathbf{x}), f_{\tilde{\theta}}(\mathbf{x})) \le L || \theta - \tilde{\theta} || ~ \forall \mathbf{x} \in \mathcal{T} ~,~ \forall \theta,\tilde{\theta}\)) then the distances are bounded \(\delta(f_{\theta_k}(\mathbf{x}), f_{\theta^*}(\mathbf{x})) \le \rho L\).

Taken together, Proposition 1 and Remark 1 provide guidelines to update \(\mathbf{G}_k\) in SPDMBN so that the introduced estimation error is bounded during initial fast learning (large \(\gamma\)) and decays towards zero in late learning (small \(\gamma\)).

3.0.0.5 SPDMBN to learn tangent space mapping at Fréchet means

Typical TSM models for classification [12] and regression [13] first use (2 ) to project \(\mathbf{Z} \in \mathcal{T} \subset \mathcal{S}_D^+\) to the tangent space at the Fréchet mean \(\mathbf{G}_\mathcal{T}\), then use (7 ) to transport the result to vary around \(\mathbf{I}\), and finally extract elements in the upper triangular part³ to reduce feature redundancy. The invertible mapping \(\mathcal{P}_{\mathbf{G}_\mathcal{T}} : \mathcal{S}_D^+ \rightarrow \mathbb{R}^{D(D+1)/2}\) is expressed as: \[\mathcal{P}_{\mathbf{G}_\mathcal{T}}(\mathbf{Z}) = \mathrm{upper} \circ \Gamma_{\mathbf{G}_\mathcal{T}\rightarrow\mathbf{I}} \circ \mathrm{Log}_{\mathbf{G}_\mathcal{T}}(\mathbf{Z}) = \mathrm{upper}(\mathrm{log}(\mathbf{G}_\mathcal{T}^{-\frac{1}{2}}\mathbf{Z}\mathbf{G}_\mathcal{T}^{-\frac{1}{2}})) \label{eq:tsm}\tag{13}\]

We propose to use a SPDMBN layer followed by a LogEig layer [37] to compute a similar mapping \(m_\phi\) (Figure 1a). A LogEig layer simply computes the matrix logarithm and vectorizes the result so that the norm is preserved. If the parametrized mean of SPDMBN is fixed to the identify matrix (\(\mathbf{G}_\phi= \mathbf{I}\)), the composition computes \[\begin{align} m_\phi(\mathbf{Z}) & = \mathrm{LogEig} \circ \mathrm{SPDMBN} (\mathbf{Z}) = \mathrm{upper} \circ \mathrm{log} \circ \Gamma_{\mathbf{G}_k\rightarrow\mathbf{I}}(\mathbf{Z})^{\frac{\nu_\phi}{\nu_k + \varepsilon}} \nonumber \\ & = \mathrm{upper} \left( \frac{\nu_\phi}{\nu_k + \varepsilon} \mathrm{log} \left( \mathbf{G}_k^{-\frac{1}{2}}\mathbf{Z}\mathbf{G}_k^{-\frac{1}{2}} \right)\right) \label{eq:tsmlayer} \end{align}\tag{14}\] where \((\mathbf{G}_k, \nu^2_k)\) are the estimated Fréchet mean and variance of the dataset \(\mathcal{T}\) at training step \(k\), and \(\phi = \{ \nu_\phi \}\) the learnable parameters. According to remark 1 \(\mathbf{G}_k\) converges to \(\mathbf{G}_{\mathcal{T}}\) and, in turn, \(m_\phi\) to a scaled version of \(\mathcal{P}_{\mathbf{G}_\mathcal{T}}\), since \(\mathrm{upper}\) is linear.

The mapping \(m_\phi\) offers several advantageous properties. First, the features are projected to a Euclidean vector space where standard layers can be applied and distances are cheap to compute. Second, distances between the projected features locally approximate \(\delta_{}\) and, therefore, inherit its invariance properties (e.g., affine mixing) [41]. This improves upon a LogEig layer [37] which projects features to the tangent space at the identity matrix. As a result, distances between LogEig projected features correspond to distances measured with the log-Euclidean Riemannian metric (LERM) [51] which is not invariant to affine mixing. Third, controlling the Fréchet variance in (14 ) empirically speeds up learning and improves generalization [46].

3.0.0.6 Domain-specific batch normalization on \(\mathcal{S}^+_D\)

Considering a multi-source UDA scenario, [38] proposed a domain-specific BN (DSBN) layer which simply keeps multiple parallel BN layers and distributes observations according to the associated domains. Formally, we consider minibatches \(\mathcal{B}_k\) that form the union of \(N_{\mathcal{B}_k} \le |\mathcal{I}_d|\) domain-specific minibatches \(\mathcal{B}_k^i\) drawn from distinct domains \(i \in \mathcal{I}_{\mathcal{B}_k} \subseteq \mathcal{I}_d\). As before, each \(\mathcal{B}_k^i\) contains \(j=1,...,M/N_{\mathcal{B}_k}\) iid observations \(\mathbf{x}_{j}\). A DSBN layer mapping \(\mathbb{R}^d\times \mathcal{I}_d \rightarrow \mathbb{R}^d\) can then be expressed as \[\mathrm{DSBN}(\mathbf{x}_j, i) = \mathrm{BN}_i(\mathbf{x}_j; \mathbf{g}_{\phi_i}, \mathbf{s}_{\phi_i}, \varepsilon, \gamma) ~, ~~~ \forall \mathbf{x}_j\in\mathcal{B}_k^i ~,~~~\forall i \in \mathcal{I}_{\mathcal{B}_k} \label{eq:dsbn}\tag{15}\] In practice, the batch size \(M\) is typically fixed. The particular choice is influenced by resource availability and the desired noise level introduced by minibatch based stochastic gradient descent. A drawback of DSBN is that for a fixed batch size \(M\) and an increasing number of source domains \(N_{\mathcal{B}_k}\), the effective batch size declines for the BN layers within DSBN. Since small batch sizes increase the noise level introduced by BN, increasing the number of domains per batch can lead to underfitting [32]. To alleviate this effect, we use the previously introduced SPDMBN layer. The proposed domain-specific BN layer on \(\mathcal{S}_D^+\) is then formally defined as \[\mathrm{SPDDSMBN}(\mathbf{Z}_j, i) = \mathrm{SPDMBN}_i(\mathbf{Z}_j; \mathbf{G}_{\phi_i}, \nu_{\phi_i}, \varepsilon, \gamma, \gamma_{train}(k)) ~, ~ \forall \mathbf{Z}_j\in\mathcal{B}_k^i \subset \mathcal{S}^+_D ~,~\forall i \in \mathcal{I}_{\mathcal{B}_k} \label{eq:spddsmbn}\tag{16}\] The layer can be readily adapted to new domains, as new SPDMBN layers can be added on the fly. If the entire data of a domain becomes available, the domain-specific Fréchet mean and variance can be estimated by solving (4 ), otherwise, the update rules in algorithm 2 can be used.

4 SPDDSMBN to crack interpretable multi-source/-target UDA for EEG data↩︎

With SPDDSMBN introduced in the previous section, we focus on a specific application domain, namely, multi-source/-target UDA for EEG-based BCIs and propose an intrinsically interpretable architecture which we denote TSMNet.

4.0.0.1 Generative model of EEG

EEG signals \(\mathbf{x}(t) \in \mathbb{R}^P\) capture voltage fluctuations on \(P\) channels. An EEG record (=domain) is uniquely identified by a subject and session identifier. After standard pre-processing steps, each domain \(i\) contains \(j = 1,...,M\) labeled observations with features \(\mathbf{X}_{ij} \in \mathcal{X} \subset \mathbb{R}^{P\times T}\) where \(T\) is the number of temporal samples. Due to linearity of Maxwell’s equations and Ohmic conductivity of tissue layers in the frequency ranges relevant for EEG [52], a domain-specific linear instantaneous mixture of sources model is a valid generative model: \[\mathbf{X}_{ij} = \mathbf{A}_{i} \mathbf{S}_{ij} + \mathbf{N}_{ij}\] where \(\mathbf{S}_{ij} \in \mathbb{R}^{Q\times T}\) represents the activity of \(Q\) latent sources, \(\mathbf{A}_{i} \in \mathbb{R}^{P\times Q}\) a domain-specific mixing matrix and \(\mathbf{N}_{ij} \in \mathbb{R}^{P\times T}\) additive noise. Both \(\mathbf{A}_{i}\) and \(\mathbf{S}_{ij}\) are unknown which demands making assumptions on \(\mathbf{A}_{i}\) (e.g., anatomical prior knowledge [53]) and/or \(\mathbf{S}_{ij}\) (e.g., statistical independence [54]) to extract interesting sources.

4.0.0.2 Interpretable multi-source/-target UDA for EEG data

As label information is available for the source domains, our goal is to identify discriminative oscillatory sources shared across domains. Our approach relies on TSM models with linear classifiers [12], as they are consistent [13] and intrinsically interpretable [17] estimators for generative models with log-linear relationships between the target \(y_{ij}\) and variance \(\mathrm{Var}\{s_{ij}^{(k)}(t)\}\) of \(k=1,...,K \le Q\) discriminative sources: \[y_{ij} = \sum_{k=1}^K b_k \mathrm{log}\left( \mathrm{Var}\{s_{ij}^{(k)}(t)\} \right)+ \varepsilon_{ij} \label{eq:encoding95model}\tag{17}\] where \(b_k \in \mathbb{R}\) summarizes the coupling between the target \(y_{ij}\) and the variance of the encoding source, and \(\varepsilon_{ij}\) additive noise. In [17] we showed that the encoding sources’ coupling and their patterns⁴ (columns of \(\mathbf{A}_{i}\)) can be recovered via solving a generalized eigenvalue problem between the Fréchet mean \(\mathbf{G}_{\mathcal{T}_i}\) and classifier patterns [55] that were back projected to \(\mathcal{S}_D^+\) with \(\mathcal{P}_{\mathbf{G}_{\mathcal{T}_i}}^{-1}\). The resulting eigenvectors are the patterns and the eigenvalues \(\lambda_k\) reflect the relative source contribution \(c_k\): \[c_k = \mathrm{max}(\lambda_k,\lambda_k^{-1})~,~~~\lambda_k = \mathrm{exp}(b_k/||\mathbf{b}||_2^2) \label{eq:pattern95contribution}\tag{18}\] To benefit from the intrinsic interpretability of TSM models, we constrain our hypothesis class \(\mathcal{H}\) to functions \(h : \mathcal{X}\times \mathcal{I}_d \rightarrow \mathcal{Y}\) that can be decomposed into a composition of a shared linear feature extractor with covariance pooling \(f_{\theta} : \mathcal{X} \rightarrow \mathcal{S}_D^+\), domain-specific tangent space mapping \(m_{\phi} : \mathcal{S}_D^+ \times \mathcal{I}_d \rightarrow \mathbb{R}^{D(D+1)/2}\), and a shared linear classifier \(g_{\psi}: \mathbb{R}^{D(D+1)/2} \rightarrow \mathcal{Y}\) with parameters \(\Theta = \{\theta, \phi, \psi\}\).

4.0.0.3 TSMNet with SPDDSMBN

Unlike previous approaches which learn \(f_{\theta}, m_{\phi}, g_{\psi}\) sequentially [13], [29]–[31], we parametrize \(h = g_{\psi} \circ m_{\phi} \circ f_{\theta}\) as a neural network and learn the entire model in an end-to-end fashion (Figure 1b). Details of the proposed architecture, denoted TSMNet, are provided in appendix 7.2.5. In a nutshell, we parametrize \(f_{\theta}\) as the composition of the first two linear convolutional layers of ShConvNet [56], covariance pooling [57], BiMap [37], and ReEig [37] layers. A BiMap layer applies a linear subspace projection, and a ReEig layer thresholds eigenvalues of symmetric matrices so that the output is SPD. We used the default threshold (\(10^{-4}\)) and found that it was never active in the trained models. Hence, after training, \(f_{\theta}\) fulfilled the hypothesis class constraints. In order for \(m_{\phi}\) to align the domain data and compute TSM, we use SPDDSMBN (16 ) with shared parameters (i.e., \(\mathbf{G}_{\phi_i} = \mathbf{G}_{\phi} = \mathbf{I}, \nu_{\phi_i} = \nu_{\phi}\)) in (14 ). Finally, the classifier \(g_{\psi}\) was parametrized as a linear layer with softmax activations. We use the standard-cross entropy loss as training objective, and optimized the parameters with the Riemannian ADAM optimizer [58].

5 Experiments with EEG data↩︎

In the following, we apply our method to classify target labels from short segments of EEG data. We consider two BCI applications, namely, mental imagery [5], [59] and mental workload estimation [60]. Both applications have high potential to aid society in rehabilitation and healthcare [18], [61], [62] but have, currently, limited practical value because of poor generalization across sessions and subjects [15], [19].

5.0.0.1 Datasets and preprocessing

The considered mental imagery datasets were BNCI2014001 [63] (9 subjects/2 sessions/4 classes), BNCI2015001 [64] (12/2-3/2), Lee2019 [65] (54/2/2), Lehner2020 [66] (1/7/2), Stieger2021 [67] (62/4-8/4) and Hehnberger2021 [68] (1/26/4). For mental workload estimation, we used a recent competition dataset [69] (12/2/3). A detailed description of the datasets is provided in appendix 7.2.1. Altogether, we analyzed a total of 603 sessions of 158 human subjects whose data was acquired in previous studies that obtained the subjects’ informed consent and the right to share anonymized data.
The python packages moabb [14] and mne [70] were used to preprocess the datasets. The applied steps comprise resampling the EEG signals to 250/256 Hz, applying temporal filters to extract oscillatory EEG activity in the 4 to 36 Hz range (spectrally resolved if required by a method) and finally extract short segments (\(\le 3 s\)) associated to a class label (details provided in appendix 7.2.2).

5.0.0.2 Evaluation

We evaluated TSMNet against several baseline methods implementing direct transfer or multi-source (-target) UDA strategies. They can be broadly categorized as component based [68], [71], Riemannian geometry aware [12], [17], [30], [72] or deep learning [25], [56], [73]. All models were fit and evaluated with a randomized leave 5% of the sessions (inter-session TL) or subjects (inter-subject TL) out cross-validation (CV) scheme. For inter-session TL, the models were only provided with data of the associated subject. When required, inner train/test splits (neural nets) or CV (shallow methods) were used to optimize hyper parameters (e.g., early stopping, regularization parameters). The dataset Hehenberger2021 was used to fit the hyper parameters of TSMNet, and is, therefore, omitted in the presented results. Balanced accuracy (i.e., the average recall across classes) was used as scoring metric. As the discriminability of the data varies considerably across subjects, we decided to report the results in the figures relative to the score of a SoA domain-specific Riemannian geometry aware method [17], which was fitted and evaluated in a 80%/20% chronological train/test split (for details see appendix 7.2.4).

5.0.0.3 Soft- and hardware

We either used publicly available python code or implemented the methods in python using the packages torch [74], scikit-learn [75], skorch [76], geoopt [77], mne [70], pyriemann [78], pymanopt [79]. We ran the experiments on standard computation PCs equipped with 32 core CPUs with 128 GB of RAM and used up to 1 GPU (24 GRAM). Depending on the dataset size, fitting and evaluating TSMNet varied from a few seconds to minutes.

5.1 Mental imagery↩︎

5.1.0.1 TSMNet closes the gap to domain-specific methods

Figure 3 summarizes the mental imagery results. It displays test set scores of the considered TL methods relative to the score of a SoA domain-specific reference method. Combining the results of all subjects across datasets (Figure 3a), it becomes apparent that TSMNet is the only method that can significantly reduce the gap to the reference method (inter-subject) or even exceed its performance (inter-session). Figure 3b displays the results resolved across datasets (for details see appendix 7.3.1). We make two important observations. First, concerning inter-session TL, TSMNet meets or exceeds the score of the reference method consistently across datasets. Second, concerning inter-subject TL, we found that all considered methods tend to reduce the performance gap as the dataset size (# subjects) increases, and that TSMNet is consistently the top or among the top methods. As a fitted TSMNet corresponds to a typical TSM model with a linear classifier, we can transform the fitted parameters into interpretable patterns [17]. Figure 4a displays extracted patterns for the BNCI2015001 dataset (inter-subject TL). It is clearly visible that TSMNet infers the target label from neurophysiologically plausible sources (rows in Figure 4a). As expected [2], the source with highest contribution has spectral peaks in the alpha and beta bands, and originates in contralateral and central sensorimotor cortex.

5.1.0.2 DSBN on \(\mathcal{S}_D^+\) drives the success of TSMNet

Since TSMNet combines several advances, we present the results of an ablation study in Table 2. It summarizes the grand average inter-session TL test scores relative to TSMNet with SPDDSMBN for n = 138 subjects. We observed three significant effects. The largest effect can be attributed to \(\mathcal{S}_D^+\), as we observed the largest performance decline if the architecture would be modified⁵ to omit the SPD manifold (4.5% with DSBN, 3% w/o DSBN). The performance gain comes at the cost of a 2.6x longer time to fit the parameters. The second largest effect can be attributed to DSBN; without DSBN the performance dropped by 3.9% (with \(\mathcal{S}_D^+\)) and 2.4% (w/o \(\mathcal{S}_D^+\)). The smallest, yet significant effect can be attributed to SPDMBN.

5.2 Mental workload estimation↩︎

Compared to the baseline methods, TSMNet obtained the highest average scores of 54.7% (7.3%) and 52.4% (8.8%) in inter-session and -subject TL (for details see appendix 7.3.1). Interestingly, the inter-session TL score of TSMNet matches the score (54.3%) of the winning method in last year’s competition [15]. To shed light on the sources utilized by TSMNet, we show patterns for a fitted model in Figure 4b. For the low mental workload class, the top contributing source’s activity peaked in the theta band and originated in pre-frontal areas. The second source’s activity originated in occipital cortex with non-focal spectral profile. Our results agree with the findings of previous research, as both areas and the theta band have been implicated in mind wandering and effort withdrawal [60].

6 Discussion↩︎

In this contribution, we proposed a machine learning framework around (domain-specific) momentum batch normalization on \(\mathcal{S}_D^+\) to learn tangent space mapping (TSM) and feature extractors in an end-to-end fashion. In a theoretical analysis, we provided error bounds for the running estimate of the Fréchet mean as well as convergence guarantees under reasonable assumptions. We then applied the framework, to a multi-source multi-target unsupervised domain adaptation problem, namely, inter-session and -subject transfer learning for EEG data and obtained or attained state-of-the art performance with a simple, intrinsically interpretable model, denoted TSMNet, in a total of 6 diverse BCI datasets (138 human subjects, 573 sessions). In the case of mental imagery, we found that TSMNet significantly reduced (inter-subject TL) or even exceeded (inter-session TL) the performance gap to a SoA domain-specific method.

Although our framework could be readily extended to online UDA for unseen target domains, we limited this study to offline evaluations and leave actual BCI studies to future work. A limitation of our framework, and also any other method that involves eigen decomposition, is the computational complexity, which limits its application to high-dimensional SPD features (e.g., fMRI connectivity matrices with fine spatial granularity). Altogether, the presented results demonstrate the utility of our framework and in particular TSMNet as it not only achieves highly competitive results but is also intrinsically interpretable. While we do not foresee any immediate negative societal impacts, we provide direct contributions towards the scalability and acceptability of EEG-based healthcare [1], [5] and consumer [18], [60] technologies. We expect future works to evaluate the impact of the proposed methods in clinical applications of EEG like sleep staging [80], [81], seizure [82] or pathology detection [83], [84].

Checklist↩︎

1. For all authors...

   1. Do the main claims made in the abstract and introduction accurately reflect the paper’s contributions and scope?

   2. Did you describe the limitations of your work?

   3. Did you discuss any potential negative societal impacts of your work?

   4. Have you read the ethics review guidelines and ensured that your paper conforms to them?

2. If you are including theoretical results...

   1. Did you state the full set of assumptions of all theoretical results?

   2. Did you include complete proofs of all theoretical results?

3. If you ran experiments...

   1. Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)?

   2. Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)?

   3. Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)?

   4. Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)?

4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...

   1. If your work uses existing assets, did you cite the creators?

   2. Did you mention the license of the assets?

   3. Did you include any new assets either in the supplemental material or as a URL?

   4. Did you discuss whether and how consent was obtained from people whose data you’re using/curating?

   5. Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content?

5. If you used crowdsourcing or conducted research with human subjects...

   1. Did you include the full text of instructions given to participants and screenshots, if applicable?

   2. Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable?

   3. Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation?

7 Supplementary Material↩︎

7.1 Proof of proposition 1↩︎

In the proof, we will use Theorem 2 of [50] which relates the distance between interpolated points along the geodesic \(\mathbf{R} \#_\gamma \mathbf{S}\), connecting points \(\mathbf{R}\) and \(\mathbf{S}\), and point \(\mathbf{T}\) to the distances between \(\mathbf{R}\), \(\mathbf{S}\) and \(\mathbf{T}\). Formally, for all \(\mathbf{R},\mathbf{S}, \mathbf{T} \in \mathcal{S}_D^+\) we have

\[\delta_{}^2(\mathbf{R} \#_\gamma \mathbf{S}, \mathbf{T}) \le (1-\gamma) \delta_{}^2(\mathbf{R}, \mathbf{T}) + \gamma \delta_{}^2(\mathbf{S}, \mathbf{T}) - \gamma (1-\gamma) \delta_{}^2(\mathbf{R}, \mathbf{S}) \label{eq:ho2013:thm3}\tag{19}\]

As a last ingredient for the proof, we use Proposition 1 of [50], which states that for a random variable \(\mathbf{U}\), following distribution \(P_{\mathcal{U}}\) defined on \(\mathcal{S}_D^+\) with Fréchet mean \(\mathbf{G}_{\mathcal{U}}\), we have for any point \(\mathbf{V} \in \mathcal{S}_D^+\):

\[\mathbb{E} \{ \delta_{}^2(\mathbf{U}, \mathbf{V}) \} \ge \underbrace{\int \delta_{}^2(\mathbf{u}, \mathbf{G}_{\mathcal{U}}) dP_{\mathcal{U}}(\mathbf{u})}_{=: \mathrm{Var}(\mathbf{U}) } + \delta_{}^2(\mathbf{V}, \mathbf{G}_{\mathcal{U}}) \label{eq:ho2013:prop1}\tag{20}\]

which means that the expected distance between \(\mathbf{V}\) and \(\mathbf{U}\) is bounded from below by the variance of \(\mathbf{U}\) and the distance between \(\mathbf{V}\) and its Fréchet mean \(\mathbf{G}_{\mathcal{U}}\).

Let us quickly repeat the definitions, assumptions and proposition 1. We define a dataset containing the latent representations generated with feature set \(\theta_k\) as \(\mathcal{Z}_{\theta_k} = \{ f_{\theta_k} (\mathbf{x}) | \mathbf{x} \in \mathcal{T} \}\), and a minibatch of \(M\) iid samples drawn from \(\mathcal{Z}_{\theta_k}\) at training step \(k\) as \(\mathcal{B}_k\). We denote the Fréchet mean of \(\mathcal{Z}_{\theta_k}\) as \(\mathbf{G}_{\theta_k}\), the estimated Fréchet mean, defined in (9 ), as \(\mathbf{G}_k\), and the estimated batch mean as \(\mathbf{B}_k\). Since the batches are drawn randomly, we consider \(\mathbf{B}_k\) and \(\mathbf{G}_k\) as random variables.

We assume that the variance \(\mathrm{Var}_{\theta_{k}}(\mathbf{B}_{k-1}) = \mathbb{E}_{\mathbf{B}_{k-1}} \{ \delta^2(\mathbf{B}_{k-1}, \mathbf{G}_{\theta_{k}}) \}\) of the previous batch mean \(\mathbf{B}_{k-1}\) with respect to the current Fréchet mean \(\mathbf{G}_{\theta_k}\) is bounded by the current variance \(\mathrm{Var}_{\theta_{k}}(\mathbf{B}_{k})\) and the norm of the difference in the parameters

\[\mathrm{Var}_{\theta_k}(\mathbf{B}_{k-1}) \le (1 + || \theta_k - \theta_{k-1} ||) \mathrm{Var}_{\theta_k}(\mathbf{B}_{k}) \label{eq:proof:assumption1}\tag{21}\]

Proposition 1 states then that the variance of the running mean \(\mathrm{Var}_{\theta_k}(\mathbf{G}_k)\) is bounded by \[\mathrm{Var}_{\theta_k}(\mathbf{G}_k) \le \mathrm{Var}_{\theta_k}(\mathbf{B}_k)\] over training steps \(k\), if \[|| \theta_k - \theta_{k-1} || \le \frac{1- \gamma^2}{(1-\gamma)^2} - 1\] holds true.

7.2 Supplementary Methods↩︎

7.2.1 Datasets↩︎

A summary of the datasets’ key attributes is listed in Table 3. The datasets contain a diverse sample of 154 human subjects, whose data was acquired in Europe (38 subjects; BNCI2014001,BNCI2015001,Lehner2021,Hehenberger2021,Hinss2021), Asia (54 subjects; Lee2019) and North America (62 subjects; Stieger2021).

7.2.2 Preprocessing↩︎

Depending on the dataset, either all or a subset of EEG channels was selected, and then resampled along the temporal dimension to a sampling rate of either 250 or 256 Hz. (see Table 3). Thereafter, an infinite impulse response (IIR) bandpass filter was used to extract EEG activity in the 4 to 36 Hz range (4th order Butterworth filter, 4 and 36 Hz cut-off frequencies, zero-phase). Some baseline methods required spectrally resolved input data. For these, we applied a bank of 8 filters with similar parameters except for the cut-off frequencies ([4, 8], [8, 12], ..., [32, 36] Hz). Finally, short epochs (=segments) were extracted (see Table 3) relative to the task cues (=labels). The labeled data were then extended with a domain index (= unique integer associated to one session of one subject).

7.2.3 Baseline methods↩︎

We considered several established and SoA baseline methods which were previously applied to inter-session/-subject TL. They can be broadly categorized as component based, Riemannian geometry aware or deep learning which we denote component, geometric and end-to-end, respectively. For the component category, we considered the popular filter-bank common spatial patterns (FBCSP+SVM) [71] and a variant [68], designed for MSMTUDA, that applies domain-specific standardization (DSS) to features before classification, denoted FBCSP+DSS+LDA. The geometric category was represented by TSM+SVM [12], a spectrally resolved variant [17] denoted FB+TSM+LR (which was also used as domain-specific baseline method). We additionally considered two MSUDA methods, denoted URPA+MDM and SPDOT+TSM+SVM here, that align SPD observations (=spatial covariance matrices) of different domains. The former uses Riemannian Procrustes Analysis (RPA) [30] to align domains, and the latter optimal transport (OT) on \(\mathcal{S}^+_D\) [72]. The end-to-end category was represented by EEGNet [73] and ShConvNet [56] two convolutional neural network architectures specifically designed for EEG data. We additionally considered variants [25] that use domain-adversarial neural networks (DANN) [85] to learn domain-invariant latent representations.

7.2.4 Domain-specific reference method↩︎

Due to the success of TSM models [10], [14], we considered a spectrally resolved model [13], [17] which consisted of a filter-bank to separate activity of canonical frequency bands. For each frequency band, PCA was used to reduce the spatial dimensionality and TSM to project the SPD features to the Euclidean vector space. Finally, all features were pooled and submitted to a penalized logistic regression classifier. For further details, see [17].

7.2.5 TSMNet↩︎

7.2.5.1 Architecture

The architecture of TSMNet is outlined in Figure 5 and detailed in Table 4. The feature extractor \(f_\theta\) comprises two convolutional layers, followed by covariance pooling [57], BiMap [37] and ReEig [37] layers. The first convolutional layer performs convolution along the temporal dimension; implementing a finite impulse response (FIR) filter bank (4 filters) with learnable parameters. The second convolutional layer applies spatio-spectral filters (40 filters) along the spatial and convolutional channel dimensions. Covariance pooling is then applied along the temporal dimension. A subsequent BiMap layer projects covariance matrices to a subspace via a bilinear mapping (i.e, \(\mathrm{BiMap}(\mathbf{Z}) = \mathbf{W}^T_\theta \mathbf{Z} \mathbf{W}_\theta\)) where the parameter matrix \(\mathbf{W}_\theta\) is constrained to have orthogonal rows (i.e., \(\mathbf{W}_\theta \in \{ \mathbf{U} \in \mathbb{R}^{I \times O} : \mathbf{U}^T \mathbf{U} = \mathbf{I}_O, I \ge O \}\)). Next, a ReEig layer rectifies all eigenvalues, lower than a threshold \(\epsilon = 10^{-4}\) (i.e, \(\mathrm{ReEig}(\mathbf{Z}) = \mathbf{U} \mathrm{max}(\mathbf{\Sigma},\epsilon \mathbf{I}) \mathbf{U}^T\) with \([\mathbf{\Sigma},\mathbf{U}] = \mathrm{eig}(\mathbf{Z})\)).
Domain-specific tangent space mapping \(m_\phi\) is implemented via combining SPDDSMBN and LogEig layers. In order for \(m_{\phi}\) to align the domain data and compute TSM, we use SPDDSMBN (16 ) with shared parameters (i.e., \(\mathbf{G}_{\phi_i} = \mathbf{G}_{\phi} = \mathbf{I}, \nu_{\phi_i} = \nu_{\phi}\)) in (14 ). The classifier \(g_{\psi}\) was parametrized as a linear layer with softmax activations.

7.2.5.2 Parameter estimation

We used the cross entropy loss as training objective, and the standard backpropagation framework with extensions for structured matrices [49] and manifold-constrained gradients [40] to propagate gradients through the layers of TSMNet. Gradients were estimated with mini-batches of fixed size (50 observations; 10 per domain; 5 distinct domains) and converted into parameter updates with the Riemannian ADAM optimizer [58] (\(10^{-3}\) learning rate, \(10^{-4}\) weight decay applied to unconstrained parameters; \(\beta_1 = 0.9, \beta_2 = 0.999\)).

For every MSMTUDA problem, comprising source and target domain sets, we split the source domains’ data into training and validation sets (80%/20% splits; randomized; stratified across domains and labels) and repeatedly iterated through the training set for 50 epochs via exhaustive minibatch sampling. As required by SPDMBN, we implemented a decaying schedule (over epochs) for the training momentum parameter \(\gamma_{train}(k)\), defined in (11 ), with \(\gamma_{train}(0) = 1\) and \(\gamma_{min} = 0.2\) attained at epoch \(K = 40\). During training, we monitored the loss on the validation data (at the end of every epoch). Post training, the model with minimal loss on the validation data was selected. For each target domain, the associated data was then passed through this model to estimate the labels. During the forward pass, the domain’s normalization statistics within the SPDMBN layers were computed by solving (4 ) with the Karcher flow algorithm [42].

7.3 Supplementary results↩︎

7.3.1 EEG data↩︎

The test set score for each considered EEG dataset is summarized in in Table 5. Significant differences between the proposed method (TSMNet with SPDDSMBN) and baseline methods are highlighted.

References↩︎