Geodesic Optimization for Predictive Shift Adaptation on EEG data

Authors: apolline  mellot
antoine  collas
sylvain  chevallier
alexandre  gramfort
denis a engemann
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 13:30:15

Abstract: Electroencephalography (EEG) data is often collected from diverse contexts involving different populations and EEG devices. This variability can induce distribution shifts in the data $X$ and in the biomedical variables of interest $y$, thus limiting the application of supervised machine learning (ML) algorithms. While domain adaptation (DA) methods have been developed to mitigate the impact of these shifts, such methods struggle when distribution shifts occur simultaneously in $X$ and $y$. As state-of-the-art ML models for EEG represent the data by spatial covariance matrices, which lie on the Riemannian manifold of Symmetric Positive Definite (SPD) matrices, it is appealing to study DA techniques operating on the SPD manifold. This paper proposes a novel method termed Geodesic Optimization for Predictive Shift Adaptation (GOPSA) to address test-time multi-source DA for situations in which source domains have distinct $y$ distributions. GOPSA exploits the geodesic structure of the Riemannian manifold to jointly learn a domain-specific re-centering operator representing site-specific intercepts and the regression model. We performed empirical benchmarks on the cross-site generalization of age-prediction models with resting-state EEG data from a large multi-national dataset (HarMNqEEG), which included $14$ recording sites and more than $1500$ human participants. Compared to state-of-the-art methods, our results showed that GOPSA achieved significantly higher performance on three regression metrics ($R^2$, MAE, and Spearman's $\rho$) for several source-target site combinations, highlighting its effectiveness in tackling multi-source DA with predictive shifts in EEG data analysis. Our method has the potential to combine the advantages of mixed-effects modeling with machine learning for biomedical applications of EEG, such as multicenter clinical trials.

Main paper: 
 introduction: 
 Machine learning (ML) has enabled advances in the analysis of complex biological signals, such as magneto-and electroencephalography (M/EEG), in diverse applications including biomarker Figure 1 : Joint shift in X and y distributions on the HarMNqEEG dataset [ 33 ] . Subset of mean PSDs (A) and age distributions (B) from three recording sites used for the empirical benchmarks.exploration [ 58 , 20 , 60 ] or developing Brain-Computer Interface (BCI) [ 57 , 15 , 1 , 2 ] . However, a major challenge in applying ML to these signals arises from their inherent variability, a problem commonly referred to as dataset shift [ 12 ] . In the case of M/EEG data this variability can be caused by differences in recording devices (electrode positions and amplifier configurations), recording protocols, population demographics, and inter-subject variability [ 36 , 14 , 23 , 29 ] . Notably, shifts can occur not only in the data X but also in the biomedical variable y we aim to predict, further complicating the use of ML algorithms.Riemannian geometry has significantly advanced EEG data analysis by enabling the use of spatial covariance matrices as EEG descriptors [ 5 , 4 , 6 , 38 , 31 , 35 , 33 , 48 , 17 , 56 ] . In [ 4 , 6 ] , the authors introduced a classification framework for BCI based on the Riemannian geometry of covariance matrices. These methods classify EEG signals directly on the tangent space using the Riemannian manifold of symmetric positive definite (SPD) matrices (S ++ d ), effectively capturing spatial information. More recently, [ 47 , 48 , 7 ] extended this framework to regression problems from M/EEG data in the context of biomarker exploration. Furthermore, [ 47 , 48 ] proved that Riemannian metrics lead to regression models with statistical guarantees in line with log-linear brain dynamics [ 10 ] and are, therefore, well-suited for neuroscience applications. Across various biomarker-exploration tasks and datasets, recent work has shown that Riemannian M/EEG representations offer parameter-sparse alternatives to non-Riemannian deep learning architectures [ 13 , 40 , 17 ] . Domain adaptation addresses the challenges posed by differences in data distributions between source and target domains, e.g., when data are recorded with different cameras in computer vision [ 55 ] , different writing styles in natural language processing [ 32 ] , or varying sensor setups in time series analysis [ 22 ] . In particular, DA considers target shift where the shift is in the outcome variable y. For classification it means source and target data share the same labels but in different proportions [ 34 ] . Target shift is also frequent in the context of multicenter neuroscience studies, as the studied population of one site may vary significantly from the studied population of another site (cf. Figure 1 ). To tackle various sources of variability in neurophysiological data like EEG, there is a need for a DA approach that can deal with a joint shift in X and y.In [ 62 ] , the authors addressed DA for EEG-based BCI using re-centering affine transformation of covariance matrices ( Section 2 ) to align data from different sessions or human participants, improving classification accuracies. Yair et al. [ 59 ] extended this with parallel transport showing its effectiveness in EEG analysis, whereas, Peng et al [ 43 ] introduced a domain-specific regularizer based on the Riemannian mean. Notably, this parallel transport approach reduces to [ 62 ] when the common reference point is the identity. In a deep learning context, Kobler et al. [ 31 ] proposed to do a per-domain online re-centering which can be seen as a domain specific Riemannian batch norm. Going beyond re-centering, Riemannian Procrustes Analysis (RPA) [ 45 ] was proposed for EEG transfer learning, using three steps: mean alignment, dispersion matching, and rotation correction. However, the rotation step is unsuitable for regression problems and RPA adapts only a single source to a target domain. Recently, [ 36 ] demonstrated the benefits of re-centering for regression problems, showing improvements in handling task variations in MEG and enhancing across-dataset inference in EEG.On the other hand, mixed-effects models (or multilevel models) have been successfully used to tackle data shifts in X and y [ 16 , 25 ] . In biomedical data, mixed-effects models are crucial due to the presence of common effects, such as disease status and age. Riemannian mixed-effects models have been used to analyze observations on Riemannian manifolds, accommodating individual trajectories with mixed effects at both group and individual levels [ 30 , 49 , 50 ] . These models adapt a base point on the manifold for each data point and utilize parallel transport for this adaptation, which is necessary for accurate trajectory modeling. However, they differ significantly from the problem we address in this work. Notably, the input data X are covariates (e.g., age or disease status) which belong to a Euclidean space and the variables y to predict belong to the manifold (e.g., MRI diffusion tensors on S ++ d ) which is the opposite of the paper's studied problem. This distinction is critical as it highlights that while both methods use the geometry of Riemannian manifolds, the nature of the predicted variables and the type of data used differ from existing Riemannian mixed-effects models.In this work, we address the challenging problem of multi-source domain adaptation with predictive shifts on the SPD manifold, focusing on distribution shifts in both the input data X and the variable to predict y. We propose a novel method called Geodesic Optimization for Predictive Shift Adaptation (GOPSA). It enables mixed-effects modeling by jointly learning parallel transport along a geodesic for each domain and a global regression model common to all domains, with the assumption that the mean ȳT of the target domain is known. GOPSA aims to advance the state of the art by: (i) addressing shifts in both covariance matrices and the outcome variable y, (ii) being tailored for regression problems, and (iii) being a multi-source test-time domain adaptation method, meaning that once trained on source data, it can generalize to any target domain without requiring access to source data or retraining a new model.We first introduce in Section 2 how to do regression from covariance matrices on the Riemannian manifold of S ++ d . We also interpret classical learning methods on S ++ d from heterogeneous domains as parallel transports combined with Riemannian logarithmic mappings. This leads us to GOPSA in Section 3 , which learns to parallel transport each domain, with algorithms at train and test times. Finally, in Section 4 , we apply GOPSA as well as different baselines on simulated data and the HarMNqEEG dataset.Notations Vectors and matrices are represented by small and large cap boldface letters respectively (e.g., x, X). The set {1, ..., K} is denoted by 1, K . S ++ is a smooth manifold [ 9 ] . A vector space is defined at each Σ ∈ S ++ d , called the tangent space, denoted T Σ S ++ d , and is equal to S d , the ambient space. Equipped with a smooth inner product at every tangent space, a smooth manifold b ecomes a Riemannian manifold. To do so, we make use of the affine invariant Riemannian metric [ 51 , 44 ] . Given Γ, Γ ′ ∈ T Σ S ++ d , this metric is ⟨Γ, Γ ′ ⟩ Σ = tr Σ -1 ΓΣ -1 Γ ′ .Riemannian mean The Riemannian distance (or geodesic distance) associated with the affineinvariant metric is δ R (Σ, Σ ′ ) ≜ log Σ -1 /2 Σ ′ Σ -1 /2F with log : S ++ d → S d being the matrix logarithm (see Appendix A.1 for a definition). This distance is used to compute the Riemannian mean Σ defined for a set{Σ i } N i=1 ⊂ S ++ d as Σ ≜ arg min Σ∈S ++ d N i=1 δ R (Σ, Σ i ) 2 .This mean is efficiently computed with a Riemannian gradient descent [ 44 , 63 ] .The idea of the covariance-based approach is to define nonlinear feature transformations into vectors that can be used as input for classical linear machine learning models. To do so, the Riemannian logarithmic mapping of Σ ′ at Σ is defined asEQUATIONThus, matrices in the Riemannian manifold S ++ d are transformed into tangent vectors.Parallel transport A classical practice to align distributions is parallel transport of covariance matrices from their mean to the identity and then apply the logarithmic mapping ( 1 ) . Parallel transport along a curve allows to move SPD matrices from one point on the curve to another point on the curve while keeping the inner product between the logarithmic mappings with any other vector transported along the same curve constant. The following lemma gives the parallel transport of Σ ′ from Σ to I d along the geodesic between these two points (See proof in Appendix A.2 ). Lemma 2.1 (Parallel transport to the identity). Given Σ, Σ ′ ∈ S ++ d , the parallel transport of Σ ′ along the geodesic from Σ to the identityI d at α ∈ [0, 1] is PT (Σ ′ , Σ, α) ≜ Σ -α /2 Σ ′ Σ -α /2 .Learning on S ++ dThe logarithmic mapping ( 1 ) at the identity is simply the matrix logarithm. Thus, a classical non-linear feature extraction [ 6 , 36 , 8 ] of a dataset {Σ i } N i=1 of Riemannian mean Σ combines parallel transport and logarithmic mapping at the identity,EQUATIONwhere uvec vectorizes the upper triangular part with off-diagonal elements multiplied by √ 2 to preserve the norm. Correcting dataset shifts by re-centering all source datasets [ 62 ] , corresponds to parallel transporting dataEQUATION)This method is the go-to approach for reducing shifts of the covariance matrix distributions coming from different domains and has been applied successfully for brain-computer interfaces [ 45 , 59 ] and age prediction from M/EEG data [ 36 ] .3 Learning to recenter from highly shifted y distributions with GOPSAIn this section, we develop a novel multi-source domain adaptation method, called Geodesic Optimization for Predictive Shift Adaptation (GOPSA), that operates on the S ++ d manifold and is capable of handling vastly different distributions of y. Our approach implements a Riemannian mixed-effects model, which consists of two components: a single parameter estimating a geodesic intercept specific to each domain and a set of parameters shared across domains.At train-time, GOPSA jointly learns the parallel transport of each of the K source domains and the regression model shared across domains. At test-time, we assume having access to the target mean response value y T and predict on the unlabeled target domain of covariance matrices (Σ T ,i ) N T i=1 . GOPSA focuses solely on learning the parallel transport of the target domain so that the mean prediction, using the regression model learned at train-time, matches y T .Parallel transport along the geodesic In Section 2 , we presented how domain adaptation is performed on S ++ d . In particular, ( 3 ) presents how to account for data shifts of each domain. However, this operator can only work if the variability between domains is considered as noise. As explained earlier, we are interested in shifts in both features and the response variable. Thus, ( 3 ) discards shift coming from the response variable and hence harms the performance of the predictive model. Based on the Lemma 2.1 , we propose to parallel transport features to any point on the geodesic between a domain-specific Riemannian mean Σ k and the identity. Indeed, GOPSA parallel transports Σ k,i on this geodesic with α ∈ [0, 1] and then applies the Riemannian logarithmic mapping ( 1 ) at the identity,EQUATIONThis allows each domain to undergo parallel transport to a certain degree, effectively moving it toward the identity.Algorithm 1: Train-Time GOPSA Input: For all k ∈ 1, K , {(Σ k,i , y k,i )} N k i=1 , initialization of γ S , step-sizes {ξ t } t≥1 for k = 1 → K do Σ k ← Riemannian mean of {Σ k,i } N k i=1 end t ←γ T ← γ T -ξ t L ′ T (γ T ) t ← t + 1 end y T ← Z T (γ ⋆ T )β ⋆ S (γ ⋆ S ) return y T 
limitations: 
 Question: Does the paper discuss the limitations of the work performed by the authors? 
train-time: 
 GOPSA aims to learn simultaneously features from ( 4 ) and a regression model. To do so, we solve the following optimization problemEQUATIONwith α S = [α 1 , . . . , α K ] ⊤ . This cost function is decomposed into three key aspects. First, covariance matrices undergo parallel transported using Lemma 2.1 to account for shifts between domains. Second, they are vectorized, and a linear regression predicts the output variable from these vectorized features. Third, the coefficients of the linear regression β S and the α S are learned jointly so that the predictor is adapted to the parallel transport and reciprocally. Besides, to enforce the constraint on α S , we re-parameterize it using the sigmoid function, which defines a bijection between R and (0, 1), thereby ensuring that the resulting α S values lie within the desired range:α k = σ(γ k ) ≜ (1 + exp(-γ k )) -1. Thus, the constrained problem ( 5 ) can be formulated as the following unconstrained optimization problemEQUATIONwith γ S = [γ 1 , . . . , γ K ] ⊤ . Let us define the matrix Z S (γ) ∈ R N S × d(d+1) /2 , with N S = K k=1 N k , as the concatenation of the source data, where each row corresponds to a feature vector:EQUATIONIn the same manner, the source labels are concatenated to y S = [y 1,1 , . . . , y K,N K ] ⊤ ∈ R N S . Given a fixed γ S , the problem ( 6) is solved with the ordinary least squares estimator. In practice, we choose to regularize the estimation of the linear regression with a Ridge penalty. Thus, ( 6) is rewritten asEQUATIONwhere β ⋆ S (γ) ∈ R d(d+1) /2 are the Ridge estimated coefficients given a fixed γ and λ > 0 is the regularization hyperparameter. The problem ( 8) is efficiently solved with any gradient-based solver [ 39 ] .The train-time of GOPSA is summarized in Algorithm 1 . The proposed training algorithm begins by calculating the Riemannian mean of covariance matrices for each domain k. It then iteratively optimizes the parameters γ S by computing the feature matrix ( 7 ) , determining Ridge regression coefficients ( 8 ) , and updating γ S using a gradient descent step on the loss function ( 8 ) until convergence. The output result is the optimized Ridge regression coefficients. For clarity of presentation, Algorithm 1 employs a gradient descent. In practice, we use L-BFGS and obtain the gradient using automatic differentiation through the Ridge solution that is plugged into the loss in ( 8 ).At test-time, we now have a fitted linear model on source data with coefficients β ⋆ S (γ ⋆ S ). The goal is to adapt a new target domain (Σ T ,i ) N T i=1 for which the average outcome y T is assumed to be known. First, let us define the matrix Z T (γ) ∈ R N T × d(d+1) /2 as the concatenation of the target dataEQUATIONThen, GOPSA adapts to this new target domain by minimizing the error between y T and its estimation computed with the fitted linear model. This minimization is performed with respect to γ T ∈ R that parametrizes the parallel transport of the target domain, i.e.,EQUATION)Finally, the predictions on the target domain areEQUATIONThe test-time procedure of GOPSA is summarized in Algorithm 2 . The algorithm begins by calculating the Riemannian mean of the target covariance matrices {Σ T ,i } N T i=1 . It then iteratively optimizes the parameter γ T by computing the feature matrix ( 9), the derivative of the loss function ( 10 ) , and updating γ T using a gradient descent step until convergence. The algorithm determines the estimated target outcomes, y T , by using the optimized γ ⋆ T on the feature matrix, combined with the pre-trained regression coefficients β ⋆ S (γ ⋆ S ). The output result is the predicted target outcomes y T . It should be noted that, once again, for clarity of presentation, Algorithm 2 employs a gradient descent, but other derivative-based optimization methods can be used.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof?Answer: [Yes] Justification: The proposed method uses the parallel transport from a given SPD matrix to the identity. The formula is presented in Lemma 2.1 , and the proof is provided in Appendix A.2 .Guidelines:• The answer NA means that the paper does not include theoretical results.• All the theorems, formulas, and proofs in the paper should be numbered and crossreferenced. • All assumptions should be clearly stated or referenced in the statement of any theorems.• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 
empirical benchmarks: 
 In this section, we built empirical benchmark to evaluate the performance of GOPSA. We first present the simulated data that we used to illustrate the relevance of our method when there is a joint distribution shift of the data and the labels. Then, we present the EEG dataset that we used to evaluate the performance of GOPSA with real data from different recording sites. Finally, we present the baseline methods that are compared with GOPSA.Simulated data To generate simulated data, we used the generative model described in [ 47 , 48 , 36 ] . The data follow the classical instantaneous mixing model:EQUATIONwhere x i (t) ∈ R d are the observed time-series, s i (t) ∈ R d are the underlying signal of the neural generators and A is the mixing matrix whose columns are the observed spatial patterns of the neural generators. Furthermore, we assume that y follow a log-linear model:EQUATIONwhere p ℓi > 0 is the variance of the ℓ-th element of the underlying signal s i (t) as introduced in [ 47 , 48 , 36 ] . From this, we generate domains (source and target) by applying shifts on X and y. To do so, we introduced a per-domain shift in the data distribution by applying an affine transformation to the covariance matricesEQUATION)with B k ∈ S ++ d and ξ ≥ 0 controlling the amplitude of the shift. Then, we shifted the label distribution by modifying the variance of the underlying signal p ℓi : (15) with ξ ≥ 0 still controlling the amplitude of the shift. Thus, the distribution of y is shifted per domain because of the log-linear relationship of ( 13 ) . It should be noted that β is kept constant across domains.p ℓi → p 1+kξ ℓiHarMNqEEG dataset The HarMNqEEG dataset [ 33 ] was used for our numerical experiments. This dataset includes EEG recordings collected from 1564 participants across 14 different study sites, distributed across 9 countries. In our analysis, we consider each study site as a distinct domain. Appendix A.5 provides detailed demographic information. The EEG data were recorded with the same montage of 19 channels of the 10/20 International Electrodes Positioning System. The dataset provides pre-computed cross-spectral tensors for each participant rather than raw data, and anonymized metadata including the age and the sex of the participants. More precisely, the shared data consists of cross-spectral matrices with a frequency range of 1.17 Hz to 19.14 Hz, sampled at a resolution of 0.39 Hz. A standardized recording protocol was enforced to ensure the consistency across EEG recording of the dataset. In addition to recording constraints, this protocol included artifact cleaning procedures. The cross-spectrum were computed using Bartlett's method (See Appendix A.3 ). Our pre-processing steps were guided by the pre-processing pipeline outlined in [ 33 ] . First, we performed a common average reference (CAR) on all cross-spectrum (See Appendix A.3 ) as different EEG references were used across domains. Subsequently, we extracted the real part of the cross-spectral tensor to obtain co-spectrum tensors containing frequency-specific covariance estimates along the frequency spectrum. Due to the linear dependence between channels introduced by the CAR, the covariance matrices are rank deficient. To address this, we applied a shrinkage regularization with a coefficient of 10 -5 to the data. Additionally, we implemented a global-scale factor (GSF) correction, which compensates for amplitude variations between EEG recordings by scaling the covariance matrices with a subject-specific factor [ 24 , 33 ] (See Appendix A.3 ). Following these pre-processing steps, we obtained a set of 49 covariance matrices for each EEG recording, with each matrix corresponding to a specific frequency bin of the EEG signal. This pre-processed co-spectrum served as the input data for our domain adaptation study.To evaluate the performance of the compared methods, we conducted experiments across several combinations of source and target sites. We selected source domains such that the union distribution of their predictive variable y encompasses a broad age range. All remaining sites were assigned as target domains. For each source-target combination we performed a stratified shuffle split approach with 100 repetitions on the target data. Stratification was based on the recording sites to ensure that each split contained a balanced proportion of participants from each site. The regularization parameter λ in Ridge regression was selected with a nested cross-validation (grid search) over a logarithmic grid of values from 10 -1 to 10 5 . To evaluate the benefit of GOPSA, we compared it against four baselines. Detailed mathematical formulations of these baselines can be found in Appendix A.4 . For each baseline method, the regression task was performed with Ridge regression.Domain-aware dummy model (DO Dummy) As GOPSA requires access to the mean ȳk of each domain, we used a domain-aware dummy model predicting always the mean ȳk of each domain.No re-center / No domain adaptation (No DA) This second baseline method involves applying the regression pipeline outlined in [ 47 , 48 ] without any re-centering. In this setup, all covariance matrices are projected to the tangent space at the source geometric mean Σ computed from all source points, no matter their recording sites.Re-center to a common reference point (Re-center and Re-scale) As introduced in Section 2 , a common transfer learning approach is a Riemannian re-centering of all domains to a common point on the manifold [ 62 , 36 ] . This baseline thus correspond to re-centering each domain k, source and target, independently by whitening them by their respective geometric mean Σ k . An extension of this approach is to perform a Riemannian re-scaling of all domains to a common dispersion, as presented in [ 45 , 36 ] . Domain-aware intercept (DO Intercept) This method consists in fitting one intercept β 0 per domain. In practice since we assume to know ȳT , we correct the predicted values so that their mean is equal to ȳT . This approach is in line with defining mixed-effects models on the Riemannian manifold [ 33 ] .Deep learning (GREEN) The GREEN model [ 40 ] source domains and 1 target domain, with shifts controlled by ξ (0 to maximum). Data are generated 100 times, with 5 sensors and 300 covariance matrices per domain. The target domain is randomly selected between the 6 domains generated as presented in Section 4 , with the remaining domains used as sources. (A) A shift is applied on the covariance matrices following ( 14). (B) A shift is applied on the variances following ( 15 ) . (C) Both shifts from ( 14) and ( 15) are applied simultaneously.used the 'G2' variant of GREEN, which starts at the covariance matrices level and includes pooling layers. This variant is designed for SPD matrices, making it an SPD network [ 26 ] . Although GREEN has been evaluated on multiple datasets for various predictive tasks, it has not yet been applied in a domain adaptation context and does not include an adaptation layer.We applied the domain-adaptation methods independently to each of the 49 frequency bins, resulting in 49 geometric means per domain, except for GREEN, which processes all frequency bands simultaneously. y T of each domain was estimated on target splits (50% of the data) that do not overlap with the evaluation target splits (50% remaining). Statistical inference for model comparisons was implemented with a corrected t-test following [ 37 ] . Experiments with 100 repetitions and all site combinations have been run on a standard Slurm cluster for 12 hours with 250 CPU cores.Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?Answer: [Yes] Justification: The empirical benchmarks ( Section 4 ) use the HarMNqEEG dataset which is available in open access. All the preprocessing steps, as well as the baselines, are extensively detailed in Appendix A.3 and Appendix A.4 , respectively. The technical details (such as gradient computation and optimization algorithm) to implement the proposed method are provided in Section 3 . 
results: 
 Simulated data Figure 2 presents the results of simulated experiments where shifts are applied on either X, y, or both (X, y) as presented in Section 4 . All methods were evaluated in three simulation scenarios: shift in X only, shift in y only, and joint shift in X and y. The intensity of the shift was controlled by ξ in all scenarios. If there is no shift in X, we observe that No DA perfectly estimates the y because the log-linear model is easily estimated across domains even when the y distribution changes ( Figure 2 B ). The performance of No DA however drops when a shift in X is introduced ( Figure 2 A and C). Re-center and Re-scale led to the same results as no scaling shift was applied in the simulation. Both were able to correct the shift in X, but performed poorly when a shift in y was added ( Figure 2 B and C). GREEN notably showed consistant performance across all scenarios, and was relatively resistant to both types of shifts given it is not designed for domain adaptation. DO Intercept and GOPSA showed the best performance across all scenarios, with an advantage for GOPSA. The interest of GOPSA is to estimate this log-linear model with shifts in (X, y) per domain ( Figure 2 C ) which other methods were not able to do. These experiments demonstrate the efficiency of the proposed method in estimating shifts in X between domains even in the presence of a shift in y, contrary to the baseline methods. Theoretically, based on the generative model of the simulated data, the data X and outcome y are linked by a log-linear relationship. This implies that, knowing the shift in X for the target domain, predictions can be made even when y distributions do not overlap between the source and target. Since GOPSA estimates the target shift in X by minimizing (ȳmean(ŷ i )) 2 , it is capable of handling such scenarios.HarMNqEEG data We computed benchmarks for five combinations of source sites and we displayed the results for the three metrics selected for performance evaluation, each colored box representing one method ( Figure 3 ). A min-max normalization was applied to each site combinations separately across methods. We first conducted model comparisons in terms of absolute performance across all baselines (A). No DA, without domain specific re-centering, performed worse than DO Dummy in terms of R 2 score and MAE. Re-center and Re-scale led to lower performances across all metrics, which can be expected as the Riemannian mean is correlated with age in our problem setting Figure 1 . Eventhough its architecture does not include an adaptation layer, GREEN As DO intercept and GOPSA showed overlapping performance distributions, we investigated their paired split-wise (nonrescaled) score differences (B). The site-specific differences of GOPSA scores minus DO Intercept are displayed with their associated p-values. For one site combination (Ba,Be,Cho,Co,Cu90,G,R), DO Intercept yielded higher R 2 scores, and no significant difference was found between the two methods for Ba,Co,G. Similarly, no significant difference was observed on Spearman's ρ results for Cu03,M,R,S. Overall, GOPSA significantly outperformed DO Intercept in five site combinations for MAE, four for Spearman's ρ and three for R 2 score. Detailed results for each source-target combination are presented in Appendix A.6 for Spearman's ρ, R 2 score, and MAE. The bottom rows correspond to the mean performance of each method of all site combinations, and their average standard deviation (see Appendix A.8 for associated boxplots). We expected GOPSA to outperform the baseline methods (e.g. DO Intercept) whenever joint (X, y) shifts occur. In our experimental benchmark, GOPSA significantly outperformed the baseline methods in some site combinations, but not all. This allows us to assume that not all site combinations show joint shifts.Model inspection Next, we investigated the impact of the different re-centering approaches on the data Figure 4 . Power spectrum densities (PSDs) were computed as the mean across sensors of the diagonals of the covariance matrices Riemannian mean for each site combination after No DA, Re-center and GOPSA (A). PSDs for No DA display the initial variability between sites without recentering (cf. Figure 1 ). Re-center resulted in flat PSDs because all data were re-centered to the identity. PSDs produced by GOPSA are flattened and more similar across sites compared to No DA without removing too much information, unlike the un-effective Re-center method (cf. Figure 3 ). The alpha values are inspected as a function of the site mean age (B). Re-center leads to alpha values all equal to one as all sites are re-centered to the identity. For GOPSA, we observed a linear relationship between alpha and the sites' mean age (R 2 = 0.99). This is a direct consequence of the optimization process in GOPSA, which thus can be regarded a geodesic intercept in a mixed-effects model. Overall, GOPSA effectively re-centered sites with younger participants closer to the identity matrix. Re-centering sites around a common point helped reduce the shift in X, while not placing all sites at the exact same reference point helped manage the shift in y, hence preserving the statistical associations between X and y.Question: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [Yes] Justification: The dataset HarMNqEEG [ 33 ] is in open access. We provide the code to reproduce the experiments from the raw data. Guidelines:• The answer NA means that paper does not include experiments requiring code. Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
conclusion: 
 We proposed a novel multi-source domain adaptation approach that adapts shifts in X and y simultaneously by learning jointly a domain specific re-centering operator and the regression model. GOPSA was specifically developed to handle joint shifts in the data distribution and the outcome distribution, as illustrated by the simulations in Figure 2 . GOPSA is a test-time method that does not require to retrain a model when a new domain is presented. GOPSA achieved state-of-the-art performance on the HarMNqEEG [ 33 ] dataset with EEG from 14 recording sites and over 1500 participants. Our benchmarks showed a significant gain in performance for three different metrics in a majority of site combinations compared to baseline methods. GOPSA can thus be used by researchers as a decision rule to infer the presence of joint shifts and, hence, serve as a tool for data exploration and model interpretation. While we focused on shallow regression models, the implementation of GOPSA using PyTorch readily supports its inclusion in more complex Riemannian deep learning models [ 26 , 56 , 11 , 40 , 31 ] . This direction seems promising given our observation that GREEN -a simple deep net combining Riemannian computation with a fully connected layer -already possessed some intrinsic robustness to data shifts. This may point at the capacity of the fully-connected layer to provide additional non-linear transformations that can accommodate the data-generating scenario in which continuous log-linear generators are modified in a discrete manner by site factors. More generally, it emphasizes the potential of complex nonlinear methods for domain adaptation, in line with a recent study on the same dataset reporting positive generalization results using a kernel method [ 28 ] . Furthermore, although this work specifically addresses age prediction, the methodology is applicable to a broader range of regression analyses. While GOPSA necessitates knowledge or estimability of the average y per domain, this requirement aligns with that of mixed-effects models [ 16 , 25 , 61 ] , which are extensively employed in biomedical statistics. By combining mixed-effects modeling with Riemannian geometry for EEG, GOPSA opens up various applications at the interface between machine learning and biostatistics, such as, biomarker exploration in large multicenter clinical trials [ 46 , 52 , 53 ] . 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes] Justification: The proposed approach opens up various applications at the interface between machine learning and biostatistics, such as biomarker exploration in large multicenter clinical trials. References are provided in Section 6 .Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 • The paper should discuss whether and how consent was obtained from people whose asset is used. The human research data used in this study was curated by [ 33 ] and conducted by academic researchers. We do not reproduce their data acquisition protocol here and refer to the original work [ 33 ] . We argue that this is appropriate as our work propose a new statistical method and does not present novel biomedical results. Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained? Answer: [Yes] Justification:The data acquisition was approved by institutional review boards as stated in [ 33 ] : "The studies involving human participants were reviewed and approved by the Ethics Committees of all involved institutions. In all cases, the participants and/or their legal guardians/next of kin provided written informed consent to participate in this study. All data were de-identified, and participants gave permission for their data to be shared as part of the informed consent process.". 
 
 Appendix: Given Σ ∈ S ++ d and its Singular Value Decomposition (SVD) Σ = U diag(λ 1 , . . . , λ d )U ⊤ , the matrix logarithm of Σ isEQUATION)Σ to the power α ∈ R isEQUATION)A.2 Proof of Lemma 2.1First, we recall that the geodesic associated with the affine invariant metric from Σ to Σ ′ isEQUATIONHence, Σ♯ α I d = Σ 1-α . From [ 59 ] , the parallel transport ofEQUATION)Hence, the parallel transport ofEQUATIONwhich concludes the proof.Bartlett estimator From [ 33 ] , the features provided in the HarMNqEEG dataset have been computed using the Bartlett's estimator by averaging more than 20 consecutive and non-overlapping segments. Thus, data consist of cross-spectral matrices with a frequency range of f min = 1.17 Hz to f max = 19.14 Hz, sampled at a resolution of ∆ω = 0.39 Hz. These cross-spectral matrices are denoted S k,i (ω) ∈ H ++dwhere k is the site, i the participant and ω ∈ {f min , f min + ∆ω, . . . , f max }.The cross-spectrum matrices S i (ω) were re-referenced from their original montages with a CAR:EQUATION)whereH ≜ I d -1 d 1 ⊤ d /d.Global Scale Factor (GSF) Co-spectrum matrices were re-scaled with an individual scalar ζ k,i that is calculated as the geometric mean of their power spectrum across sensors and frequencies:EQUATIONwhere N ω ≜ fmax-fmin ∆ω + 1. The GSF correction is then applied to the co-spectrum (the real part of the cross-spectrum) for all frequencies ω:EQUATIONThe Σ k,i (ω) ∈ S ++ d are the features used the Section 4 .The four baseline methods used in this work are detailed in the following. For every methods we have access to K labeled source domains, each including N k covariance matrices and their corresponding variables of interest (Σ k,i , y k,i ) N k i=1 . The method DO Dummy and the mixed-effects model baseline DO Intercept both access the mean value y T of the target domain variable to predict. We remind that as the dataset used in the empirical benchmarks is constituted of several frequency bands, each method is applied to each frequency band independently and then computed vectors are concatenated.Domain-aware dummy model (DO Dummy) The DO Dummy simply returns the mean value y T for every predictions of a given target domain.No re-center / No domain adaptation (No DA) The covariance matrices are used as input of the regression pipeline without any re-centering. First, the geometric mean of the source matrices is computed:EQUATION)Then, source feature vectors are extracted with:EQUATION)Finally, the target feature vectors are extracted from the target data (Σ T ,i ) N T i=1 with:EQUATION)Re-center to a common reference point (Re-center) Domains are re-centered to a common reference point, here we decided to use the identity. First, the geometric mean of each domain is computed:EQUATION)Then, feature vectors are extracted using the specific geometric mean of each domain:EQUATION)Covariance matrices of the target domain (Σ T ,i ) N T i=1 are also re-centered to the identity using their geometric mean :EQUATION)For both No DA and Re-center, the regression task was performed using a Ridge regression, which included an intercept:EQUATION)where z k,i is computed with ( 25 ) or ( 28 ) . The predicted values were computed as:EQUATIONwhere z T ,i is computed with ( 26 ) or ( 30 ) .Domain-aware intercept (DO Intercept) In addition to the K labeled source domains, we assume to have access to the mean of the variable to predict of the target domain ȳT . At train-time, we fit a Ridge regression with a specific intercept for each domainEQUATIONThen, at test-time, we fit a new intercept β 0,T using the target features:EQUATIONThe fitted intercept isEQUATIONand the predictions are A.8 Boxplots of each source-target sites for the three metricsEQUATIONThe following figures represent the performance scores that are displayed in subsection A.6 . We also claim this method is state of the art on the HarMNqEEG dataset. These two claims are asserted in Section 3 and in the Results paragraph of Section 4 .Guidelines:• The answer NA means that the abstract and introduction do not include the claims made in the paper. • The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers. • The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings. • It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.Justification: A limitation of the proposed method is to require the mean outcome value y T of the target set. This limitation is stated in the Contributions paragraph of Section 1 and discussed in Section 6 .Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. , with an open-source dataset or instructions for how to construct the dataset). (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 


Decision: This paper presents a method for tackling domain adaptation challenges in EEG data analysis, specifically addressing shifts in both the feature space (represented by SPD matrices) and outcome variables. The approach is novel and well justified. The experimental evaluation seems not sufficient in the initial reviews, and the authors' rebuttal with additional experiments has successfully addressed this issue. The motivation of the study should be more clearly presented in the final version.

Reviews: review: paper_summary: paper_summary: This study presents an approach to learning robust models under joint shifts in X and y (an important issue in healthcare). Focusing on EEG signals, a covariance-matrix-based learning framework is developed to address this challenge. Empirical results on EEG-specific benchmarks demonstrate its overall efficacy.
strength_weakness: strength_weakness: - Lack of discussion/comments on: 1) other relevant solutions for the same problem setting since joint shifts in X and y are an issue for many other biosignal modalities, 2) expected failure modes and future directions based on this approach.
questions: questions: - Q: Is the ComBat harmonization tool (used in imaging and genetics to correct batch effects) a relevant solution within the same problem setting? "ComBat is a batch adjustment method that removes additive and multiplicative differences between sites due to the use of different scanning devices" - https://cran.r-project.org/web/packages/combat.enigma/index.html.
- Q: What happens when source and target distributions for the response variable don't overlap or only partially? (regarding the choice made in line 238).
- Q: What happens if we only have a bad/noisy estimate of the target domain's mean value (in practice)? Is there any intuition on how things would break?
limitations: limitations: See weaknesses.
paper_summary: paper_summary: The authors proposed GOPSA a new approach for alignment of EEG datasets from multiple subjects and sites. The method respects  the Riemannian manifold of the covariance matrices and learns the parallel transport length parameter simultaneously with the regression model used for solving the downstream task. While simple  and intuitive the approach is well justified, rooted in the previous research and  adequately evaluated using a large public EEG dataset and the associated age prediction task.
strength_weakness: strength_weakness: 1. No performance analysis on a simulated data
2. No interpretation of the obtained decision rule.
3. Table 1 - not all the difference in performance of GOPSA with the other competing techniques seem significant.
questions: questions: 1. Could the authors provide some interpretation of the age prediction model? What frequency bands appear pivotal? Which electrodes?
2. Could the authors add the statistical test to Table 1?   
3. Could the authors discuss and compare their approach against the technique used in several DL studies applied to the multisubject datasets when a specific subject adaptation trainable layer is used to interface  his or her data with the "oracle" classification engine?
4. The authors use covariance matrices as features, however it seems that the transport learnt can be applied to the actual multichannel data directly - could the authors discuss this possibility in the light of forming large aligned datasets in the channel x time form rather than their covariance matrices?
limitations: limitations: The authors adequately addressed the limitations of their study, see also Weaknesses section of this review.
paper_summary: paper_summary: This paper presents a method for tackling domain adaptation challenges in EEG data analysis, specifically addressing shifts in both the feature space (represented by SPD matrices) and outcome variables $y$. The proposed method is designed on top of Riemannian mixed-effects model and is tailored for regression problems. Also, a key highlight of the proposed method is its capability to generalize from the source domain to any target domain without the need for retraining.
strength_weakness: strength_weakness: * Motivation Clarity. 
My main concern is that the necessity of investigating domain adaptation (DA) methods on the SPD manifold is not convincingly presented. While SPD manifold representations are prevalent in EEG analysis, I did not see what specific issues were caused by domain shift. Or what unique issues does DA on the SPD manifold present that the proposed method wants to investigate and resolve? A clear problem setup is helpful for understanding the position of the research.

* The proposed method is on top of a Riemannian mixed-effects model (lines 142-143), which has been used to tackle data shifts in $X$ and $y$ as introduced in Related work (lines 68-74). While the authors mention that there is an opposing setting in previous studies (lines 75-78), it is unclear which specific limitation this paper focuses on and what the challenge is. 

* The main characteristics/advantages of the proposed method are less explained. The paper overall explains very well about the derivation of the algorithm, but it is hard to understand the technical significance compared to existing methods and in what aspects it is better. While the work offers a thoughtful generalization of Riemannian mixed-effects model, it seems like it does not include any concepts that are wholly new. 

* It would be helpful to elaborate on why domain shifts occur in both the data and biomedical variables. A minor question is: Is this an EEG-specific issue or common in biomedical data?
questions: questions: Please refer to the weakness section.
limitations: limitations: I understand this is a professional paper on EEG study. However, I am currently unclear about its overall significance, particularly regarding the motivation or research question for investigating DA methods on the SPD manifold. However, my final decision is open to change pending the authors' rebuttal and further discussions.
paper_summary: paper_summary: This paper proposes a novel method, Geodesic Optimization for Predictive Shift Adaptation (GOPSA), for predictive regression modeling with multi-source domain adaptation. The proposed method employs a domain-specific re-centering operator and a regression model using EEG data for age prediction. The method is innovative.
strength_weakness: strength_weakness: 1. This method assumes knowledge of the average value of the target label $\bar y_T$, which is problematic due to the potential information leakage from the test set. As stated in lines 146-149, the method relies solely on $\bar y_T$ for adjusting regression model predictions, making it difficult to evaluate the generalizability of the method on unseen data. 

2. The paper does not compare the proposed method with other alignment techniques such as re-scale and rotation correction [1].

3. There is no comparison with other Riemannian-based methods or state-of-the-art deep learning techniques for age prediction using M/EEG.

4. It may be inappropriate to claim 'achieved state-of-the-art performance' on a dataset without comparing it with other existing methods.



[1] Mellot, Apolline, et al. "Harmonizing and aligning M/EEG datasets with covariance-based techniques to enhance predictive regression modeling." Imaging Neuroscience 1 (2023): 1-23.
questions: questions: The performance differences between the proposed GOPSA and the baseline, domain-aware (DO) intercept appear to be not large enough in several source-target site combinations. The T-test seems correct (100 repetitions per site combination). However, did you check the assumptions of the T-test: (1) whether the differences follow a normal distribution and (2) whether the variances are equal?
limitations: limitations: Only one dataset was used in this study. Experiments on different datasets should be conducted to validate the generalizability of the method.

Comments: comment: I would like to thank the authors for all the additional experiments and the effort they put into addressing most of my concerns, as well as those of the other reviewers. Although there was not enough time to include an additional dataset to further test the generalizability, I am willing to increase the rating from 5 to 7.
Thank you to the authors for the response and efforts. The additional experiments improve the quality of the work and address my concerns. I am increasing my rating to 7 and hope that the motivation behind the study will be more clearly presented in the final version.
Thank you authors for the rebuttal and taking the effort to run additional experiments. I think the points in this rebuttal and others would add value for readers in the discussion and future work sections of the final version. I maintain my initial recommendation.