Deep Optimal Transport for Domain Adaptation on SPD Manifolds

A matrix $S\in\mathbb{R}^{n\times n}$ is symmetric and positive-definite if it is symmetric and all its eigenvalues are positive, i.e., $S=S^{T}$, and $x^{T}\cdot S\cdot x>0$, for all $x\in\mathbb{R}^{n}\setminus\{{\mathbf{0}}\}$. The Frobenius inner product and Frobenius norm on $m\times n$ matrices $A$ and $B$ are defined as $\langle\,A,B\,\rangle_{\mathcal{F}}:=\operatorname{Tr}{(A^{\top}\cdot B)}$ and $||\,A\,||_{\mathcal{F}}^{2}:=\langle\,A,A\,\rangle_{\mathcal{F}}$ respectively.

The Log-Euclidean metric, a widely used Riemannian metric on SPD manifolds, introduces a vector space structure to the space of SPD matrices, endowing it with numerous exceptional properties [14, 15]. These include invariance under inversion, logarithmic multiplication, orthogonal transformations, and scaling. Formally, a commutative Lie group structure, known as the tensor Lie group, is endowed on $\mathcal{S}_{++}$ using the logarithmic multiplication $\odot$ defined by $P_{1}\odot P_{2}:=\exp\big{(}\log P_{1}+\log P_{2}\big{)}$, where $P_{1},P_{2}\in\mathcal{S}_{++}$, and $\exp$ and $\log$ are matrix exponential and logarithm respectively. The bi-invariant metric exists on the tensor Lie group, which yields the following scalar product $g^{LEM}_{P}(v,w):=\langle\,D_{P}\log v,D_{P}\log w\,\rangle_{\mathcal{F}}$, where tangent vectors $v,w$ at matrix $P$, and $D_{P}\log$ represent the differential $D\log$ at $P$. We denote an SPD manifold equipped with the Log-Euclidean metric as ($\mathcal{S}_{++}$, $g^{LEM}$). The sectional curvature of ($\mathcal{S}_{++}$, $g^{LEM}$) is flat, and thus parallel transport on it simplifies to an identity. Riemannian distance $d_{g^{LEM}}(P_{1},P_{2})=||\log{P_{1}}-\log{P_{2}}||_{\mathcal{F}}$ between $P_{1}$ and $P_{2}$ on this SPD manifold. The Log-Euclidean Fr$\acute{e}$chet mean $\overline{w}(\mathcal{B})$ of a batch $\mathcal{B}$ of SPD matrices $\{S_{i}\}_{i=1}^{|\mathcal{B}|}$ is defined as follows:

In this study, another Riemannian metric on SPD manifolds, affine-invariant Riemannian metric [16], is introduced for comparison with $g^{LEM}$. SPD manifold equipped with this metric is denoted as ($\mathcal{S}_{++}$, $g^{AIRM}$). Formally, for any $P$, it is defined as $g^{AIRM}_{P}(v,w):=\langle\,P^{-\frac{1}{2}}\,v\,P^{-\frac{1}{2}},P^{-\frac{1}{2}}\,w\,P^{-\frac{1}{2}}\,\rangle_{\mathcal{F}}$, where tangent vectors $v,w\in\mathcal{T}_{P}\mathcal{S}_{++}$. ($\mathcal{S}_{++}$, $g^{AIRM}$) also exhibit nice properties, such as affine invariance, scale invariance, and inversion invariance. Since these properties are not the primary focus of this study, please refer to details in [17].

Optimal Transport (OT) seeks the minimal cost transport plan between source and target distributions [18]. In this study, we consider the case where these distributions are available only through a finite number of samples. Given two measurable spaces $X$ and $\bar{X}$, an empirical distribution of the measure $\mu\in X$ and $\nu\in\bar{X}$ are written as $\mu:=\sum_{i\in I}p_{i}\delta_{x_{i}}$ and $\nu:=\sum_{j\in J}q_{i}\delta_{\bar{x}_{i}}$, where discrete variable $X$ takes values $\{x_{i}\}_{i\in I}\in X$, countable index set $I$ is with weights $\sum_{i\in I}p_{i}=1$. $\bar{X}$ takes values $\{\bar{x}_{j}\}_{j\in J}\in\bar{X}$, countable index set $J$ is with weights $\sum_{j\in J}q_{j}=1$. Dirac measure $\delta_{s}$ denotes a unite point mass at the point $s\in X\text{ or }\bar{X}$. The discrete Kantorovich formalization of the OT problem is given as $\min_{\gamma\in\prod(\mu,\nu)}\big{\langle}\gamma,c(x,\bar{x})\big{\rangle}_{\mathcal{F}}$, where transportation plan $\gamma\in\prod(\mu,\nu):=\{\gamma\in\mathbb{R}_{\geq 0}^{|X|\times|\bar{X}|}\big{|}\gamma\mathbbm{1}_{|X|}=\mu,\text{ and }\gamma^{\top}\mathbbm{1}_{|\bar{X}|}=\nu\}$, and $\mathbbm{1}_{d}$ is a $d-$dimensional all-one vector.

Assuming that the transfer between source and target distributions can be estimated using optimal transport, Joint Distribution Optimal Transport (JDOT) [8] extends the OT-DA framework by jointly optimizing both the optimal coupling and the prediction function $f$ for the two distributions, i.e., $\min_{\gamma\in\prod(\mu,\nu)}\big{\langle}\gamma,c^{\prime}(x,y;\,\bar{x},f(\bar{x}))\big{\rangle}_{\mathcal{F}}$, where $y$ is the label set for samples $x$, $f$ is the prediction function, and the joint cost function is a linear combination of distance function and cross-entropy loss, i.e., $c^{\prime}(x,y;\bar{x},f(\bar{x})):=\alpha_{1}d_{\ell^{2}}(x,\bar{x})+\alpha_{2}\mathcal{L}(y,f(\bar{x}))$. Moreover, DeepJDOT employs neural networks to map raw data into feature space and thus extends JDOT by setting the new joint objective function $\mathcal{L}:=\mathcal{L}(y,f(g(x)))+\big{\langle}\gamma,c^{\prime}\big{(}g(x),y;\,g(\bar{x}),f(g(\bar{x})\big{)}\big{\rangle}_{\mathcal{F}}$, where $g$ represents the neural networks [9].

Extensive research has explored various classes of Riemannian metrics on SPD manifolds for various applications. These explorations have led to the introduction of numerous statistical methods to address computational issues on SPD manifolds, drawing inspiration from manifold-valued neuroimaging data [19, 20]. During the feature engineering era, Barachant et al. utilized the Riemannian distance between EEG spatial covariance matrices for motor imagery classification, resulting in a significant impact on various brain-computer interface applications [21, 22, 23].

In the deep learning era, Ju et al. employed geometric deep learning on SPD manifolds to extract discriminative information from the time-space-frequency domain for motor imagery classification, ultimately enhancing classification performance [24, 25, 26, 27, 28]. The neural network layers utilized in these works are derived from SPDNet [29], which is specifically designed to process SPD matrices while preserving their SPD properties throughout the layers during non-linear learning. These layers include:

• •

  Bi-Map Layer: This layer transforms the covariance matrix $S$ by applying the Bi-Map transformation $WSW^{T}$, where the transformation matrix $W$ must have full row rank.

• •

  ReEig Layer: This layer introduces non-linearity to SPD manifolds using $U\max{(\epsilon I,\Sigma)}U^{T}$, where $S=U\Sigma U^{T}$, $\epsilon$ is a rectification threshold, and $I$ is an identity matrix.

• •

  LOG Layer: This layer maps elements on SPD manifolds to their tangent space using $U\log{(\Sigma)}U^{T}$, where $S=U\Sigma U^{T}$. The logarithm function applied to a diagonal matrix takes the logarithm of each diagonal element individually.

we first calculate the gradient of any squared distance function on general Riemannian manifolds.

We say a function $\psi:\mathcal{M}\longmapsto\mathbb{R}\cup\{-\infty\}$ is c-concave if it is not identically $-\infty$ and there exists $\varphi:\mathcal{M}\longmapsto\mathbb{R}\cup\{\pm\infty\}$ such that

where $c(p,q)$ is a nonnegative cost function that measures the cost of transporting mass from $q$ to $p\in\mathcal{M}$. Brenier’s polar factorization on general Riemannian manifolds is given as follows:

On the SPD manifold ($\mathcal{S}_{++}$, $g^{LEM}$), let’s consider source measure $\mu$ and target measure $\nu\in\mathbb{P}(\mathcal{S}_{++})$ that are accessible only through finite samples $\{S_{i}\}_{i=1}^{N}\sim\mu$ and $\{\bar{S}_{j}\}_{j=1}^{M}\sim\nu$. Each sample contributes to its respective empirical distribution with weights $\frac{1}{N}$ and $\frac{1}{M}$. We assume there exists a linear transformation on ($\mathcal{S}_{++}$, $g^{LEM}$) between two discrete feature space distributions of the source and target domains that can be estimated with the discrete Kantorovich problem, i.e., there exists

where $\prod(\mu,\nu):=\{\gamma\in\mathbb{R}_{\geq 0}^{N\times M}\big{|}\gamma\mathbbm{1}_{N}=\frac{1}{N},\text{ and }\gamma^{\top}\mathbbm{1}_{M}=\frac{1}{M}\}$. The pseudocode for the computation of the transport plan is presented in Algorithm 1.

In the following paragraph, we will employ discrete transport plans to construct c-concave functions and utilize Lemma 1 and 2 to obtain continuous transport plans. According to Algorithm 1, the new coordinates for target samples on SPD manifolds are given by the following expression:

We define $\varphi(\bar{S}):=+\infty$ for $\bar{S}\notin\{\bar{S}^{\dagger}_{j}\}_{j=1}^{M}$, and set the cost function $c(S,\bar{S}):=\frac{1}{2}\cdot d_{g^{LEM}}^{2}(S,\bar{S})$. This leads to the formulation of the corresponding c-concave function as follows:

For a fixed $S$, we suppose $j^{\ast}\in\{1,...,M\}$ achieves the minimum. By invoking Lemma 1 and 2, we promptly derive the unique optimal transport map as follows,

which implies the following theorem:

Theorem 3 encompasses an extension of [6, Theorem 3.1] pertaining to ($\mathcal{S}_{++}$, $g^{LEM}$) within the context of the OT-DA framework. Significantly, our framework is designed to incorporate any supervised neural networks-based transformation, distinguishing it from the conventional use of affine transformations, and employs the squared Riemannian distance. In Appendix VI-A, we illuminate the connection between our proposed approach and another similar OT-DA framework on SPD manifolds proposed in [10] that utilizes the squared $\ell_{2}$ Euclidean distance in its cost function.

We propose a novel geometric deep learning approach, Deep Optimal Transport (DOT) on ($\mathcal{S}_{++}$, $g^{LEM}$), to address the joint distribution adaptation problems involving both marginal distribution shifts and conditional distribution shifts. We assume samples from the source and target distributions can be restored via OT. The proposed architecture is depicted in Figure 1.

We model the space of EEG spatial covariance matrices in ($\mathcal{S}_{++}$, $g^{LEM}$) and employ DOT to address the issue that EEG statistical characterization is shifted during the feedback phase. Using the event-related desynchronization/synchronization methods to classify motor imagery classes, EEG frequency components in specific frequency bands are found which provided the best discrimination between movement imagination [37]. Hence, it is reasonable to regard features in source and target feature space holds the same, i.e., $X_{S}=X_{T}$, after transforming to a latent space using the Bi-Map layer, and thus model the space of EEG spatial covariance matrices using SPD manifolds with the Log-Euclidean metric. Moreover, the Bi-Map layer provides a transitive group action on SPD manifolds, allowing any transition between two points on manifolds. In Fig. 2, we illustrate how the DOT-based MI-EEG classifier functions in a domain adaptation problem.

Specifically, EEGs are represented by $X\in\mathbb{R}^{n_{C}\times n_{T}}$, where $n_{C}$ and $n_{T}$ denote the number of channels and timestamps, respectively. The corresponding EEG spatial covariance matrix, $S:=S(\Delta f\times\Delta t)\in\mathcal{S}_{++}^{n_{C}}$, represents an EEG segment within the frequency band $\Delta f$ and time interval $\Delta t$. Let $\mathcal{F}$ and $\mathcal{T}$ represent the frequency and time domains, respectively. A given segmentation $Seg(\mathcal{F},\mathcal{T})$ is written as $\{\Delta f\times\Delta t|\Delta f\subset Range(\mathcal{F})\text{, and }\Delta t\subset Range(\mathcal{T})\}$, on the time-frequency domain of EEG signals. To ensure the theoretical convergence of this classifier in a real-world context, the following three assumptions are necessary:

Assumption 1 is a necessary condition of modeling the OT-DA problem. Assumption 2 is straightforward because when we are uncertain about which EEG frequency band has a higher weight for a specific task, we tend to assume that each frequency band has equal priority for the classification task. Assumption 3 ensures that features extracted from different bands do not overlap.

Three commonly used EEG-based motor imagery datasets are used for evaluation, accessed from two publicly available data packages MOABB and BNCI Horizon 2020¹¹1  The MOABB package can be accessed at https://github.com/NeuroTechX/moabb, and the BNCI Horizon 2020 project is available at http://bnci-horizon-2020.eu/database/data-sets.. The selection criteria require each subject in the dataset to have two-session recordings collected on different days. All digital signals were filtered using Chebyshev Type II filters at intervals of 4 Hz. The filters were designed to have a maximum loss of 3 dB in the passband and a minimum attenuation of 30 dB in the stopband. The following are the descriptions of these three datasets.

1). Korea University (KU) Dataset: The KU dataset comprises EEG signals obtained from 54 subjects during a binary-class imagery task. The EEG signals were recorded using 62 electrodes at a sampling rate of 1,000 Hz. For evaluation, we selected 20 electrodes located in the motor cortex region. The dataset was divided into two sessions, S1 and S2, each consisting of a training and a testing phase. Each phase included 100 trials, with an equal number of right and left-hand imagery tasks. The EEG signals were epoched from second 1 to second 3.5 relative to the stimulus onset, resulting in a total epoch duration of 2.5 seconds.

2). BNCI2014001 Dataset: The BNCI2014001 dataset (also known as BCIC-IV-2a) involved 9 participants performing a motor imagery task with four classes: left hand, right hand, feet, and tongue. The task was recorded using 22 Ag/AgCl electrodes and three EOG channels at a sampling rate of 250 Hz. The recorded data were filtered using a bandpass filter between 0.5 and 100 Hz, with a 50 Hz notch filter applied. The study was conducted over two sessions on separate days, with each participant completing a total of 288 trials (six runs of 12 cue-based trials for each class). The EEG signals were epoched from the cue onset at second 2 to the end of the motor imagery period at second 6, resulting in a total epoch duration of 4 seconds.

3). BNCI2015001 Dataset: BNCI2015001 involves 12 participants performing a right-hand versus both-feet imagery task. Most of the study consisted of two sessions conducted on consecutive days, with each participant completing 100 trials for each class, resulting in a total of 200 trials per participant. Four subjects participated in a third session, but their data should have been included in the evaluation. The EEG data were acquired at a sampling rate of 512 Hz and filtered using a bandpass filter ranging from 0.5 to 100 Hz, with a notch filter at 50 Hz. The recordings began at second 3 with the prompt and continued until the end of the cross period at second 8, resulting in a total duration of 5 seconds.

The experimental settings for domain adaptation on three evaluated datasets are described as follows, illustrated in Fig. 3: 1). Semi-supervised Domain Adaptation: For the KU dataset, S1 was used for training, and the first half of S2 was used for validation, while the second half was reserved for testing. Early stopping was adopted during the validation phase. 2). Unsupervised Domain Adaptation: For the BNCI2014001 and BNCI2015001 datasets, S1 was used for training, and S2 was used for testing. No validation set was used, and a maximum epoch value was set for stopping in practice.

We specified that the transfer procedure must always occur from S1 (source domain) to S2 (target domain), i.e., $S1\rightarrow S2$ for each dataset. It is reasonable because there is no meaning if negative transfer is used in the learning process. The former and latter sessions could be interpreted as the calibration and feedback phases.

The geometric methods in this study can be broadly divided into four categories. These include categories of parallel transport, deep parallel transport, optimal transport, and deep optimal transport. Each category includes various transfer methods, performing on three types of base machine learning classifiers: the support vector machine, minimum distance to Riemannian mean, and SPDNet. We begin by introducing these base classifiers and discussing transfer methods in each category.

1. 1.

   Support Vector Machine (SVM) [38]: This is a well-known machine learning classifier used for classification and regression tasks, which aims to find an optimal hyperplane that maximally separates data points of different classes, allowing it to handle linear and non-linear classification problems effectively.

2. 2.

   Minimum Distance to Riemannian Mean (MDM) [21]: This employs the geodesic distance on SPD manifolds for classification. More specifically, a centroid is estimated for each predefined class, and subsequently, the classification of a new data point is determined based on its proximity to the nearest centroid.

3. 3.

   SPDNet [29]: This has been introduced in the preliminary part. In the experiments, we employ a one-depth Bi-Map layer with both the input and output dimensions equal to the number of electrodes on the primary motor cortex: 20 (KU), 22 (BNCI2014001), and 13 (BNCI2015001).

As computations on the SPD manifold are related to Riemannian metrics, we will emphasize the specific Riemannian metrics associated with each method. The different transfer methods across categories according to various Riemannian metrics can be seen in the first three columns of Table II.

We synthesize a source domain dataset consisting of 50 covariance matrices of size $2\times 2$, drawing a Gaussian distribution on the SPD manifold proposed in a generative approach [43] with a median of 1 and a standard deviation of 0.4. We utilize the following transformation matrix in the generation,

and data in the target domain can be expressed as $WSW^{\top}$.

Subfigure (a) exhibits red points sourced from the source domain, with the blue points originating from the target domain. In Subfigure (b), the blue points represent the relocation of the red points via the EMD method from OT. The associations made by this transformation are represented by thin blue lines. Remarkably, the position of these transformed blue points is extremely proximate to the points of the target domain in Subfigure (a). This can be attributed to the fact that the OT-DA method essentially computes a weighted average of the source and target domain data points, based on weights learned to maximize the transportation costs.

The transformed blue points in Subfigures (c) and (d) were derived from mapping the red points via neural network weights learned through the DOT and DeepJDOT methods, respectively. Corresponding point pairs are linked with thin blue lines. These figures reveal the shared fundamental characteristic of the DOT and DeepJDOT methods that both project the source and target domains onto a subspace learned via a network, which, in these two figures, is represented by a straight line. In differential geometry, this is a submanifold, which refers to a lower-dimensional manifold embedded within a higher-dimensional one. In this case, the two-dimensional matrices are represented by a line, essentially projecting the higher dimensional data onto a lower dimensional space. It reduces the complexity of the data while preserving information.

As shown in Table II, we classify all approaches into three categories according to Riemannian metrics. The first category does not require a Riemannian metric for computation, the second employs $g^{AIRM}$, while the third uses $g^{LEM}$.

The experimental results indicate that metrics strongly influence the results. Specifically, in the table, EMD+SVM achieved 50.89 ($g^{AIRM}$) and 68.33 ($g^{LEM}$) on the BNCI2014001 dataset. Moreover, several observations are introduced as follows: 1). Participant Number: The improvement is related to the number of participants. On datasets with fewer participants, the enhancement brought by transfer methods is more significant, with several scenarios showing changes of more than 5%. 2). Base Classifier Results: The extent of improvement is linked to the base classifier’s performance. Lower initial results from the base classifier provide more opportunities for transfer methods to achieve significant gains. 3). Methodological Factors: Improvements are associated with the classifier neural network architecture, parameter selection, and data preprocessing. Therefore, the conclusions of this study may not be universally applicable.

Furthermore, we summarize the observations of all methods affected by various Riemannian metrics as follows:

• •

  ($\mathcal{S}_{++}$, $g^{AIRM}$) includes categories of OT, PT, and DPT. Combining the results from all three datasets, the RCT method performs well, especially when paired with the MDM or SPDNet classifier. The EMD method shows the largest variation across datasets. The second step (ROT) in the RPA method, which is limited to semi-supervised scenarios, performs worse than the first step (RCT) of RPA. The RieBN method provides slight improvements but does not outperform RCT across the board.

• •

  ($\mathcal{S}_{++}$, $g^{LEM}$) includes two categories of OT and DOT. The EMD method, when paired with the SVM classifier, performs significantly better compared to its pairing with the SPDNet classifier. SPDSW and logSW methods have obvious improvements compared to their base classifier (SVM), among which SPDSW performs better than logSW by 1% to 2% on all three datasets. The MDA and CDA methods perform similarly, except for the BNCI2014001 dataset, the improvements on the other two datasets are between 1% to 2%, while the MDA method is slightly better than the CDA method or a mixture of the MDA and CDA methods. Notably, the use of pseudo-labels in the CDA method affects its results. The DeepJDOT method’s performance across three datasets falls short of our proposed approach. The primary reason for this could be attributed to our decision to replace the task of computing transportation costs for each data point with merely computing transportation costs for the centers of each frequency band’s dataset. Minimizing the transportation costs does not necessarily imply improved classification accuracy, thus the proposed approach reduces computational cost while enhancing the model’s potential.

Besides, neural network-based approaches generally achieve better classification outcomes compared to other categories of approaches, such as EMD, SPDSW, and logSW. To compute the CDA loss in DOT, we used the pseudo-label method, with pseudo labels derived from the Graph-CSPNet [26], which currently performs the best for subject-specific tasks. Specifically, the results of SPDSW and logSW on the BNCI2014001 dataset were obtained using the data processing procedure and source codes provided by their GitHub repository, outperforming ours. This is because the performance of the SVM methods heavily depends on a selection of handcrafted features.

We claim that Theorem 3 with cost function $c(x,\bar{x})=||x-\bar{x}||_{\ell^{2}}^{2}$ is equilvant to [6, Theorem 3.1]. It was first observed and proved by Yair et al.[10] who use the matrix vectorization operator to establish an isometry between the affine transformation and the Bi-Map transformation. We follow their arguments and revise a bit of their proof.

To expose the relationship between the Bi-Map transformation in Theorem 3 and affine transformation in Theorem 4, we introduce the matrix vectorization operator $vec(\cdot)$ to stack the column vectors of matrix $A=(a_{1}\big{|}a_{2}\big{|}\cdots\big{|}a_{n})$ below on another as follows,

Then, we have the following lemma.

According to Lemma 5 and 6, the Bi-Map transformation $W\cdot S\cdot W^{\top}$ turns to be affine transformation $A=W\otimes W$, and thus Theorem 3. Hence, Theorem 3 endowed with the cost function $c\big{(}vec(x),vec(y)\big{)}=||vec(x)-vec(y)||_{\ell^{2}}^{2}$ is equilvant to Theorem 4.

Parallel transport on Riemannian manifolds transfers a tangent vector at one point along a curve on the manifold to another point, such that the transported vector remains parallel to the original vector [30]. Under two Riemannian metrics on SPD manifolds $g^{AIRM}$ and $g^{LEM}$, expressions of parallel transport exhibit distinctly and determine two distinct modeling approaches for the problems [44, 10]. Specifically, for any $S_{1}$ and $S_{2}$ on ($\mathcal{S}_{++}$, $g^{AIRM}$), parallel transport from $S_{1}$ to $S_{2}$ for $s\in T_{S_{1}}\mathcal{S}_{++}$ can be expressed as $\Gamma_{S_{1}\rightarrow S_{2}}(s)=EsE^{\top}$, where $E=(S_{1}\cdot S_{2}^{-1})^{\frac{1}{2}}$. For any $S_{1}$ and $S_{2}$ on ($\mathcal{S}_{++}$, $g^{LEM}$), parallel transport is given by $\Gamma_{S_{1}\rightarrow S_{2}}(s)=s$.

Keep in mind that DOT seeks a neural networks-based transformation $\varphi:X_{S},X_{T}\longmapsto Z\subset\mathcal{S}_{++}$ such that $\mathbb{P}_{S}(\varphi(X_{S}))=\mathbb{P}_{T}(\varphi(X_{T}))$ and $\mathbb{Q}_{S}(Y_{S}|\varphi(X_{S}))=\mathbb{Q}_{T}(Y_{T}|\varphi(X_{T}))$, where $\varphi(X_{S})$ and $\varphi(X_{T})$ are transformation $\varphi$ on the source and target feature spaces. Neural networks-based transformation $\varphi$ extracts discriminative information from feature spaces. After the extraction, it is reasonable to view the feature space of the source domain as the same as it of the target domain, as they will be transformed under neural networks with the same weights, so are their transformed feature spaces, i.e., $X_{T}=X_{S}$ and $\varphi(X_{T})=\varphi(X_{S})\subset Z$, and thus we utilize the Log-Euclidean metric to formulate the space of covariance matrice.