SPDIM: Source-Free Unsupervised Conditional and Label Shift Adaptation in EEG

## Abstract

Abstract The non-stationary nature of electroencephalography (EEG) introduces distribution shifts across domains (e.g., days and subjects), posing a significant challenge to EEG-based neurotechnology generalization.
Without labeled calibration data for target domains, the problem is a source-free unsupervised domain adaptation (SFUDA) problem.
For scenarios with constant label distribution, Riemannian geometry-aware statistical alignment frameworks on the symmetric positive definite (SPD) manifold are considered state-of-the-art.
However, many practical scenarios, including EEG-based sleep staging, exhibit label shifts.
Here, we propose a geometric deep learning framework for SFUDA problems under specific distribution shifts, including label shifts.
We introduce a novel, realistic generative model and show that prior Riemannian statistical alignment methods on the SPD manifold can compensate for specific marginal and conditional distribution shifts but hurt generalization under label shifts.
As a remedy, we propose a parameter-efficient manifold optimization strategy termed SPDIM.
SPDIM uses the information maximization principle to learn a single SPD-manifold-constrained parameter per target domain.
In simulations, we demonstrate that SPDIM can compensate for the shifts under our generative model.
Moreover, using public EEG-based brain-computer interface and sleep staging datasets, we show that SPDIM outperforms prior approaches.

## 1 Introduction

Electroencephalography (EEG) measures multi-channel electric brain activity from the human scalp (Niedermeyer & da Silva, ) and can reveal cognitive processes (Pfurtscheller & Da Silva, ), emotion states (Suhaimi et al., ), and health status (Alotaiby et al., ).
Neurotechnology and brain-computer interfaces (BCI) aim to extract patterns from the EEG activity that can be utilized for various applications, including rehabilitation and communication (Wolpaw et al., ).
Despite their capabilities, they currently suffer from a low signal-to-noise ratio, low specificity, and non-stationarities manifesting as distribution shifts across days and subjects (Fairclough & Lotte, ).

For EEG-based neurotechnology, distribution shifts have been traditionally mitigated by collecting labeled calibration data and training domain-specific models (Lotte et al., ), limiting neurotechnology utility and scalability (Wei et al., ).
As an alternative, domain adaptation (DA) learns a model from one or multiple source domains that performs well on different (but related) target domain(s), offering principled statistical learning approaches with theoretical guarantees (Ben-David et al., ; Hoffman et al., ).
Within the BCI field, DA primarily addresses cross-session and cross-subject transfer learning (TL) problems (Wu et al., ), aiming to achieve robust generalization across domains (e.g., sessions and subjects) without supervised calibration data.

BCIs that generalize across domains without requiring labeled calibration data are one of the grand challenges in EEG-based BCI research (Fairclough & Lotte, ; Wolpaw et al., ).
Since target domain data is typically not available during training, the problem corresponds to a source-free unsupervised domain adaptation (SFUDA) problem (Liu et al., ; Yang et al., ).
For this problem class, Riemannian geometry-aware statistical alignment frameworks (Barachant et al., ) operating with symmetric, positive definite (SPD) matrix-valued features are considered state-of-the-art (SoA) in cross-domain (Roy et al., ; Mellot et al., ) generalization.
They offer several advantageous properties, including invariance to linear mixing, that are suitable for EEG data (Congedo et al., ).
Additionally, prior works proposed consistent (Sabbagh et al., ) and inherently interpretable (Kobler et al., ) estimators for generative models with a log-linear relationship between the power of latent sources and the labels.
To facilitate generalization across domains, statistical alignment frameworks aim to align first (Zanini et al., ; Yair et al., ) and second (Rodrigues et al., ; Kobler et al., ) order moments, denoted Fréchet mean and variance on Riemannian manifolds.
Once the moments are aligned, a model trained on source domains typically generalizes to related target domains (Zanini et al., ; Wei et al., ; Kobler et al., ; Mellot et al., ).

Figure: Figure 1: Visualization of SPDIM framework. a, EEG data is obtained by a domain-specific linear mixing from source signals. b, visualization of the latent representations learned by the feature extractor. c, visualization of the latent alignment without (top) and with (bottom) SPDIM under label shifts. d, density plots show that SPDIM better aligns the distribution.
Refer to caption: extracted/5956360/figures/fig_overview.png

Although infrequently studied, many applications, including EEG-based sleep staging, exhibit label shifts (Thölke et al., ).
Under label shifts, aligning the moments of the marginal feature distributions can increase the generalization error (Bakas et al., ).
To address various sources of distribution shifts in EEG, an SFUDA approach that can deal with additional label shifts is required.
Machine learning literature offers several frameworks for SFUDA problems with label shifts (Li et al., ; Liang et al., ), but few have been applied to EEG data.
For example, Li et al. () employed the information maximization (Shi & Sha, ) objective for cross-domain generalization.
Within Riemannian geometry methods, Mellot et al. () studied an EEG-based age regression problem and proposed a framework to facilitate generalization across populations with different prior distributions.

Here, we propose a geometric deep learning framework to tackle SFUDA classification problems under distribution shifts, including label shifts.
We introduce a realistic generative model with a log-linear relationship between the covariance of latent sources and the labels (Figure 1a).
We provide theoretical analyses showing that prior Riemannian statistical alignment methods on the SPD manifold can compensate for the conditional distribution shifts introduced in our model but hurt generalization under additional label shifts (Figure 1b, c (top)).
As a remedy, we propose a parameter-efficient manifold optimization strategy termed SPDIM.
SPDIM employs the information maximization principle to learn a domain-specific SPD-manifold-constrained bias parameter to compensate over-corrections introduced via aligning the Fréchet mean (Figure 1b, c (bottom)).

## 2 Preliminaries

### 2.1 SFUDA Scenario

Let $\mathrm{y}$ and $\mathrm{x}$ denote random variables representing the true labels and associated features, and let $p_{j}(\mathrm{y})$ and $p_{j}(\mathrm{x}|\mathrm{y})$ be the prior and conditional probability distributions for domain $j$.
In the transfer learning scenario considered here, we assume that the class priors $p_{j}(\mathrm{y})$ can be different from each other, as well as specific shifts in the conditional distributions $p_{j}(\mathrm{x}|\mathrm{y})$.
We consider a labeled source dataset $\mathcal{D}_{s}=\{(x_{i},y_{i},j_{i})|x_{i}\in\mathcal{X}_{s},y_{i}\in\mathcal
{Y}_{s},j_{i}\in\mathcal{J}_{s}\}_{i=1}^{L_{s}}$ and an unlabeled target dataset $\mathcal{D}_{t}=\{(x_{i},j_{i})|x_{i}\in\mathcal{X}_{t},y_{i}\in\mathcal{Y}_{t
},j_{i}\in\mathcal{J}_{t}\}_{i=1}^{L_{t}}$, where $y_{i}$ and $j_{i}$ (or $j(i)$ if used as suffix) indicate the associated label and domain for each example $i$.
We additionally assume that both datasets share the feature space (i.e., $\mathcal{X}_{s}{=}\mathcal{X}_{t}{=}\mathcal{X}$), contain the same classes (i.e., $\mathcal{Y}_{s}{=}\mathcal{Y}_{t}{=}\mathcal{Y}$), and comprise examples from different domains (i.e., $\mathcal{J}_{s}\cap\mathcal{J}_{t}{=}\varnothing$).
The goal is to transfer the knowledge learned from $\mathcal{D}_{s}$ to $\mathcal{D}_{t}$ via first learning a source decoder $h_{s}$ within hypothesis class $\mathcal{H}$ and then use $h_{s}$ and the unlabeled target dataset $\mathcal{D}_{t}$ to learn $h_{t}\in\mathcal{H}$.

### 2.2 Riemannian Geometry on 𝒮 D + subscript superscript 𝒮 𝐷 \mathcal{S}^{+}_{D} caligraphic_S start_POSTSUPERSCRIPT + end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_D end_POSTSUBSCRIPT

The SPD manifold $\mathcal{S}^{+}_{D}=\{C\in\mathbb{R}^{D\times D}:C^{T}{=}C,C\succ 0\}$ together with an inner product on its tangent space $\mathcal{T}_{C}\mathcal{S}^{+}_{D}$ at each point $C\in\mathcal{S}^{+}_{D}$ forms a Riemannian manifold.
Tangent spaces have Euclidean structure with easy-to-compute distances, which locally approximate Riemannian distances on $\mathcal{S}^{+}_{D}$ (Absil et al., ).
In this work, we consider the affine invariant Riemannian metric (AIRM) $g_{C}^{\text{AIRM}}(S_{1},S_{2}){=}\mathrm{Tr}(C^{-1}S_{1}C^{-1}S_{2})$ as the inner product, which gives rise to the following distance (Bhatia, ):

$$ $\delta(C_{1},C_{2})=\|\log(C_{1}^{-\frac{1}{2}}C_{2}C_{1}^{-\frac{1}{2}})\|_{F}$ (1) $$

where $C_{1}$ and $C_{2}$ are two SPD matrices, $\mathrm{Tr}(\cdot)$ denotes the trace, $\log(\cdot)$ the matrix logarithm, and $\|\cdot\|_{F}$ the Frobenius norm.
For a set of points $\mathcal{C}=\{C_{i}\in\mathcal{S}_{+}^{D}\}_{i\leq M}$, the Fréchet mean is defined as the minimizer of the average squared distances:

$$ $\bar{C}=\arg\min_{G\in\mathcal{S}_{D}^{+}}\nu_{G}(\mathcal{C})=\arg\min_{G\in \mathcal{S}_{D}^{+}}\frac{1}{M}\sum_{j=1}^{M}\delta^{2}(G,C_{i})$ (2) $$

For $M=2$, there exists a closed form solution:

$$ $\bar{C}=\arg\min_{G\in\mathcal{S}_{D}^{+}}\nu_{G}(\{C_{1},C_{2}\})=\left.C_{1} \#_{t}C_{2}\right\rvert_{t=\frac{1}{2}}=\left.C_{1}^{\frac{1}{2}}\left(C_{1}^{ -\frac{1}{2}}C_{2}C_{1}^{-\frac{1}{2}}\right)^{t}C_{1}^{\frac{1}{2}}\right \rvert_{t=\frac{1}{2}}$ (3) $$

where parameter $t\in[0,1]$ smoothly interpolates along the geodesic $C_{1}\#_{t}C_{2}$ (i.e., the shortest path) connecting both points.

The logarithmic map $\text{Log}_{\bar{C}}:\mathcal{S}^{+}_{D}\rightarrow\mathcal{T}_{\bar{C}}
\mathcal{S}^{+}_{D}$ and exponential map $\text{Exp}_{\bar{C}}:\mathcal{T}_{\bar{C}}\mathcal{S}^{+}_{D}\ \rightarrow
\mathcal{S}^{+}_{D}$ project points between the manifold and the tangent space at point $\bar{C}$:

$$ $\displaystyle\text{Log}_{\bar{C}}(C_{i})=\,$ $\displaystyle{\bar{C}}^{\frac{1}{2}}\log({\bar{C}}^{-\frac{1}{2}}C_{i}{\bar{C} }^{-\frac{1}{2}}){\bar{C}}^{\frac{1}{2}}$ (4) $\displaystyle\text{Exp}_{\bar{C}}(S_{i})=\,$ $\displaystyle{\bar{C}}^{\frac{1}{2}}\exp({\bar{C}}^{-\frac{1}{2}}S_{i}{\bar{C} }^{-\frac{1}{2}}){\bar{C}}^{\frac{1}{2}}$ (5) $$

To transport points $S_{i}\in\mathcal{T}_{\bar{C}}\mathcal{S}^{+}_{D}$ from the tangent space at $\bar{C}$ to the tangent space at $\bar{C}_{\phi}$, parallel transport on $\mathcal{S}_{+}^{D}$ can be used as:

$$ $\Gamma_{\bar{C}\rightarrow\bar{C}_{\phi}}(S_{i})=P^{\top}S_{i}P,\quad P=\left( \bar{C}^{-1}\bar{C}_{\phi}\right)^{\frac{1}{2}}$ (6) $$

While parallel transport is generally defined for tangent space vectors (Absil et al., ), for $(\mathcal{S}^{+}_{D},g_{\cdot}^{\text{AIRM}})$ it can be directly applied without explicitly computing tangent space projections (i.e., $\text{Exp}_{\bar{C}_{\phi}}\circ\Gamma_{\bar{C}\rightarrow\bar{C}_{\phi}}\circ
\text{Log}_{\bar{C}}=\Gamma_{\bar{C}\rightarrow\bar{C}_{\phi}}$) (Brooks et al., ; Yair et al., ).

If $\bar{C}_{\phi}$ lies along the geodesic connecting $\bar{C}$ with the identity matrix $I_{D}$, there exists a step-size $t\in\mathbb{R}$ so that $\bar{C}_{\phi}=\bar{C}\#_{t}I_{D}$, and (6) simplifies to Mellot et al. ():

$$ $\Gamma_{\bar{C}\rightarrow I_{D}}(C_{i};t)=\Gamma_{\bar{C}\rightarrow\bar{C}\# _{t}I_{D}}(C_{i})=\bar{C}^{-\frac{t}{2}}C_{i}\bar{C}^{-\frac{t}{2}}$ (7) $$

## 3 Methods

### 3.1 Generative Model

In the case of EEG, the features $x_{i}\in\mathcal{X}\subseteq\mathbb{R}^{P\times T}$ comprise epochs of multivariate time-series data with $P$ spatial channels and $T$ consecutive temporal samples.
Propagation of brain activity to the EEG electrodes, located at the scalp, is typically modeled as a linear mixture of sources (Nunez & Srinivasan, ):

$$ $x_{i}=A_{j(i)}z_{i}$ (8) $$

where $A_{j(i)}\in\{A\in\mathbb{R}^{P\times P}:\text{rank}(A){=}P\}$ is a domain-specific forward model and $z_{i}\in\mathbb{R}^{P\times T}$ the activity of latent sources.
Utilizing the uniqueness of the polar decomposition for invertible matrices, we constrain the model to

$$ $A_{j}=Q\mathrm{exp}(P_{j})$ (9) $$

where $Q\in\{Q\in\mathbb{R}^{P\times P}:Q^{T}Q{=}I_{P}\}$ is a orthogonal matrix modeling rotations and $\mathrm{exp}(P_{j})\in\mathcal{S}^{+}_{P}$ a SPD matrix modeling domain-specific scalings.

Like (Sabbagh et al., ; Kobler et al., ; Mellot et al., ), we consider zero-mean (i.e., $\mathbb{E}\{z_{i}\}=0_{P}~{}\forall\,i$) signals, and a log-linear relationship between the spatial covariance $E_{i}=\mathrm{Cov}(z_{i})\in\mathcal{S}^{+}_{P}$ of the latent sources $z_{i}$ and the target $y_{i}$.
Specifically, we model the source covariance matrices as:

$$ $\log(E_{i})=S_{i},~{}~{}\mathrm{upper}(S_{i})=s_{i}$ (10) $$

where the invertible mapping $\mathrm{upper}(S)\in\mathbb{R}^{P(P+1)/2}$ extracts the upper triangular elements of a symmetric matrix while preserving its norm (i.e., $\|\mathrm{upper}(S)\|_{2}{=}\|S\|_{F}$), and latent log-space features $s_{i}\in{R^{P(P+1)/2}}$ are generated as:

$$ $s_{i}=B\tilde{y}_{i}+\varepsilon_{i}$ (11) $$

where $B\in\{B\in\mathbb{R}^{P(P+1)/2\times|\mathcal{Y}|}:\mathrm{rank}(B){=}|
\mathcal{Y}|\}$, $\tilde{y}_{i}=[\delta(y_{i},k)-p_{j}(\mathrm{y}{=}y_{k})]_{k=1,..,|\mathcal{Y}|}$ represents one-hot-coded, zero-mean labels, and $\varepsilon_{i}\in\mathbb{R}^{P(P+1)/2}$ zero-mean additive noise.

###### Proposition 1

Given the specified generative model and a set of examples $\mathcal{E}_{j}=\{(E_{i},y_{i},j)|E_{i}\in\mathcal{S}_{P}^{+},y_{i}\in\mathcal
{Y}\}_{i\leq M_{j}}$ of domain $j\in\mathcal{J}$, we have that the Fréchet mean of $\mathcal{E}_{j}$, defined in (2), converges to the identity matrix $I_{P}$ with $M_{j}\rightarrow\infty$ for all domains $j\in\mathcal{J}$.

Our proof, provided in appendix A.1, relies on the uniqueness of the Fréchet mean for $(\mathcal{S}^{+}_{P},g_{\cdot}^{\text{AIRM}})$ and that $\tilde{y}_{i}$ and $\varepsilon_{i}$ in (11) are zero-mean.

For the generated multi-variate time-series features $x_{i}\in\mathcal{X}$ the empirical covariance matrix $C_{i}$ is:

$$ $C_{i}:=\mathrm{Cov}(x_{i})=\frac{1}{T}x_{i}x_{i}^{T}\in\mathcal{S}_{P}^{+}$ (12) $$

Due to the linear relationship between observed features and latent sources (8), we obtain a direct relationship to the latent source covariance matrices $E_{i}$:

$$ $C_{i}=A_{j(i)}E_{i}A_{j(i)}^{T}$ (13) $$

Since $E_{i}$ encodes the target $y_{i}\in\mathcal{Y}$ and $A_{j}$ is invertible, $C_{i}$ are sufficient descriptors to decode $y_{i}$.

###### Remark 1

Note that for each domain $j$ the covariance matrices of the generated data $\{C_{i}:j_{i}=j\}$ are not necessarily jointly diagonalizable.
Depending on the structure of $B$ and $\varepsilon_{i}$, the proposed generative model reduces to a jointly diagonalizable model if all the off-diagonal elements in $\log(E_{i})$ are zero $\forall i$.

###### Proposition 1

###### Remark 1

### 3.2 Decoding Framework

Given source domain data $\mathcal{D}_{s}$ and a hypothesis class $\mathcal{H}$, we aim to learn a decoder function $h_{s}\in\mathcal{H}$, and - once the unlabeled target data $\mathcal{D}_{t}$ is revealed - use $\mathcal{D}_{t}$ and $h_{s}$ to learn $h_{t}\in\mathcal{H}$.
Following (Kobler et al., ), we constrain the hypothesis class $\mathcal{H}$ to functions $h:\mathcal{X}\times\mathcal{J}\rightarrow\mathcal{Y}$ that can be decomposed into a composition of a shared feature extractor $f_{\theta}:\mathcal{X}\rightarrow\mathcal{S}_{D}^{+}$ (Figure 1b), latent alignment $m_{\phi}:\mathcal{S}_{D}^{+}\times\mathcal{J}\rightarrow\mathbb{R}^{D(D+1)/2}$ (Figure 1c), and a shared linear classifier $g_{\psi}:\mathbb{R}^{D(D+1)/2}\rightarrow\mathcal{Y}$ (Figure 1d) with parameters $\Theta=\{\theta,\phi,\psi\}$.
Within this section, we focus on theoretical considerations for $m_{\phi}$ under our generative model.

#### Tangent space mapping (TSM) to recover log ⁡ ( E i ) subscript 𝐸 𝑖 \log(E_{i}) roman_log ( italic_E start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT )

Considering a set of labeled data $\mathcal{D}=\{(x_{i},y_{i}):x_{i}\in\mathcal{X},y_{i}\in\mathcal{Y},j_{i}{=}j~
{}\forall i\}$ obtained from a single domain $j$, TSM (Barachant et al., ) provides an established (Lotte et al., ; Jayaram & Barachant, ) decoding approach to infer $y_{i}$.
TSM requires SPD-matrix valued representations.
For the considered generative model, covariance features $C_{i}$, as defined in (5), are a natrual choice (i.e., $f_{\theta}\triangleq\mathrm{Cov}$).
In a nutshell, TSM first estimates the Fréchet mean $\bar{C}$ of $\mathcal{C}=\{\mathrm{Cov}(x_{i}):(x_{i},y_{i})\in\mathcal{D}\}$ , projects each ${C}_{i}$ to the tangent space at $\bar{C}$, and finally transports the data to vary around $I_{P}$. Formally,

$$ $\tilde{m}_{\phi}(C_{i}):=\mathrm{upper}\circ\Gamma_{\bar{C}\rightarrow I}\circ \mathrm{Log}_{\bar{C}}(C_{i})=\mathrm{upper}\left(\log\left(\bar{C}^{-\frac{1} {2}}\,C_{i}\,\bar{C}^{-\frac{1}{2}}\right)\right),~{}\phi=\{\bar{C}\}$ (14) $$

where the resulting representations $\tilde{m}_{\phi}(C_{i})\in\mathbb{R}^{P(P+1)/2}$ have a linear relationship to $\log(E_{i})$ (Sabbagh et al., ; Kobler et al., ) and through (10) and (11) also to the labels $y_{i}$.

#### RCT+TSM compensates marginal and conditional shifts

In the context of functional neuroimaging data, the recentering (RCT) transform (Zanini et al., ; Yair et al., ) and its extensions (Rodrigues et al., ; He & Wu, ; Kobler et al., ; Mellot et al., ) address source-free UDA problems.
Combined with TSM, RCT+TSM essentially applies (14) independently to each domain $j$.
The outputs $\tilde{m}_{\phi(j(i))}(C_{i})$, where $\phi(j)=\{\bar{C}_{j}\}$, are treated as domain-invariant and passed on to the shared classifier $g_{\psi}:\mathbb{R}^{D(D+1)/2}\rightarrow\mathcal{Y}$.

Although RCT+TSM is an established method to address SFUDA problems for neuroimaging data (Lotte et al., ; Wei et al., ; Roy et al., ), there is a lack of understanding of what kind of distribution shifts RCT+TSM can compensate.

###### Proposition 2

For the generative model, specified in section 3.1, RCT+TSM compensates conditional distribution shifts introduced by the invertible linear map $A_{j}$, defined in (8), and recovers domain-invariant representations if there are no label shifts (i.e., $p_{j}(\mathrm{y}{=}y){=}p(\mathrm{y}{=}y)~{}\forall y\in\mathcal{Y},~{}\forall
j
\in\mathcal{J}_{s}\cup\mathcal{J}_{t}$).

A detailed proof is provided in appendix A.2.
Starting with $\tilde{m}_{\phi(j(i))}(C_{i})$, as defined in (14) and utilizing the unique polar decomposition (9) of $A_{j}$ into rotational $Q$ and scaling $\mathrm{exp}(P_{j})$ transformations along with proposition 1 we obtain:

$$ $\log\left(\bar{C}_{j(i)}^{-\frac{1}{2}}C_{i}\bar{C}_{j(i)}^{-\frac{1}{2}} \right)=Q\log(E_{i})Q^{T}=Q\mathrm{upper}^{-1}\left(B\left(1_{y_{i}}-\bar{y}_{ j(i)}\right)+\varepsilon_{i}\right)Q^{T}$ (15) $$

where the vector $\bar{y}_{j}=[p_{j}(\mathrm{y}{=}y)]_{y=1,...,|\mathcal{Y}|}$ contains the class priors for domain $j$.
If there are no label shifts the class priors are constant $p_{j}(\mathrm{y}{=}y)=p(\mathrm{y}{=}y)~{}\forall y\in\mathcal{Y},~{}\forall j
\in\mathcal{J}_{s}\cup\mathcal{J}_{t}\Longrightarrow\bar{y}_{j}=\bar{y}$, and (15) only contains domain-invariant terms on the right hand side.
Thus, RCT+TSM compensates the conditional shifts introduced by $\mathrm{exp}(P_{j})$. $\square$

###### Remark 2

If all sources in (8) can be partitioned into relevant and irrelevant sources $z_{i}{=}[z_{i}^{\parallel},z_{i}^{\bot}],z_{i}^{\parallel}{=}f(y_{i}),z_{i}^{
\bot}{\neq}f(y_{i})$ that are independent from each other, the source covariance matrices $E_{i}$ have a block-diagonal structure.
Then, RCT+TSM compensates marginal shifts in $z_{i}^{\bot}$ and conditional shifts in $z_{i}^{\parallel}$ introduced by $\mathrm{exp}(P_{j})$.

###### Proposition 2

###### Remark 2

#### Alignment under label shifts

We aim to extend RCT+TSM to extract label shift invariant representations.
Specifically, we aim to apply additional transformations on $(\mathcal{S}^{+}_{D},g_{\cdot}^{\text{AIRM}})$ that attenuate the effect of class priors $\bar{y}_{j}$ in (15).
We first rewrite (15) as:

$$ $\displaystyle\bar{C}_{j(i)}^{-\frac{1}{2}}C_{i}\bar{C}_{j(i)}^{-\frac{1}{2}}$ $\displaystyle=\exp\left(Q\log(E_{i})Q^{T}\right)=\exp\left(Q\mathrm{upper}^{-1 }\left(B\left(1_{y_{i}}-\bar{y}_{j(i)}\right)+\varepsilon_{i}\right)Q^{T}\right)$ (16) $\displaystyle=\exp\left(Q\underbrace{\mathrm{upper}^{-1}\left(B1_{y_{i}}+ \varepsilon_{i}\right)}_{\log(\tilde{E}_{i})}Q^{T}-Q\underbrace{\mathrm{upper} ^{-1}\left(B\bar{y}_{j(i)}\right)}_{\bar{P}_{j(i)}}Q^{T}\right)$ (17) $$

where we split $\log(E_{i})$ into domain-invariant $\log(\tilde{E}_{i})$ and label shift $\bar{P}_{j(i)}$ terms.
To separate both terms into products of matrices, we utilize $\exp(\alpha(A+B))=\exp(\alpha B/2)\exp(\alpha A)\exp(\alpha B/2)+\mathcal{O}(
\alpha^{3})$ (Higham, , Theorem 10.5), resulting in:

$$ $\bar{C}_{j(i)}^{-\frac{1}{2}}C_{i}\bar{C}_{j(i)}^{-\frac{1}{2}}=Q\exp\left( \bar{P}_{j(i)}\right)^{-\frac{1}{2}}Q^{T}Q\tilde{E}_{i}Q^{T}Q\exp\left(\bar{P} _{j(i)}\right)^{-\frac{1}{2}}Q^{T}+\mathcal{O}(||\log(E_{i})||_{F}^{3})$ (18) $$

with the approximation error $\mathcal{O}(||\log(E_{i})||_{F}^{3})$ decaying cubically for $||\log(E_{i})||_{F}<1$.
In this form, it is straightforward to see that an additional bilinear transformation on the left hand side in (18) with an SPD matrix can approximately compensate the effect of $\exp(\bar{P}_{j(i)})$.
We denote this parameter as domain-specific bias parameter $\Phi_{j(i)}\in\mathcal{S}_{D}^{+}$, and generalize RCT+TSM to:

$$ $m_{\phi(j(i))}(C_{i})=\mathrm{upper}\circ\log\left(\Phi_{j(i)}^{\frac{1}{2}} \bar{C}_{j(i)}^{-\frac{1}{2}}\,C_{i}\,\bar{C}_{j(i)}^{-\frac{1}{2}}\Phi_{j(i)} ^{\frac{1}{2}}\right),~{}\phi(j(i))=\{\bar{C}_{j(i)},\Phi_{j(i)}\}$ (19) $$

Note that if $\exp(\bar{P}_{j})$ and $\exp({P}_{j})$ share the same eigenvectors, they commute and lie on the same geodesic connecting $\exp({P}_{j})$ with $I_{D}$.
Consequently, the combined effect of $\exp(\bar{P}_{j})$ and $\exp({P}_{j})$ is constrained to the geodesic connecting $\bar{C}_{j}$ with $I_{D}$.
Then, the solution space can be constrained to $\Phi_{j}\in\{\Phi:\Phi\in\mathcal{S}_{D}^{+},\Phi=\bar{C}_{j}\#_{\varphi_{j}}I
_{D},\varphi_{j}\in\mathbb{R}\}$, and (7) used to simplify (19) to:

$$ $m_{\phi(j(i))}(C_{i})=\mathrm{upper}\circ\log\circ\Gamma_{\bar{C}_{j(i)} \rightarrow\bar{C}_{j(i)}\#_{\varphi_{j(i)}}I_{D}}(C_{i}),~{}\phi(j(i))=\{\bar {C}_{j(i)},\varphi_{j(i)}\}$ (20) $$

where the geodesic step-size parameter $\varphi_{j(i)}\in\mathbb{R}$ needs to be learned.

#### Latent alignment with domain-specific SPD batch norm.

Parametrizing the feature extractor $f_{\theta}:\mathcal{X}\rightarrow\mathcal{S}_{\tilde{P}}^{+}$ as a neural network naturally extends the decoding framework to neural networks with SPD matrix-valued features (Huang & Gool, ).
In this end-to-end learning setting, SPD batch norm (SPDBN)(Brooks et al., ; Kobler et al., ) and domain-specific batch norm (Kobler et al., ) layers can be utilized to implement $m_{\phi}$.

### 3.3 SPD Manifold Information Maximization

We utilize the labeled source domain dataset $\mathcal{D}_{s}$ to learn the shared feature extractor $f_{\theta}$, $g_{\psi}$, as well as $m_{\phi}$ for the source domains $j\in\mathcal{J}_{s}$. Specifically, we used the cross-entropy loss as the training objective, employing the PyTorch framework (Paszke et al., ) with extensions for structured matrices (Ionescu et al., ) and manifold-constrained gradients (Absil et al., ) to propagate gradients through the layers. Gradients were estimated using fixed-size mini-batches (50 observations; 10 per domain across 5 domains) and updated parameters with the Riemannian ADAM optimizer (Bécigneul & Ganea, ) ($10^{-3}$ learning rate, $10^{-4}$ weight decay, $\beta_{1}=0.9$, $\beta_{2}=0.999$). We split the source domains’ data into training and validation sets (80% / 20% splits, randomized, stratified by domain and label) and iterated through the training set for 100 epochs using exhaustive minibatch sampling.

For the target domains, we keep $f_{\theta}$ and $g_{\psi}$ fixed and learn $\phi(j)~{}\forall j\in\mathcal{J}_{t}$.
Specifically, for each domain $j\in\mathcal{J}_{t}$, we use (2) to compute the Fréchet mean $\bar{C}_{j(i)}$ of $\mathcal{C}_{j}=\{f_{\theta}(x_{i}):(x_{i},j(i))\in\mathcal{D}_{t},j(i)=j\}$.
To estimate $\Phi_{j(i)}$ in an unsupervised fashion, we propose SPD manifold information maximization (SPDIM).
SPDIM employs the information maximization (IM) loss (Shi & Sha, ) to ensure that target outputs are individually certain and globally diverse.
The IM loss is a popular training objective for SFUDA and test-time adaptation frameworks (Liang et al., ; ).

Depending on the choice of bias parameter, we distinguish between SPDIM(bias) defined in (19) and SPDIM(geodesic) defined in equation (20).
In practice, when the target domain data is revealed, we initialize $\Phi_{j(i)}$ with $I_{D}$ and $\varphi_{j(i)}$ with $1$.
We then minimize the following $\mathcal{L}_{\text{CEM}}$ and $\mathcal{L}_{\text{MEM}}$ that together constitute the
$\mathcal{L}_{\text{IM}}$ loss:

$$ $\mathcal{L}_{\text{IM}}=\mathcal{L}_{\text{CEM}}+\mathcal{L}_{\text{MEM}}$ (21) $$

$$ $\mathcal{L}_{\text{CEM}}=-\mathbb{E}_{(x_{i},j_{i})\in\mathcal{D}_{t}}\left\{ \sum_{k=1}^{|\mathcal{Y}|}\delta_{k}\left(\frac{h_{t}(x_{i},j_{i})}{T}\right) \log\delta_{k}\left(\frac{h_{t}(x_{i},j_{i})}{T}\right)\right\}$ $$

$$ $\mathcal{L}_{\text{MEM}}=\sum_{k=1}^{|\mathcal{Y}|}\hat{p}_{k}\log\hat{p}_{k}, ~{}\hat{p}_{k}=\mathbb{E}_{(x_{i},j_{i})\in\mathcal{D}_{t}}\left\{\delta_{k} \left(\frac{h_{t}(x_{i},j_{i})}{T}\right)\right\}$ $$

where $\delta_{k}$ is the k-th element of the softmax output, $\mathcal{L}_{\text{CEM}}$ is conditional entropy minimization, $\mathcal{L}_{\text{MEM}}$ is marginal entropy maximization and factor $T$ is temperature scaling.
IM balance is more effective than conditional entropy minimization because minimizing only the conditional entropy may lead to mode collapse with all test data allocated to one class (Grandvalet & Bengio, ). We additionally employed a temperature scaling factor
$T$ to adjust the model’s prediction confidence on the target data, a common technique for calibrating probabilistic models (Li et al., ; Guo et al., ). Temperature scaling uses a single scalar parameter $T>0$ for all classes, and it increases the softmax output entropy when $T>1$ and decreases it when $0<T<1$.

## 4 Experiments

We conducted simulations and experiments with public EEG motor imagery and sleep stage datasets to evaluate our proposed framework empirically.
We followed Kobler et al. () to fit the source decoder $h_{s}$.
Then, we use $h_{s}$ together with $\mathcal{D}_{t}$ to evaluate SPDIM(bias) against SPDIM(geodesic) and other baseline SFUDA methods, categorized as:

Recentering (RCT) As discussed in section 3.2, this baseline transports the latent covariance features from their domain-specific Fréchet mean to the identity matrix.

Information maximization (IM) Similar to Li et al. (), this method minimizes the $\mathcal{L}_{\text{IM}}$ loss in (21) to fine-tune all model parameters (i.e., IM(all)), $g_{\psi}$ (i.e., IM(classifier)) or merely the bias parameter of $g_{\psi}$ (i.e., IM(classifier bias)).

No domain adaptation (w/o) Methods in this category treat the problem as a standard supervised learning problem; they do not utilize the domain labels $\mathcal{J}_{s}\cup\mathcal{J}_{t}$ during training and testing.
In models that perform TSM, all data are projected to the tangent space at the Fréchet mean of the source domain data.

We used publicly available Python code for baseline methods and implemented custom methods using the packages torch (Paszke et al., ), scikit-learn (Pedregosa et al., ), braindecode (Schirrmeister et al., ), geoopt (Kochurov et al., ) and pyRiemann Barachant et al. ().
We conducted the experiments on standard computation PCs with 32-core CPUs, 128 GB of RAM, and a single GPU.

### 4.1 Simulations

To examine the effectiveness of SPDIM, we simulated binary classification problems under our generative model.
We directly generated $2\times 2$ covariance matrices $C_{i}$.
To do so, we first generated log-space features $s_{i}\in\mathbb{R}^{3}$, defined in (11), using the scikit-learn function make_classification with 1 encoding source, 2 noise sources, and varied the class separability parameter to simulate different signal-to-noise ratios.
The data were then normalized to have zero mean and unit variance.
To obtain $E_{i}\in\mathcal{S}_{D}^{+}$, we applied $\mathrm{upper}^{-1}$ and $\mathrm{Exp}_{I_{D}}$, as defined in (5).
At this level, the data were split into source domain and target domain data, with each domain containing 250 observations per class.
Finally, the data were projected to the channel space using domain-specific mixing matrices $A_{j}$, as defined in (9).
To vary the degree of label shift, we artificially varied the label ratio, defined as the proportion of the minority class to the majority class, in the target domain via randomly dropping samples.
We used balanced accuracy as the evaluation metric and examined the performance of each method with different class separability under different imbalance ratios.

The results are summarized in Figure 2.
The simulation results demonstrate that as the label shifts become more severe (i.e., lower label ratio), all methods suffer from performance decreases.
Still, SPDIM (bias) outperforms the considered baseline approaches in maintaining a higher score for the same label ratio and signal-to-noise ratio (i.e., class separability).
Even in severe cases (0.1 label ratio), SPIM (bias) improves upon RCT by more than 5% for moderate to high signal-to-noise ratios.
Figure 1 b-d visualizes the representations obtained with RCT (top) and SPDIM (bottom) for a representative example.

Figure: Figure 2: SPDIM simulation results for 2 classes, $2\times 2$ covariance features and 1 source and 1 target domain. Balanced accuracy scores (higher is better) across target domain label ratios (majority to minority class) on the x-axis and class separability $p(s|y)$ across panels. Source domain labels were balanced.
Refer to caption: x1.png

### 4.2 EEG Motor Imagery Data

We considered 4 public motor imagery datasets: BNCI2014001 (Tangermann et al., ) (9 subjects/2 sessions/4 classes/22 channels), BNCI2015001 (Faller et al., ) (12/2-3/2/13), Zhou2016 (Zhou et al., ) (4/3/3/14), and BNCI2014004 (Leeb et al., ) (9/5/2/3).
We used MOABB (Jayaram & Barachant, ; Chevallier et al., ) to pre-process the continuous time-series data and extract labeled epochs.
Pre-processing included resampling EEG signals to 250 or 256 Hz, applying temporal filters to capture frequencies between 4 and 36 Hz, and extracting 3-second epochs linked to specific class labels.
Following Kobler et al. (), we use TSMNet as the model and treat sessions as domains, and use a leave-one-group-out cross-validation (CV) scheme to fit and evaluate the methods.
To evaluate cross-session transfer, we fitted and evaluated models independently per subject and treated the session as the grouping variable.
To evaluate cross-subject transfer, we treated the subject as the grouping variable.
After running pilot experiments with the BNCI2014001 dataset, we set the temperature scaling factor in (21) to $T=2$ for binary classification problems and $T=0.8$ otherwise.
Other hyper-parameters (e.g., early stopping) were fit with a single stratified (domain and labels) inner train/validation split.
As in the simulation experiment, we artificially introduced label shifts in the target domains.
We used balanced accuracy as the metric to examine the performance of each method at label ratios of 1.0 and 0.2 for the source and target domains, respectively.

Grand average results across all 4 datasets are summarized in Figure 3, and grouped by the transfer learning scenario (cross-session, cross-subject).
To attenuate large variability across subjects, we report scores relative to the score obtained with RCT w/o label shifts (i.e., 1.0 label ratio).
Although no method can perfectly compensate the artificially introduced label shifts, SPDIM(bias) is consistently at the top (cross-subject) or among the top (cross-session) performing methods.
Significance testing (n=34 subjects), summarized in Supplementary Table A3, revealed that performance for SPDIM(bias) is significantly higher than w/o, RCT, IM(classifier), and IM(all) in the cross-subject setting as well as RCT, IM(classifier), and IM(all) in the cross-session scenario.

Figure: Figure 3: Motor-imagery results (0.2 label ratio). Average of test-set scores (balanced accuracy; higher is better; standard-deviation in brackets) relative to TSMNet+RCT (Kobler et al., ) applied to balanced target domain data (i.e., 1.0 imbalance ratio). For extended results per dataset, see Supplementary Table A3.
Refer to caption: x2.png

### 4.3 EEG-based Sleep Staging

The aim of this experiment is to demonstrate the effectiveness of our method on datasets with inherent label shifts.
Due to the inherent variability of sleep, most sleep stage datasets exhibit label shifts across subjects (Eldele et al., ).
Sleep stage classification plays a key role in assessing sleep quality and diagnosing sleep disorders (Perez-Pozuelo et al., ).
Yet, automated frameworks lack accuracy and suffer from poor generalization across domains, resulting in accuracy drops compared to expert neurologists.

We considered 4 public sleep stage datasets: CAP (Terzano et al., ; Goldberger et al., ), Dreem (Guillot et al., ), HMC (Alvarez-Estevez & Rijsman, ; ), and ISRUC (Khalighi et al., ).
A detailed description is provided in Supplementary Table A1.
We consider sleep stages following the AASM (Berry et al., ) standard (W, N1, N2, N3, REM).
If data was originally scored following the K&M (Wolpert, ) standard (W, S1, S2, S3, S4, REM), we merged stages S3 and S4 into a single stage N3.
EEG data pre-processing followed (Guillot & Thorey, ) and was implemented with MNE-python (Gramfort et al., ).
First, all EEG channels were retained, and an IIR band-pass filter ranging from 0.2 to 30 Hz was applied.
The signals were then resampled to 100 Hz, and non-overlapping 30-second epochs were extracted together with the associated labels.
To attenuate the effects of gross outliers, each recording was scaled to have unit inter-quartile range and a median of zero, and values exceeding 20 times the inter-quartile range were clipped Perslev et al. ().
Lastly, subjects with corrupted data (e.g., mismatched labels and epochs) were excluded, resulting in a total of 426 remaining subjects.

Since label shifts occur in the source domains of sleep stage data, the considered models are trained with a balanced mini-batch sampler, which is a popular method to compensate for label shifts during training in deep learning (Cao et al., ).
The balanced sampler over-sampled minority classes to ensure that the label distribution per domain is balanced within each mini-batch.

To evaluate the methods, we employed a 10-fold grouped cross-validation scheme, ensuring that each group (i.e., subject) appears either in the training set (i.e., source domains) or the test set (i.e., target domains).
As before, we set the temperature scaling factor $T=0.8$ and used stratified (labels and domains) inner train/validation splits for early stopping.
In addition to the baseline methods introduced before, we included two additional benchmark deep learning architectures that were specifically proposed for sleep staging: Chambon (Chambon et al., ) and Usleep (Perslev et al., ).
For both architectures, we use the implementation provided in braindecode (Schirrmeister et al., ).

Table 1 summarizes the results across datasets along with the grand average results of published baseline methods.
Extended results for all considered baseline methods are listed in Supplementary Table A2.
Interestingly, TSMNet+RCT (Kobler et al., ) obtained similar or higher results compared to non-SFUDA baselines.
TSMNet+SPDIM(bias) significantly outperforms all other methods for the patient and healthy subject groups.
Overall, the margin to TSMNet+RCT was approx. 5% in the patient group, which indicates that our proposed method has great potential for clinical applications.

Table: Table 1: Sleep-staging results per dataset. Average of test-set scores (balanced accuracy; higher is better; standard-deviation in brackets). Permutation-paired t-tests were used to identify significant differences between TSMNet+SPDIM (*proposed*) and baseline methods (1e4 permutations, 4 tests, t-max correction). Significant differences are highlighted ($\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}$ $p\leq 0.05$, $\mathbin{\vbox{\hbox{\scalebox{0.75}{$\bullet$}}}}$ $p\leq 0.01$, $\mathbin{\vbox{\hbox{\scalebox{1.0}{$\bullet$}}}}$ $p\leq 0.001$). Extended results are provided in Supplementary Table A2.

Ablation Study
Table 2 summarizes grand average test scores relative to SPDIM(bias).
We highlight three observations.
First, all considered ablations lead to a significant performance drop of at least 3% compared to SPDIM(bias), suggesting the combined importance of IM paired with the manifold-constrained bias parameter and RCT.
Second, fine-tuning the classifier bias parameter yields approximately the same performance as SPDIM(geodesic), extending the finding of (Mellot et al., ) from regression to classification scenarios.
Third, IM SFUDA methods obtained the top 3 scores, but two variants (IM(classifier) and IM(full)) failed to improve upon the non-SFUDA method, which underscores the importance of selecting the right parameter to prevent the IM loss from over-fitting.

**Table 2: Sleep-staging ablation results. Balanced accuracy scores (higher is better) relative to the proposed method for the combined sleep staging data. Averages and standard deviation summarize the individual (n=426 subjects) test-set scores. Student’s t values and adjusted p values indicate the effect strength (permutation-paired t-tests, 1e4 permutations, 6 tests, t-max correction).**
| SFUDA | $\Delta$ balanced accuracy (%) |  |
| --- | --- | --- |
|  | mean (std) | t-val (p-val) |
| SPDIM(bias) | - | - |
| SPDIM(geodesic) | -3.0 (3.6) | -17.4 (0.0001) |
| RCT | -4.8 (4.5) | -22.1 (0.0001) |
| w/o | -3.7 (8.5) | -9.0 (0.0001) |
| IM(classifier bias) | -3.0 (3.8) | -16.1 (0.0001) |
| IM(classifier) | -4.8 (5.4) | -18.5 (0.0001) |
| IM(all) | -27.1 (13.5) | -41.5 (0.0001) |

## 5 Discussion

We proposed a geometric deep learning framework, denoted SPDIM, to address SFUDA problems with conditional and label shifts and demonstrated its utility in highly relevant EEG-based neurotechnology application scenarios.
We first introduced a realistic generative model, and provided theoretical analyses showing that prior Riemannian statistical alignment methods that align the Fréchet mean can compensate for the conditional distribution shifts introduced in our generative model, but hurt generalization under additional label shifts.
As a remedy, we proposed SPDIM to learn a domain-specific SPD manifold-constrained bias parameter to compensate over-corrections introduced via aligning the Fréchet means by employing the information maximization principle.
In simulations and experiments with real EEG data, SPDIM consistently achieved the highest scores among the considered baseline methods.

## 6 Acknowledgments

Motoaki Kawanabe and Reinmar J Kobler were partially supported by the Innovative Science and Technology Initiative for Security Grant Number JPJ004596, ATLA, Japan.

## 7 Author Contributions

Shanglin Li contributed to the study under the co-supervision of Motoaki Kawanabe and Reinmar J. Kobler.
Shanglin Li and Reinmar J. Kobler developed the decoding methods.
Reinmar J. Kobler contributed the theoretical analysis.
Shangling Li and Reinmar J. Kobler performed the simulation experiments.
Shanglin Li conducted the experiments with EEG data.
The draft of the manuscript was written by Shanglin Li and Reinmar Kobler.
All authors read and approved the final manuscript.

## Appendix A Appendix

### A.1 Proof of Proposition 1

Proposition 1 stats that given the generative model specified in section 3.1 and a set of examples $\mathcal{E}_{j}=\{(E_{i},y_{i},j)|E_{i}\in\mathcal{S}_{P}^{+},y_{i}\in\mathcal
{Y}\}_{i\leq M_{j}}$ of domain $j\in\mathcal{J}$, we have that the Fréchet mean of $\mathcal{E}_{j}$, defined in (2), converges with $M_{j}\rightarrow\infty$ to the identity matrix $I_{P}$ for all domains $j\in\mathcal{J}$.

Since $(\mathcal{S}^{+}_{P},g_{\cdot}^{\text{AIRM}})$ forms a Cartan-Hadamard manifold with global non-positive sectional curvature, a unique Fréchet mean exists (Bhatia, ). At the global minimum, we have:

$$ $\displaystyle 0\overset{!}{=}$ $\displaystyle\nabla_{E}\nu_{E}(\mathcal{E}_{j})=\nabla_{E}\frac{1}{M_{j}}\sum_ {i=1}^{M_{j}}\delta^{2}(E,E_{i})=\frac{1}{M_{j}}\sum_{i=1}^{M_{j}}\nabla_{E} \delta^{2}(E,E_{i})=-\frac{1}{M_{j}}\sum_{i=1}^{M_{j}}\mathrm{Log}_{E}(E_{i})$ (22) $\displaystyle=$ $\displaystyle-\frac{1}{M_{j}}\sum_{i=1}^{M_{j}}E^{\frac{1}{2}}\log(E^{-\frac{1 }{2}}E_{i}E^{-\frac{1}{2}})E^{\frac{1}{2}}=-\frac{1}{M_{j}}E^{\frac{1}{2}} \left(\sum_{i=1}^{M_{j}}\log(E^{-\frac{1}{2}}E_{i}E^{-\frac{1}{2}})\right)E^{ \frac{1}{2}}$ (23) $\displaystyle 0\overset{!}{=}$ $\displaystyle\sum_{i=1}^{M_{j}}\log(E^{-\frac{1}{2}}E_{i}E^{-\frac{1}{2}})$ (24) $$

where we used the derivative of the Riemannian distance (Moakher, ; Pennec, ) in (22).
For $E=I_{P}$, (24) simplifies to:

$$ $\displaystyle\sum_{i=1}^{M_{j}}\log(E_{i})\overset{!}{=}0\underset{M_{j} \rightarrow\infty}{\Longrightarrow}\mathbb{E}\{\log(E_{i})\}$ $\displaystyle\overset{!}{=}0$ (25) $\displaystyle\mathbb{E}\{B\tilde{y}_{i}+\varepsilon_{i}\}=B\mathbb{E}\{\tilde{ y}_{i}\}+\mathbb{E}\{\varepsilon_{i}\}$ $\displaystyle\overset{!}{=}0$ (26) $$

which holds true because by definition $\tilde{y}_{i}$ and $\varepsilon_{i}$ are zero-mean.
Therefore, $E_{j}=I_{P}$ is the global minimizer of $\nu_{E}(\mathcal{E}_{j})$ for all domains $j\in\mathcal{J}$. $\square$

### A.2 Proof of Proposition 2

Proposition 2 states that for the generative model, specified in section 3.1, RCT+TSM compensates conditional distribution shifts introduced by the invertible linear map $A_{j}$, defined in (8), and recovers domain-invariant representations if there are no target shifts (i.e., $p_{j}(\mathrm{y}{=}y)=p(\mathrm{y}{=}y)~{}\forall y\in\mathcal{Y},~{}\forall j
\in\mathcal{J}_{s}\cup\mathcal{J}_{t}$).

Due to the congruence invariance of the Fréchet mean (Bhatia, ) we have

$$ $\bar{C}_{j}=\arg\min_{G\in\mathcal{S}^{+}_{P}}\nu_{G}(\mathcal{C}_{j})=A_{j} \bar{E}_{j}A_{j}^{T}$ (27) $$

where $\bar{E}_{j}$ is the Fréchet mean of $\mathcal{E}_{j}$.
Utilizing proposition 1 (i.e., $\bar{E}_{j}=I_{P}~{}\forall\,j$) and plugging in (9) for $A_{j}$, (27) simplifies to:

$$ $\displaystyle\bar{C}_{j}=A_{j}I_{P}A_{j}^{T}=Q\exp(P_{j})\exp(P_{j})Q^{T}$ $\displaystyle=Q\exp(P_{j})^{2}Q^{T}$ (28) $\displaystyle\Rightarrow\bar{C}_{j}^{-\frac{1}{2}}$ $\displaystyle=Q\exp(P_{j})^{-1}Q^{T}$ (29) $$

The RCT+TSM transform, as defined in (14), recenters the data $C_{i}\in\mathcal{C}_{j}$ for each domain $j$.
Excluding the invertible mapping $\mathrm{upper}$, RCT+TSM computes:

$$ $\displaystyle\log\left(\bar{C}_{j}^{-\frac{1}{2}}C_{i}\bar{C}_{j}^{-\frac{1}{2 }}\right)$ $\displaystyle=\log\left(Q\exp(P_{j})^{-1}Q^{T}C_{i}Q\exp(P_{j})^{-1}Q^{T}\right)$ (30) $\displaystyle=\log\left(Q\exp(P_{j})^{-1}A_{j}E_{i}A_{j}^{T}\exp(P_{j})^{-1}Q^ {T}\right)$ (31) $\displaystyle=\log\left(Q\exp(P_{j})^{-1}Q^{T}Q\exp(P_{j})E_{i}\exp(P_{j})Q^{T }Q\exp(P_{j})^{-1}Q^{T}\right)$ (32) $\displaystyle=\log\left(QE_{i}Q^{T}\right)=Q\log(E_{i})Q^{T}$ (33) $$

where we used (13) for $C_{i}$, the fact that $Q$ is an orthogonal matrix and $\log(ACA^{-1})=A\log(C)A^{-1}$ for non-singular $A$ and $C\in\mathcal{S}_{P}^{+}$.
Plugging (11) for $\log(E_{i})$, we obtain a direct relationship to the label $y_{i}$ via $\tilde{y}_{i}$

$$ $\displaystyle\log\left(\bar{C}_{j}^{-\frac{1}{2}}C_{i}\bar{C}_{j}^{-\frac{1}{2 }}\right)$ $\displaystyle=Q\mathrm{upper}^{-1}\left(B\tilde{y}_{i}+\varepsilon_{i}\right)Q ^{T}$ (34) $\displaystyle=Q\mathrm{upper}^{-1}\left(B\left(1_{y_{i}}-\bar{y}_{j}\right)+ \varepsilon_{i}\right)Q^{T}$ (35) $$

where the vector $\bar{y}_{j}=[p_{j}(\mathrm{y}{=}y_{k})]_{k=1,...,|\mathcal{Y}|}$ contains the class priors for domain $j$.
If there are no target shifts the class priors are constant $p_{j}(\mathrm{y}{=}y)=p(\mathrm{y}{=}y)~{}\forall y\in\mathcal{Y}$ for all domains $j\in\mathcal{J}_{s}\cup\mathcal{J}_{t}$, and consequently $\bar{y}_{j}=\bar{y}$.
Then (35) simplifies to:

$$ $\log\left(\bar{C}_{j}^{-\frac{1}{2}}C_{i}\bar{C}_{j}^{-\frac{1}{2}}\right)=Q \mathrm{upper}^{-1}\left(B\left(1_{y_{i}}-\bar{y}\right)+\varepsilon_{i}\right )Q^{T}$ (36) $$

which contains only domain-invariant terms on the right hand side.
Thus, RCT+TSM compensates the conditional shifts introduced by $A_{j}$. $\square$

### A.3 Sleep stage dataset details

**Table A1: Sleep stage dataset details. Overview of datasets and their subgroups (SG), including the number of recordings, subjects, channel numbers, patient status, scorer (the chosen scorer if there are multiple scorers), and scoring rule.**
| Dataset | Recordings | Subjects | Channel numbers | Patients | Scorer | Scoring rule |
| --- | --- | --- | --- | --- | --- | --- |
| ISRUC-SG1 | 100 | 100 | 6 | ✓ | Scorer1 | AASM |
| ISRUC-SG2 | 16 | 8 | 6 | ✓ | Scorer1 | AASM |
| ISRUC-SG3 | 10 | 10 | 6 | ✗ | Scorer1 | AASM |
| Dreem-SG1 | 22 | 22 | 12 | ✗ | Scorer1 | AASM |
| Dreem-SG2 | 50 | 50 | 8 | ✓ | Scorer1 | AASM |
| HMC | 154 | 154 | 4 | ✓ | - | AASM |
| CAP-SG1 | 36 | 36 | 13 | ✓ | - | K&M |
| CAP-SG2 | 34 | 34 | 9 | ✓ | - | K&M |
| CAP-SG3 | 22 | 22 | 5 | ✓ | - | K&M |

### A.4 EEG-based Sleep Staging Results

Table: Table A2: Sleep-staging results per dataset. Average and standard deviation of test-set scores (balanced accuracy; higher is better) across public sleep staging datasets. For all IM and SPDIM variants, the parameters that were tuned to the test-data with the IM loss are indicated in brackets. Permutation-paired t-tests were used to identify significant differences between our *proposed* (i.e., TSMNet+SPDIM(bias)) and baseline methods (1e4 permutations, 8 tests, t-max correction). Student’s t values summarize the effect strength. Significant differences are highlighted ($\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}$ $p\leq 0.05$, $\mathbin{\vbox{\hbox{\scalebox{0.75}{$\bullet$}}}}$ $p\leq 0.01$, $\mathbin{\vbox{\hbox{\scalebox{1.0}{$\bullet$}}}}$ $p\leq 0.001$).

### A.5 EEG-based Motor Imagery BCI

Table: Table A3: Motor imagery BCI results per dataset. Average and standard deviation of test-set scores (balanced accuracy; higher is better) across public motor imagery BCI datasets. For all IM and SPDIM variants, the parameters that were tuned to the test-data with the IM loss are indicated in brackets. Permutation-paired t-tests were used to identify significant differences between the *proposed* (i.e., TSMNet+SPDIM(bias)) and baseline methods (1e4 permutations, 7 tests, t-max correction). Student’s t values summarize the effect strength. Significant differences are highlighted ($\mathbin{\vbox{\hbox{\scalebox{0.5}{$\bullet$}}}}$ $p\leq 0.05$, $\mathbin{\vbox{\hbox{\scalebox{0.75}{$\bullet$}}}}$ $p\leq 0.01$, $\mathbin{\vbox{\hbox{\scalebox{1.0}{$\bullet$}}}}$ $p\leq 0.001$).

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 8
presentation: 4
questions: - Q: Does this framework treat one subject or one EEG recording as one source/target domain containing both/multiple class labels?
- Q: How does this framework for "latent space alignment" compare/relate to non-reimannian approaches for SFUDA for EEGs/multivariate timeseries? See [1] for a recent example. The "test-time adaptation" (Section 3.2) studies listed in [2] might also be relevant.
- Q: What factors other than dataset size and label shifts could account for the high variability/stdev in Table 1? In most cases, handling label shift (either with RCT or SPDIM) decreases variability compared to "w/o", but its still seems high.
- Minor comments: 1) pixel resolution of Figure 1 can be improved, 2) typo in citations at line 218 and 236., 3) line 443 remove "standard-deviation in brackets"
- The anonymous code link is broken?

[1] He, Huan, et al. "Domain adaptation for time series under feature and label shifts." International Conference on Machine Learning. PMLR, 2023.

[2] Garg, Saurabh, et al. "Rlsbench: Domain adaptation under relaxed label shift." International Conference on Machine Learning. PMLR, 2023.
soundness: 3
strengths: - Rigorous and clear presentation of technical details and full analytic workflow.
- This work is a great example of theory-guided methods design for EEG.
- Impactful choice of research problem - performance of EEG models under label shifts will remain a ubiquitous concern, both clinically and in the BCI space.
summary: This study focuses on the realistic issue of label shifts in EEG across subjects and/or sessions (relative class proportions in target domains when source domains are class-balanced). Using theoretical analysis, it extends the SotA statistical alignment framework for handling distribution shifts in EEG to also include label shifts. The proposed SPDIM includes a domain-specific bias parameter estimated from unlabeled target data that reduces over-corrections done by the current SotA framework. Results on synthetic data and real-world EEGs demonstrate the value of SPDIM over baselines.
weaknesses: - (line 166) Q: Is the assumption of number of latent brain sources = number of observed scalp channels = P necessary or realistic?
- No discussion of study limitations and/or future directions.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: - Is the domain-specific formard model $A_{j}$ learned from features of a specific domain, or is it predefined? 
- How are domain specific parameters $\Phi_{j(i)}$ and the geodesic step-size parameters $\varphi_{j(i)}$ learned according to the proposed information maximization principle described in Section 3.3?
- Is there any relationship between the adaptation performance and predefined hyperparameters, such as the rank of $A$ and the number of domains within $\mathcal{D}_{s}$?
soundness: 3
strengths: - The motivation is clear and easy to follow. SPDIM aims to address adaptation under label shifts, a common challenge in real-world EEG datasets. 
  Theoretical analysis further explains the causes of deviations under label shifts.

- Simulation experiments qualitatively validate the benefits of SPDIM in the presence of label shifts. 
   Cross-subject and cross-session experiments on motor and sleep-staging EEG datasets illustrate its superiority over existing alignment methods based on the SPD manifold.
summary: The study addresses a source-free unsupervised domain adaptation problem and proposes SPDIM, a framework based on the SPD manifold. 
SPDIM compensates for label shifts using proposed generative models, which prior Riemannian statistical alignment methods do not effectively handle. 
Additionally, SPDIM applies the information maximization principle to learn domain-specific parameters. 
Simulation experiments demonstrate its superiority under various levels of label shift, and empirical analysis on real EEG datasets shows that it outperforms previous approaches.
weaknesses: - Some notations in equations seem confusing. For example, the index $j$ under $\sum$ may need to be $i$ in Eq. (2). The invertible mapping $upper$ is defined on $S$, but $upper^{-1}$ appears in Eq.(10). 
Additionally, $j_i$ and $j$ use the same letter but with different meanings, which could lead to ambiguity. The notation $Q$ in Eq.(15) seems to appear without prior introduction. 

- Some aspects of the method require further clarification. As mentioned in Line 249, the right-hand side of Eq. (15) is claimed to contain only domain-invariant terms. However, from my perspective, $C_i$ depends on the domain-specific 
matrix $A_{j}$, as suggested by Eq. (13). According to Proposition 1, $ \bar{C} _ {j(i)} $ converges to $I_P$. These indicate that $Q$ is linked to $A_{j}$, which may not be domain-invariant. Additionally, the relationship between the information maximization approach introduced in Section 3.3 and SPDIM (bias) / SPDIM (geodesic) is unclear. 

- To better demonstrate SPDIM’s effectiveness, it would be beneficial to compare it with additional statistical alignment methods beyond those based on the SPD manifold. This would provide a more comprehensive evaluation against existing approaches.

---

confidence: 2
contribution: 3
initial_rating: 5
presentation: 1
questions: Plz go and check weaknesses
soundness: 3
strengths: 1.	The introduction of an SPD-manifold-constrained bias parameter is an advancement for tackling SFUDA in EEG.
2.	The framework has been applied effectively across different tasks, showcasing broad applicability.
3.	SPDIM outperforms conventional methods, showing its resilience under varying label distributions.
summary: The "SPDIM" paper introduces a novel geometric deep learning framework aimed at enhancing source-free unsupervised domain adaptation (SFUDA) for EEG data under both conditional and label shifts. By leveraging the symmetric positive definite (SPD) manifold and employing a parameter-efficient manifold optimization strategy, the proposed method, SPDIM, addresses significant generalization challenges in EEG data processing, especially where traditional Riemannian geometry methods fall short due to label shifts. SPDIM shows promising improvements across multiple EEG-based applications, including brain-computer interface tasks and sleep staging, demonstrating its efficacy over prior alignment frameworks
weaknesses: 1.	The motivation behind addressing label shifts and domain gaps with SPDIM is somewhat implicit, without clearly laying out why these challenges necessitate the proposed framework. 
2.	The paper contains an extensive number of equations and mathematical formulations in the main text, which can make the methodology difficult to follow. 
3.	Although the paper compares SPDIM with several baselines, a broader set of comparisons, especially with newer unsupervised or semi-supervised EEG methods, could provide further insights into SPDIM’s performance and robustness.
4.	While SPDIM improves accuracy under domain shifts, the model’s interpretability remains limited. 

I will reconsider my assessment after checking rebuttual.

---

confidence: 5
contribution: 3
initial_rating: 5
presentation: 3
questions: - You mention there are conditional shifts in EEG data (p_j(x|y) changes between domains). Can you relate this with your modelization?
- What is D in the Remark 1?
- Does the Propostion 2 still hold when M_j does not tend to the infinite?
- You train your model on the target domain (in an unsupervised manner). Did you train/test split the target domain?
- How easy to train are the methods you use? e.g. USleep is rarely used as a baseline in other sleep staging papers. Providing infos the lr scheduler, batch size, … would be valuable.
- I am surprised that the spatial covariance is enough to classify sleep stages. Usually, the temporal information is used but not the spatial one. Can you comment on this?

A few typos:
- D and P are both used for the data dimension
- There are “?” in lines 218 and 236.
- Q and U are both used for domain-invariant par of the mixing matrix.
soundness: 3
strengths: - The paper is well written.
- The modelisation is original and insightful. I really enjoyed reading the modelisation part of the paper.
- The developed methods are tested on 3 setups: synthetic data, Motor-Imagery and Sleep-Staging.
summary: The paper introduces SPDIM, a framework for source-free unsupervised domain adaptation (SFUDA) in EEG-based applications, which are challenged by distribution shifts across sessions or subjects. SPDIM leverages the geometry of symmetric positive definite (SPD) matrices to handle conditional and label shifts, aligning EEG data across domains without requiring labeled target data. The approach introduces a domain-specific SPD-manifold bias to counteract label shifts, and optimizes alignment using an information maximization principle, which prevents mode collapse by ensuring class diversity and prediction confidence. Experimental results on EEG datasets for motor imagery and sleep staging show that SPDIM outperforms baseline SFUDA methods, demonstrating robust generalization across domains even under significant label distribution changes.
weaknesses: - At the time of reviewing the paper, the code is not available: “The repository is not found.” is returned by anonymous.4open.science
- A modelisation per domain of EEG data was proposed in [1] which could be worth citing in your introduction. Indeed, the authors mention there exists a linear mapping per domain to get domain-invariant tangent vectors (and without assumption on the mixing matrix (9)).
- The experiment on motor imagery is limited since you artificialy unbalance the labels. Finding real world data which are naturally unbalanced would add value to the paper.
- The mean accuracy of the 2 proposed methods are within the standard deviation of the recenter for the motor imagery application.
- On the sleep-staging setup, you do not compare with adaptation methods expect recenter. You should compare at least to STMA or TMA (Spatio-Temporal Monge Alignment) which is presented in [2].
- The presentation of the results is not homogeneous between the two applications. In particular, it is strange to me to call an “ablation study” a comparison with other methods.

[1] Collas, Antoine, Rémi Flamary, and Alexandre Gramfort. "Weakly supervised covariance matrices alignment through Stiefel matrices estimation for MEG applications." arXiv preprint arXiv:2402.03345 (2024).

[2] Gnassounou, T., Collas, A., Flamary, R., Lounici, K., & Gramfort, A. (2024). Multi-Source and Test-Time Domain Adaptation on Multivariate Signals using Spatio-Temporal Monge Alignment. arXiv preprint arXiv:2407.14303.


I put a rating of 5 but I am open to increasing it.

---

confidence: 3
contribution: 2
initial_rating: 6
presentation: 3
questions: Please see the weakness.
soundness: 3
strengths: This paper investigates the label shift problem in SFUDA, with a very strong motivation and significant practical relevance.

The method proposed in this paper has a certain theoretical foundation, and the derivation of some propositions may provide inspiration for solving the label shift problem.
summary: Due to the difficulty of the SPD method in handling label shift issues, this paper proposes a geometric deep learning framework, SPDIM, for SFUDA problems under specific distribution shifts, including label shifts. SPDIM employs the information maximization principle to learn a single SPD-manifold-constrained parameter per target domain. Using public EEG-based brain-computer interface and sleep staging datasets, we demonstrate that SPDIM outperforms prior approaches.
weaknesses: **The clarifications and revisions have addressed most of my concerns.**

---

Although this paper focuses on EEG SFUDA problems, the proposed method does not appear to be specifically designed for EEG but seems to be a more general approach applicable to any label shift scenario. From the perspective of EEG research, the method lacks specificity for EEG data, while from the perspective of SFUDA research, the paper only validates the method on EEG data, lacking more reliable experimental verification.

The experiments are not solid. The paper does not clearly present the experimental setup, such as the hyperparameters of the models, the partitioning method of the source and target domains, etc. Additionally, the EEG decoding methods compared in the experiments are not sufficiently strong. The paper does not compare some classic EEG decoding models, such as EEGNet and EEG Conformer, nor does it compare some sleep staging models, such as DeepSleepNet. The domain adaptation methods only compare Information Maximization (IM), and such insufficient comparisons are not enough to prove the superiority of the proposed method.

The writing of this paper still has some room for improvement. For example: Figure 1 has low resolution, and the four sub-figures in Figure 2 lack sub-titles.

Meta-review: This paper was considered "Rigorous" with "clear presentation", and "a great example of theory-guided methods design for EEG". by reviewer PCcy and also endorsed by reviewers CvWN and ocKU, in particular after a valuable discussion with ocKU that surely helped clarify the paper and increase it's experimental part.

Based on reviews and further discussions, this paper is considered a relevant and good contribution for the ICLR community and in particular the ML researchers working on neural signal decoding.

Decision: Accept (Poster)