Apolline Mellot

Antoine Collas

Pedro Luiz

Coelho Rodrigues

Denis Engemann

Alexandre Gramfort

Pedro L C Rodrigues

Harmonizing and aligning M/EEG datasets with covariance-based techniques to enhance predictive regression modeling

Keywords: MEG/EEG, machine learning, dataset shift, domain adaptation, Riemannian geometry, brain age

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Magneto-and electroencephalography (M/EEG) are brain recording methods with a high temporal resolution on the order of milliseconds, offering a unique and noninvasive neuroscience method enabling basic research and * Correspondence: apolline.mellot@inria.fr, 1 Rue Honoré d'Estienne d'Orves, 91120 Palaiseau clinical applications [START_REF] Hari | MEG-EEG Primer[END_REF]. While quantitative approaches to analyzing M/EEG signals have historically focused on detecting statistical effects, the field has progressively embraced machine learning (ML) approaches whose success is evaluated through predictive modeling. In the context of brain health, classification models are widely used for various applications, e.g., for epileptic seizure detection [START_REF] Tzallas | Epileptic seizure detection in EEGs using time-frequency analysis[END_REF], Brain

Computer Interface (BCI) [START_REF] Lotte | A review of classification algorithms for EEG-based brain-computer interfaces: a 10 year update[END_REF], or automatic sleep staging [START_REF] Chambon | A deep learning architecture for temporal sleep stage classification using multivariate and multimodal time series[END_REF][START_REF] Perslev | U-sleep: resilient high-frequency sleep staging[END_REF]. Even though the regression context has been less explored in the literature, it has been shown to be successful for biomarker learning, e.g., focusing on brain age as an application [START_REF] Al Zoubi | Predicting age from brain EEG signals-a machine learning approach[END_REF][START_REF] Engemann | Combining magnetoencephalography with magnetic resonance imaging enhances learning of surrogate-biomarkers[END_REF][START_REF] Sun | Brain age from the electroencephalogram of sleep[END_REF]. In the following, we focus on methods for regression modeling in the particular case of statistical discrepancies between datasets, for example, due to different populations, acquisition devices, or tasks performed during the recording. In other words, we aim to fit a regression model on one dataset and apply it on another.

Different approaches have been explored to predict cognitive-behavioral or biomedical outcomes from M/EEG data. Methods like Common Spatial Filtering (CSP) [START_REF] Koles | The quantitative extraction and topographic mapping of the abnormal components in the clinical EEG[END_REF] or Source Power Comodulation (SPoC) [START_REF] Dähne | SPoC: a novel framework for relating the amplitude of neuronal oscillations to behaviorally relevant parameters[END_REF] build on top of supervised spatial filtering for dimensionality reduction and unmixing of overlapping, yet physiologically distinct, signal generators. Lately, deep learning based techniques have been the focus of interest as they can learn good feature representation directly from the raw signal, hence potentially simplifying processing pipelines [START_REF] Roy | Deep learning-based electroencephalography analysis: a systematic review[END_REF][START_REF] Schirrmeister | Deep learning with convolutional neural networks for EEG decoding and visualization[END_REF]. Independently, an alternative approach has emerged from the BCI community which, like spatial filter methods, summarizes M/EEG data by covariance matrices. But instead of decomposing covariance matrices into filters, this approach uses mathematical tools motivated by the Riemannian geometry of the space of symmetric positive definite (SPD) matrices [START_REF] Barachant | Multiclass brain-computer interface classification by Riemannian geometry[END_REF][START_REF] Barachant | Classification of covariance matrices using a Riemannian-based kernel for BCI applications[END_REF] to define non-linear feature transformations that facilitate statistical learning with linear models. These techniques perform remarkably well given their simplicity [START_REF] Congedo | A new generation of brain-computer interface based on riemannian geometry[END_REF][START_REF] Nguyen | Inferring imagined speech using EEG signals: a new approach using riemannian manifold features[END_REF][START_REF] Sabbagh | Predictive regression modeling with MEG/EEG: from source power to signals and cognitive states[END_REF] and are competitive with methods exploiting anatomical information or endto-end deep learning approaches [START_REF] Engemann | A reusable benchmark of brain-age prediction from M/EEG restingstate signals[END_REF].

As the field of Riemannian geometry applied to M/EEG is beginning to expand and consolidate, many opportunities remain unexplored. In this work, we focus on investigating the utility of the Riemannian framework for defining dataset-harmonizing transformations.

The recent emergence of large public databases and advances in ML have led to promising prediction models. Yet, these models can be sensitive to shifts in the data distribution and may perform poorly when applied to datasets from other clinical or research contexts. We refer to these gaps between datasets as dataset shifts [START_REF] Dockès | Preventing dataset shift from breaking machine-learning biomarkers[END_REF][START_REF] Quinonero-Candela | Dataset shift in machine learning[END_REF]. This issue has also been referred to as batch effects [START_REF] Li | Harmonized-multinational qeeg norms (harmnqeeg)[END_REF]. In this work, we focused on domain adaptation techniques that attempt to deal with these shifts. We aim for a predictive model to not only perform well on the data it has been trained on, called source domain, but also when applied to data from a distinct statistical distribution, called target domain. Many domain adaptation methods exist, ranging from simple approaches minimizing the difference between the second-order statistics of source and target domains [START_REF] Sun | Correlation alignment for unsupervised domain adaptation[END_REF], to more sophisticated models measuring the distance between deep representations of the source and target domains based on optimal transport [START_REF] Damodaran | DeepJDOT: Deep joint distribution optimal transport for unsupervised domain adaptation[END_REF]. In the context of brain data analysis, Canonical Correlation Analysis (CCA) and multiway CCA (MCCA) have been largely

applied to find common brain activity and combine data across subjects when the same stimulus is presented to them [START_REF] De Cheveigné | Multiway canonical correlation analysis of brain data[END_REF][START_REF] Dmochowski | Extracting multidimensional stimulus-response correlations using hybrid encoding-decoding of neural activity[END_REF][START_REF] Lankinen | Intersubject consistency of cortical meg signals during movie viewing[END_REF][START_REF] Lankinen | Consistency and similarity of meg-and fmri-signal time courses during movie viewing[END_REF]. The experimental setups we are interested in do not meet these assumptions, for example, we consider recordings at rest and from different datasets with different subjects for which CCA or MCCA are not adapted. As we wish to work with M/EEG covariance matrices as basic signal representations for machine learning, we focus on techniques that explicitly use the geometry of SPD matrices to model the statistical distributions of distinct source and target datasets. One first ap-proach proposed in the BCI context is to re-center the distributions to a common point of the SPD space [START_REF] Li | Transfer learning based on hybrid riemannian and euclidean space data alignment and subject selection in brain-computer interfaces[END_REF][START_REF] Yair | Parallel transport on the cone manifold of SPD matrices for domain adaptation[END_REF][START_REF] Zanini | Transfer learning: A riemannian geometry framework with applications to brain-computer interfaces[END_REF]. To get a better alignment of the distributions, [START_REF] Bleuzé | Transfer learning for the riemannian tangent space: Applications to brain-computer interfaces[END_REF][START_REF] Maman | Domain adaptation using riemannian geometry of SPD matrices[END_REF][START_REF] Rodrigues | Riemannian procrustes analysis: Transfer learning for brain-computer interfaces[END_REF] propose to complement this re-centering step by adding a re-scale and a rotation correction. These covariance-based methods were initially developed to solve classification problems and are not necessarily applicable to regression without modification. In addition, most of them require labels in the target domain [START_REF] Bleuzé | Transfer learning for the riemannian tangent space: Applications to brain-computer interfaces[END_REF][START_REF] Rodrigues | Riemannian procrustes analysis: Transfer learning for brain-computer interfaces[END_REF] for alignment.

We focus on unsupervised alignment methods that can be readily used for regression modeling.

In this paper, we develop a model-based approach for tackling dataset shifts in M/EEG data in which we consider re-centering, re-scaling, and rotation techniques from previous research on classification [START_REF] Bleuzé | Transfer learning for the riemannian tangent space: Applications to brain-computer interfaces[END_REF][START_REF] Maman | Domain adaptation using riemannian geometry of SPD matrices[END_REF][START_REF] Rodrigues | Riemannian procrustes analysis: Transfer learning for brain-computer interfaces[END_REF] to regression contexts, while assuming that no labeled data is available in the target domain. We build on top of the conceptual framework from [START_REF] Sabbagh | Predictive regression modeling with MEG/EEG: from source power to signals and cognitive states[END_REF] The article is organized as follows. In Section 2, we extend the generative model from [START_REF] Sabbagh | Predictive regression modeling with MEG/EEG: from source power to signals and cognitive states[END_REF] to express and decompose dataset shifts into distinct factors, which motivates the three steps we use to compensate for dataset shifts: re-centering, re-scaling, and rotation correction. In Sections 3 and 4, we assess the robustness of these alignment steps using simulations and real-world M/EEG data.

Methods

To describe dataset shifts that can occur with M/EEG signals, we extend the generative model of M/EEG regression tasks from [START_REF] Sabbagh | Predictive regression modeling with MEG/EEG: from source power to signals and cognitive states[END_REF] where the prediction outcome is continuous. A canonical example that we will use in this work is brain age prediction [START_REF] Engemann | A reusable benchmark of brain-age prediction from M/EEG restingstate signals[END_REF][START_REF] Xifra-Porxas | Estimating brain age from structural MRI and MEG data: Insights from dimensionality reduction techniques[END_REF]. This model has also been applied to event-level regression of muscular activity with electromyogram and MEG recordings [START_REF] Sabbagh | Predictive regression modeling with MEG/EEG: from source power to signals and cognitive states[END_REF]. We describe and discuss the parameters of the generative model to understand which mechanisms can explain dataset shifts. Finally, we present various alignment strategies aiming to draw a geometrical analysis of the possible shifts and compensate for them. x i (t) ∈ R P is written as:

x i (t) = A ′ s i (t) + A ′′ n i (t) , (1)
where s i (t) ∈ R Q is the underlying signal generating this observation with Q ≤ P , and n i (t) ∈ R (P -Q) ∼ N (0, σ 2 n I P -Q ) causes a contamination due to noise. We denote A ′ = [a 1 , . . . , a Q ] ∈ R P ×Q the mixing matrix whose columns are the spatial patterns of the neural generators, and A ′′ = [a Q+1 , . . . , a P ] ∈ R P ×(P -Q) the matrix of the spatial noise patterns. Note that in this model, the noise is not considered independent across sensors but spatially correlated, as is typically the case with environmental or physiological artifacts present in M/EEG data. This model can be rewritten by combining the generator patterns and the noise in a single invertible matrix A = [a 1 , . . . , a Q , a Q+1 , . . . , a P ] ∈ R P ×P which is the concatenation of A ′ and A ′′ . The model is thus given by

x i (t) = Aη i (t) , (2)
where η i (t) ∈ R P denotes the concatenation of s i (t) and n i (t). In this model, we assume A to not depend on i nor t. It is also assumed that the statistical generators

s i (t) = {s i,j (t), j = 1 . . . Q} are zero-mean, uncorrelated,
and independent from the noise. In other words, we assume that the noise generated by artifacts is completely independent of brain activity. In the following, j will denote the generator's index.

We now consider the covariances C i of M/EEG signals X i ∈ R P ×T with T the number of time samples:

C i = X i X ⊤ i T ∈ R P ×P . ( 3
)
The covariance of M/EEG signals holds the sensors' variance on its diagonal. In our statistical model and with the previous assumptions, the covariance of the statistical generators is a diagonal matrix whose elements are the variances of each generator

E t s i (t)s ⊤ i (t) = diag(p i ) with p i ∈ R Q ,
also referred to below as "powers". Thus, we can conveniently summarize the M/EEG covariances as follows:

C i = AH i A ⊤ (4)
where

H i = E t η i (t)η ⊤ i (t) ∈ R P ×P
is a block matrix of diag(p i ) on the upper Q×Q part, and the noise covariance is in the lower (P -Q) × (P -Q) block. We here assume that E t s i (t)n ⊤ i (t) = 0, meaning that the matrix H i is block diagonal.

For regression modeling from M/EEG, it is natural to model the outcome y i as a linear combination of a function of the generators' power p i,j = E t [s 2 i,j (t)] ∈ R:

y i = β 0 + Q j=1 β j f (p i,j ) + ϵ i , (5)
where β j are regression coefficients, f is a known function, and ϵ i ∼ N (0, σ 2 ϵ ) is an additive random perturbation. For example, ageing (y) could impact brain activity in distinct brain networks (s) to different extents (β j , . . . , β Q ). This could lead for example to a log-decay or log-increase of brain activity per year, hence, motivating a logarithmic function f = log, which is a wide-spread function describing the scaling of various facets of brain structure and function [START_REF] Buzsáki | The log-dynamic brain: how skewed distributions affect network operations[END_REF]) including neural firing rates, axonal diameters, synaptic weights, and, importantly power and frequency scaling.

Replacing the generator power with the empirical average of the squared generators, the model is given by:

y i = β 0 + Q j=1 β j log 1 T T t=1 s 2 i,j (t) + ϵ i . (6)
Model violations. The assumption that A does not depend on the observation (subject) is not valid when working with actual M/EEG data. Each subject has a different head morphology, which results in slight variations in their respective mixing matrices:

A i = A+E i with E i ∈ R P ×P .
When subscript i is omitted below, A represents the average head morphology of the subjects. In our simulations below, we will assume for simplicity that E i is drawn from

N (0, σ 2 A I P ).
Riemaniann geometry basics. We are working with covariance matrices that belong to the space S ++ P of symmetric positive matrices. These matrices lie in a Riemannian manifold that can be equipped with an appropriate distance [START_REF] Förstner | A metric for covariance matrices[END_REF]:

δ R (C 1 , C 2 ) = log(C -1 1 C 2 ) F = P k=1 log 2 λ k 1 2 , (7)
where λ k , k = 1, . . . P are the real eigenvalues of C -1 1 C 2 . Note that the matrix logarithm of an SPD matrix C is computed via its eigenvalue decomposition with the log function applied to its eigenvalues: log(C) = Udiag (log(λ 1 ), . . . , log(λ P )) U ⊤ .

(8)

Similarly, the power α ∈ R + * of an SPD matrix is:

C α = Udiag(λ α 1 , . . . , λ α P )U ⊤ . ( 9
)
We can then find the geometric mean, or Riemannian mean, of a set of covariances as the minimum of the function:

C = arg min

C∈S ++ P N i=1 δ 2 R (C, C i ) . ( 10
)
Regression method. The approach we focus on in this work involves learning linear models from covariance matrices [START_REF] Barachant | Multiclass brain-computer interface classification by Riemannian geometry[END_REF][START_REF] Barachant | Classification of covariance matrices using a Riemannian-based kernel for BCI applications[END_REF]. [START_REF] Sabbagh | Manifold-regression to predict from MEG/EEG brain signals without source modeling[END_REF], 2020) show that this Riemann-based model is robust to different preprocessing choices and to model violation. This model also stands out in terms of performance when applied for regression tasks to M/EEG data in various settings.

In this framework, covariances C i are used as input of the model. The covariance matrices are vectorized to get a feature vector:

z i = uvec(S ⊙ log( C-1 2 C i C-1 2 )) ∈ R P (P +1

Possible data shifts

Each parameter of the model described in equations ( 2) and ( 6) can vary for different reasons. We are interested in fitting a regression algorithm to a source dataset X (S) =

x (S)

i , i = 1, . . . , N (S) to later predict outcomes on a distinct target dataset X (T ) = x (T ) i , i = 1, . . . , N (T ) both recorded with the same P sensors. The datasets are not necessarily composed of the same number of subjects. Applying our previous notations, we can describe the source dataset as follows:

         x (S) i (t) = A ′ (S) i s (S) i (t) + A ′′ (S) i n (S) i (t) = A (S) i η (S) i (t) y (S) i = β 0 + Q j=1 β j log 1 T T t=1 (s (S) i,j (t)) 2 + ϵ (S) i ( 12
)
where

A (S) i = A (S) +E (S) i with E (S) i ∼ N 0, σ (S) A 2 I P , n (S) i ∼ N 0, σ (S) n 2 I P -Q and ϵ i ∼ N 0, σ (S) ϵ 2 .
We remind that the statistical generator powers are defined as diag(p

(S) i ) = E t s (S) i (t)s (S)⊤ i (t) .
The same equations can be written for the target dataset by replacing the exponent (s) by (t). We now list physical reasons that could induce differences between source and target datasets and link them to parameter changes of the corresponding generative models (12).

1. If we consider two different populations, the head morphology may vary, and the subject-averaged mixing matrices A (S) and A (T ) would differ.

2. Having different populations in both datasets would also imply that they will not have the same mixing matrices distribution: σ (S)

A ̸ = σ (T )
A . 3. When data are recorded with different devices, the recording conditions and noise might not be the same, resulting in different signal-to-noise ratio (SNR):

σ (S) n ̸ = σ (T ) n .
4. Clinical outcomes e.g., neuropsychological testing scores can be noisy. This noise could differ from one dataset to another: σ

(S) ϵ ̸ = σ (T ) ϵ .
Because of all those possible causes of variability in the model parameters, machine learning approaches may fail to provide good predictions across datasets. In this work, we focus on shifts that only affect the data, and we assume that the regression coefficients β j are the same for source and target. In particular, we are interested in understanding changes related to the mixing matrix and the variance of the statistical generators. The variability of these parameters across subjects and datasets affects the observed signals and results in variability in the data distribution.

Below, we discuss which statistical methods could help reduce these different shifts between the data distributions of two different datasets.

Alignment methods

We aim to learn a regression model from one dataset, the source domain, that will perform well on another, the target domain. As we focus on shifts affecting the data distribution, we investigate domain adaptation methods that align the source and the target distributions using geometrical transformations. The methods we chose for understanding and reducing dataset shifts are articulated in three alignment steps: re-centering, equalizing dispersion, and rotation correction. This choice was inspired by transfer learning methods used in Brain-Computer Interfaces (BCI) application and, more specifically, by the Riemannian Procrustes Analysis (RPA) of [START_REF] Rodrigues | Riemannian procrustes analysis: Transfer learning for brain-computer interfaces[END_REF]. These steps can be used independently, usually by only re-centering the data, or combined. In the following, we detail these alignment functions in a general manner.

Step 1: re-centering. The most commonly used method of transfer learning on symmetric positive definite (SPD) matrices is to re-center each dataset in a common reference point on the Riemannian manifold [START_REF] Li | Transfer learning based on hybrid riemannian and euclidean space data alignment and subject selection in brain-computer interfaces[END_REF][START_REF] Zanini | Transfer learning: A riemannian geometry framework with applications to brain-computer interfaces[END_REF]. This reference can be chosen as one of the domains' geometric mean or an arbitrary point on the manifold. Here we propose to re-center each domain to the Identity by whitening them with their respective geometric mean C(S) and C(T ) :

h rct C(S) C (S) i = C(S)-1 2 C (S) i C(S)-1 2 (13) h rct C(T ) C (T ) i = C(T )-1 2 C (T ) i C(T )-1 2 . ( 14
)
Differently put, re-centering applies separate whitening for source versus target data. This helps avoid errors in the tangent space projection when the average covariance is different for the source and target, e.g., because the mixing matrices are different, as in Figure 1 (B). This is a

Riemannian equivalent of the centering step in classical z-scoring.

Step 2: equalizing the dispersion. In this second step, the idea is to re-scale the covariances distribution around their mean C as illustrated in Figure 1(C). We first compute the mean dispersion d of the covariances as the sum of the square distance between each matrix of the set and their geometric mean C over the number of samples in the dataset:

d (S) = 1 N N i=1 δ 2 R C (S) i , C(S) (15)
d (T ) = 1 N N i=1 δ 2 R C (T ) i , C(T ) , (16)
with δ R the Riemannian distance defined in Equation 7.

Then, we re-scale all covariances with 1 d so that the distribution has a dispersion of 1:

h str d (S) C (S) i = C (S) 1 d (S) i (17) h str d (T ) C (T ) i = C (T ) 1 d (T ) i . ( 18
)
This is a Riemannian equivalent of the re-scaling step in classical z-scoring. By analogy, univariate rescaling of two groups modifies the data so that, e.g. a t-test would find its assumptions of equal variances met while not detecting any difference of means between the two datasets.

Step 3: rotation correction. Until now, we have applied correction measures that process source and target data independently. This is not the case in this third step which implies shared information between source and target datasets. The rotation correction is the most delicate of the three steps. It requires estimating many more parameters than the others, and the source and target feature spaces must be the same size (P (S) = P (T ) = P ). In the literature, several methods for rotation estimation exist.

In the following, we detail two of the methods we selected:

1. The first rotation correction we implemented is inspired by [START_REF] Maman | Domain adaptation using riemannian geometry of SPD matrices[END_REF]. The covariances are first vectorized by mapping them in the tangent space (at identity after re-centering). Then we compute the Singular Value Decomposition (SVD) of these tangent vectors S) and T ) :

z i ∈ R P (P +1)/2 , Z (S) = z (S) i , i = 1 . . . N (
Z (T ) = z (T ) i , i = 1 . . . N (
Z (S) = U (S)⊤ S (S) V (S) (19)
Z (T ) = U (T )⊤ S (T ) V (T ) . ( 20
)
The columns of the U ∈ R P (P +1)/2×P (P +1)/2 matrices are the left singular vectors ordered from largest to smallest singular values. The SVD is done separately on source and target distributions so the resulting singular vectors will unlikely have the same direction. As we desire for corresponding singular vectors between source and target to have an acute angle, we reorient them with a sign correction applied to the columns of U (T ) :

u (T ) j = sign u (S) j u (T ) j u (T ) j , ∀j (21)
Finally, the U matrices are used for rotation correction:

h rot U (S) Z (S) = U (S)⊤ Z (S) (22) h rot U (T ) Z (T ) = U (T )⊤ Z (T ) . ( 23
)
We will refer to this rotation correction method as unpaired.

2. The second way to estimate the rotation that we used is inspired by [START_REF] Bleuzé | Transfer learning for the riemannian tangent space: Applications to brain-computer interfaces[END_REF]. In this paper, they consider a classification question and propose to correct the rotation between source and target distributions by matching their respective classes' mean.

This is done by solving the Procrustes problem arg min

R∈O(N (T ) ) R Z(T ) - Z(S) F , (24)
where O(N (T ) ) is the orthogonal group, Z(S) the concatenation of the classes' mean tangent vector from the source domain, and similarly for Z(T ) with the available labeled data of the target domain. Then, to correct the rotation, the target tangent vectors are transformed using the solution of the Procrustes problem

h rot R Z(T ) = R Z(T ) . (25)
As we wish to be in a regression context without access to target labels, we modified this method by solving the Procrustes problem on all the tangent vectors arg min

R∈O(N (T ) ) RZ (T ) -Z (S) F . (26)
In practice, this solution is found by computing the SVD of the product of the source and target tangent vectors

Z (T ) Z (S)⊤ = U ⊤ SV , (27)
and is R = VU ⊤ . In this step, we include the information on which source point should be matched to which target point. It means the source and the target dataset should have the same number of observations (N (S) = N (T ) = N ) and be composed of "matching" observations (for example, the same set of subjects but different tasks/recording conditions/devices). We will refer to this method as paired.

This step allows us to align source and target distributions in a shared space. The rotation correction is helpful when the mixing matrices are different between the domains (Fig. 1 (D)).

Benchmarks

We conducted a first benchmark with simulated data to evaluate how the alignment steps can compensate for changes in the generative model parameters. We then applied the alignment steps to MEG data from the same dataset. In this benchmark, we used the different tasks performed by the subjects as domains. We first evaluated the methods using the same subjects in source and target but performing different tasks, and then used some subjects performing one task as the source domain and other subjects performing another task as the target domain. Finally, we evaluated whether the alignment makes it possible to generalize from one M/EEG dataset to another. The regression problem of these benchmarks is age prediction. It means we have no noise on the outcome as assumed in Section 2.2. More details about these experimental set-ups can be found in Section 3.2.3. In the following, we detail precisely the different setups.

Simulations

In the simulation study, we generated simulated data with the generative model presented in Section 2.1. We set the dimension of the matrices to P = 20 to have a matrix size coherent with real data, and the number of statistical generators to Q = P , in other words, we considered signals without noise. Mixing matrices A were generated as Gaussian random matrices in R P ×P from N (0, 1). Instead of generating signals s, we directly computed their powers p as random numbers from a uniform distribution in [0, 1). The same powers were used for both the source and the target sets. We then constructed the covariance matrices

C (S) = C (S) i , i = 1, . . . , N and C (T ) = C (T )
i , i = 1, . . . , N , and the outcome to predict as in equations ( 4) (with

H i = diag(p i ) because P = Q)
and ( 6) (with ϵ i = 0, ∀i). We designed several shift scenarios by altering either the mixing matrices or the powers in order to evaluate the alignment methods. In practice, the shifts were created by first generating the source data and then building the target data as a modified version of the source data. By doing so, one point of the source set corresponds to one point of the target set. This way, it is possible to evaluate the paired rotation correction.

More details about the shift scenarios are presented in the following paragraphs.

Simulation scenarios

We detail the changes we introduced for each scenario between source and target distributions. As stated in Section 2.2, we focused on modeling shifts involving changes in the mixing matrix or the variance of the statistical generators. The parameters that are not mentioned were the same for source and target. For the first three scenarios, we aimed to find transformations/shifts between source and target to which each alignment step is robust. Table 2 summarizes the scenarios and the associated parameter changes in source and target.

Translation. In this scenario, we wish to assess how robust the alignment methods are when the source and target mixing matrices are different. Specifically, we built the target mixing matrix as

A (T ) = B α A (S) with B ∈ S ++ P and α ∈ R + * .
Here, one considered that the mixing matrix perturbation is done by an SPD matrix to decompose the case of A (T ) ̸ = A (S) into translation and rotation effects. The benchmark extends the previous simulations from [START_REF] Sabbagh | Predictive regression modeling with MEG/EEG: from source power to signals and cognitive states[END_REF]. More explanations

Scenarios Translation

Scale

Translation and rotation Noise on mixing matrix Source A (S) fixed p

(S) i fixed A (S) fixed, A (S) i = A (S) + E (S) i with E (S) i ∼ N 0, σ (S) A 2 I P σ (S) A = 10 -2 fixed Target A (T ) = (B) α A (S)
with α > 0 are provided about this decomposition in the translation and rotation paragraph. The parameter α controls the strength of the perturbation and thus how

p (T ) i = (p (S) i ) σp with σ p > 0 A (T ) = mA t + (1 -m)A (S) with A t ̸ = A (S) fixed and m ∈ [0, 1] A (T ) i = A (T ) + E (T ) i with E (T ) i ∼ N 0, σ (T ) A 2 I P σ (T ) A > 0
A (T ) is different from A (S) (if α = 0, A (T ) = A (S) ).
Scale. We wanted to create a scenario in which the source and target distributions have different dispersions. In this scenario, we constructed the target covariances with an exponent on the powers: p

(T ) i = p (S) i σp with σ p > 0.
The parameter σ p controls how different the dispersions are. This modification was only applied to the data, so the outcome values y were unchanged.

Translation and rotation. For this scenario, we built the source and target data from completely different mixing matrices and thus generalized the translation scenario. To evaluate how alignment methods performed for a growing difference between the source and the target mixing matrices, we defined a parameter m such as S) . A t was fixed and generated as a random matrix in R P ×P from N (0, 1). In this manner, we created an interpolation between A t and A (S) to generate A (T ) : if m = 0, A (T ) = A (S) and if m = 1,

A (T ) = mA t + (1 -m)A (
A (T ) = A t .
In this scenario, A (T ) ̸ = A (S) but we constructed the source and target covariances with the same H i matrices following equation ( 4). Thus we can write:

H i = [A (S)⊤ ] -1 C (S) i [A (S) ] -1 (28)
We can replace this H i expression in the target covariances to get: Noise on mixing matrix. We finally introduced individual noise in the mixing matrix to get a more realistic scenario: and A (T ) were the same. This scenario was inspired by the simulation study of [START_REF] Sabbagh | Predictive regression modeling with MEG/EEG: from source power to signals and cognitive states[END_REF] in which the same level of noise on the mixing matrix was added in the train and the test sets. Here, we explored the situa-tion in which the noise levels in the train (source) and test (target) mixing matrices were different.

C (T ) i = A (T ) H i A (T )⊤ (29) = A (T ) [A (S)⊤ ] -1 C (S) i [A (S) ] -1 A (T )⊤ (30) = DC (S) i D ⊤ (31) with D = A (T ) [A (S)⊤ ] -1 .
A (S) i = A (S) + E (S) i assuming that E (S) i ∼ N 0, σ (

Alignment, vectorization, and regression

Once the covariance matrices were generated according to a given scenario, data of both domains were aligned with the methods detailed in Section 2.3. Then, we vectorized the matrices in the tangent space as in (11) with C(S) as a reference point for both domains. To avoid numerical issues, we removed low-variance features (see Appendix A for more details). The remaining features were then standardized to get features with zero mean and unit variance.

To predict from the standardized vectors in these simulations, for simplicity, we used Ridge regression with its regularization term set to 1. This model was trained on the source data, and predictions were made on the target data. We evaluated these predictions with R 2 scores.

Results are presented in Section 4.1 and discussed in Section 5.

M/EEG empirical benchmarks

In the following empirical benchmarks, we focused on one MEG and two EEG datasets for evaluating our alignment methods with real-world data. We first describe these datasets and their preprocessing, then explain how we computed the covariance matrices of the signals. We followed the same preprocessing and processing steps as in [START_REF] Engemann | A reusable benchmark of brain-age prediction from M/EEG restingstate signals[END_REF] for the 'filterbank-riemann' pipeline. Finally, we detail the design of each benchmark.

Datasets

Cam-CAN MEG data. The Cambridge Center of Aging and Neuroscience (Cam-CAN) dataset [START_REF] Taylor | The cambridge centre for ageing and neuroscience (cam-CAN) data repository: Structural and functional MRI, MEG, and cognitive data from a cross-sectional adult lifespan sample[END_REF] consists of MEG recordings from a healthy population covering a wide age range. These data were recorded for each subject during resting state with eyes closed, an audio-visual (passive) task with visual and auditory stimuli presented separately, and a sensorimotor (smt) task with the same stimuli as the previous task combined with a manual response. All data were collected with the same 306-channel VectorView MEG system (Elekta Neuromag, Helsinki) with a sampling rate of 1 kHz.

Sample description:

We included 646 subjects (319 female, 327 male) with all three recordings. Their age distribution is from 18.5 to 88.9 years with an average of 54.9 ± 18.4 years and an almost uniform spread over the age range. There was no exclusion of participants. The set of subjects of each benchmark only depends on the availability of recordings for the source and the target tasks and the success of the preprocessing and the feature extraction.

Thus some subjects with only two tasks recorded are not included in all benchmarks leading to small variations of the subject sample between benchmarks.

Preprocessing: We applied a FIR band-pass filter between 0.1 and 49 Hz to all data. We decimated the signals with a factor of 5 to get a sampling frequency of 200 Hz. To compensate for environmental noise, we performed a temporal signal-space-separation (tSSS) method [START_REF] Taulu | Applications of the signal space separation method[END_REF] with a chunk duration of 10 seconds and a correlation threshold of 98%. We only picked channels corresponding to magnetometers (after tSSS signals from magnetometers and gradiometers are mixed and linearly related).

TUAB EEG data. The Temple University (TUH) EEG Corpus [START_REF] Harati | The TUH EEG CORPUS: A big data resource for automated EEG interpretation[END_REF]

Preprocessing:

A band-pass filter between 0.1 and 49 Hz was applied to the data and resampled to 200 Hz.

To keep a maximum of data, recordings with eyes closed and eyes open were pooled before feature extraction.

Data processing and feature extraction

After preprocessing, each filtered recording was segmented in 10 s epochs without overlap. Epochs were then filtered into seven frequency bands as defined in Table 3 as in [START_REF] Engemann | A reusable benchmark of brain-age prediction from M/EEG restingstate signals[END_REF]. We performed artifact rejection by thresholding extreme peak-to-peak amplitudes on single epochs using the local autoreject method [START_REF] Jas | Autoreject: Automated artifact rejection for MEG and EEG data[END_REF]. Subsequently, we computed covariance matrices from the set of artifact-free epochs with the Oracle Approximating Shrinkage (OAS) estimator [START_REF] Chen | Shrinkage algorithms for MMSE covariance estimation[END_REF]. The ensuing regression pipeline, including all alignment steps, is illustrated in Figure 2.

For MEG signals, the tSSS method reduces noise by projecting them in a subspace mainly containing the signal, leading to rank-deficient covariance matrices. As a result, it is not possible to correctly apply our alignment methods directly, as rank-deficient covariance matrices are not SPD matrices. To extract valid SPD matrices, we follow the approach from [START_REF] Sabbagh | Manifold-regression to predict from MEG/EEG brain signals without source modeling[END_REF] and apply Principal Component Analysis to reduce the dimensionality of the covariance matrices, which renders them full rank. We denote the resulting dimension R. We use as filters the eigenvectors W (S) ∈ R R×P corresponding to the R highest eigenvalues of the mean source covariance C(S) .

Matrices of both domains are transformed as:

C (S) i,filtered = W (S) C (S) i W (S)⊤ ∈ R R×R (32) C (T ) i,filtered = W (S) C (T ) i W (S)⊤ ∈ R R×R (33)
For Cam-CAN data, we set R = 65. This spatial filter is applied to each frequency band separately. We did not apply this procedure to EEG data.

Preprocessing Covariance computation

Covariance computation

Projection and vectorization

Projection and vectorization

Concatenation Regression

Re-center Re-scale

Filtered epochs

Covariance matrices

Tangent vectors

Re-center Re-scale summarize the M/EEG recording by the covariance matrix after performing artifact cleaning (Section 3.2.1). The covariances computation, alignments steps, projection to the tangent space, and vectorization steps are done separately for seven frequency bands of Table 3. Alignment steps detailed in Section 2.3 are computed from the covariance distribution across all subjects. The re-center and re-scale steps are performed separately for source and target datasets. The Procrustes steps combine information across source and target datasets. Finally, the seven resulting tangent vectors are concatenated to form one vector per subject used for regression.

Empirical benchmarks

Cam-CAN (MEG): same subjects. For this benchmark, we used the experimental tasks of the Cam-CAN dataset for defining the different domains. Here, the source subjects were the same as the target subjects. Only the experimental task (e.g., audiovisual VS audiovisual + motor) changed from one domain to the other. All subjects included underwent MEG recordings for both the source and the target domain. Therefore, importantly, the mixing matrix and the age distributions were the same for the source and target. As we dealt only with healthy participants, we aimed to minimize the error in age prediction when learning on one task (source domain) and predicting on another (target domain). To estimate standard deviation, we did a bootstrap with 100 repetitions.

Cam-CAN (MEG): different subjects. In this second benchmark on the Cam-CAN data, we again defined the different domains in terms of the experimental MEG tasks performed by the subjects. Yet, the critical difference with the previous benchmarks is that the source subjects and the target subjects were distinct persons. To implement this analysis, we randomly divided all Cam-CAN subjects into subsets of 80% forming the source subjects, and the left-out 20 % forming the target subjects. A stratification was performed by age decade to maintain similar age distributions between splits. We repeated this split with 100 different random initializations.

TUAB & LEMON (EEG): different datasets. In this benchmark, we gauged the performance of alignment methods when the source and target domains are two different datasets.

Here, the source domain was composed of data from TUAB, and the target one of data from LEMON. These datasets were not recorded with the same device. However, they had 15 channels in common. We picked the same channels on both datasets to define covariance matrices of the same shape and similar information. The target set was kept fixed, and we implemented a bootstrap procedure on the source subjects to estimate standard deviations. In this

z i = diag(log(W SPoC C i W ⊤ SPoC )) (34)

Alignment, vectorization, and regression

The matrices from both domains were first aligned with the methods described in Section 2.3. We projected the aligned data in the tangent space at C(S) to get tangent vectors. Tangent vectors from all frequency bands were concatenated. Then we applied ridge regression after standardizing (z-scoring) all the features. To select the regularization hyperparameter, we used a generalized (leave-oneout) cross-validation [START_REF] Golub | Generalized crossvalidation as a method for choosing a good ridge parameter[END_REF] on a logarithmically spaced grid of 100 points from 10 -5 to 10 10 . For quantifying prediction performance, we use the R 2 score.

Software

We processed the M/EEG data with the open-source package MNE-Python [START_REF] Gramfort | MEG and EEG data analysis with MNE-Python[END_REF], MNE-BIDS [START_REF] Appelhoff | MNE-BIDS: Organizing electrophysiological data into the BIDS format and facilitating their analysis[END_REF], and the associated MNE-BIDS-Pipeline (https://mne.tools/mne-bids-pipeline/). The covariance matrices computation and the predictive pipeline were done with the coffeine library [START_REF] Sabbagh | Predictive regression modeling with MEG/EEG: from source power to signals and cognitive states[END_REF] (https://github.com/coffeine-labs/coffeine). The covariance matrices were manipulated in the alignment methods with the Pyriemann package [START_REF] Barachant | pyriemann/pyriemann: v0[END_REF]. Results analyses were performed with the Scikit-Learn software [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF].

Results

We now present the results we obtained from the simulation and M/EEG benchmarks. We first computed the baseline performance of the Riemannian framework without alignment. We then added one alignment step at a time to evaluate its impact on prediction performance. For example, the re-scale method corresponds to re-centering and re-scaling (step 1 + step 2). Likewise, steps 1 and 2 were always performed before rotation correction (Procrustes unpaired and Procrustes paired).

We implemented an element-wise z-scoring to provide comparisons between this widely used method [START_REF] Apicella | EEG-based measurement system for monitoring student engagement in learning 4.0[END_REF][START_REF] Chen | Personal-zscore: Eliminating individual difference for EEG-based cross-subject emotion recognition[END_REF], which is not ideal as it ignores the covariance manifold that can be described and

Simulations

We simulated data according to the four shift scenarios detailed in Section 3.1. Figure 3 presents the results for each alignment method and each scenario. ) with B ∈ S + P . All methods that included re-centering the distributions on the same reference point performed well. (B) displays the performance achieved when the dispersion of covariances differs between source and distributions (σp ̸ = 1). Here, the re-scaling step was essential to align the distributions correctly. (C) In this scenario, A (S) ̸ = A (T ) , which led to a translation and a rotation of the target set compared to the source set. Re-centering was not insufficient, and a rotation correction was needed to achieve good performance. Interestingly, while Procrustes paired performed well, the unpaired correction broke as the difference between the mixing matrices increased. (D) In this scenario, different levels of individual noise were added to the mixing matrices of both domains. For low σ (T ) A values, all methods except the unpaired rotation correction performed similarly with R 2 scores decreasing slowly. For higher values, the scores dropped, and correcting the rotation with the paired method performed best.

A (T ) = B α A (S
The top left Panel (A) illustrates the scores for each alignment method on data generated following the translation scenario. The value of α controls the shift. As expected, the farther apart the source and target mixing matrices were (higher α values), the worse the performance on the unaligned domains method became (in blue). The z-score baseline (light orange) failed even earlier than using Riemannian geometry without alignment. Methods including a re-centering step (green, dark orange, pink, and brown) did not suffer from this shift. This suggests that whitening the source and target distributions by their respective geometric mean mostly compensated for the mixing matrix being perturbed by an SPD matrix. It allowed the regression model to access the log of the powers with little distortion, hence, allowing the linear model to infer the correct function.

Panel (B) presents the scale scenario in which the log of the target powers were scaled by a parameter σ p in the signals. When σ p = 1, the source and target distributions were exactly the same C (S) = C (T ) . In this case, methods including the re-scaling step adjusted for this shift and made accurate predictions, whereas the performance of other methods deteriorated as σ p increased. Re-centering helped to achieve better predictions compared to no alignment. The z-score method performed slightly better than not aligning the distribution but was still worse than recentering.

The third Panel (C) corresponds to the translation and rotation scenario. Here, the target mixing matrix was modified by interpolating between the source mixing matrix A (S) and another randomly generated matrix A t . The parameter m controls where the target mixing matrix is located between these two other matrices, thus how different A (S) and A (T ) were. The only method reaching perfect predictions, irrespective of the value of m, was Procrustes paired. However, the unpaired Procrustes method failed where m > 0.5 and even fell behind re-centering. When no rotation correction was applied, re-centering helped to compensate for slight differences between the mixing matrices, but the performance dropped as this difference increased. As expected, a re-centering step and a rotation correction were needed to correct a shift consisting of translation and rotation.

In the last scenario noise on mixing matrix, displayed in Panel (D), we introduced noise in both source and target mixing matrices to simulate individual differences between subjects. σ (S) A was set at 10 -2 , and even when σ

(T ) A = σ (S) A we had A (T ) i ̸ = A (S)
i . Procrustes unpaired performed worst in this scenario. The unpaired rotation correction was not robust to noise on the mixing matrix. All the other methods performed similarly for low values of σ (T )

A . When σ (T ) A > 10 -1 , all methods deteriorated. The z-score method again showed lower R 2 scores than all other methods. Results suggest that the best solution for this scenario is the paired rotation correction.

The paired rotation correction method performed best in all scenarios but requires the source and target sets to be the same size and have corresponding/paired points.

When this is not the case (for datasets with different subjects, for example), re-centering and re-scaling should be the best solution for improving performance. The unpaired rotation estimation seems particularly unstable when the induced shift is too big or when there is noise.

M/EEG data

We now examine the performance of alignment methods with M/EEG data. As with simulated data, the unpaired Procrustes method was not sufficiently robust and led to chance-level performance. We, therefore, do not report it in the following figures.

Cam-CAN (MEG): same subjects. We first focused on a regression problem for which there was an individual noise on the mixing matrices, but their distribution was the same in source and target sets because the subjects were the same. It also implies that the age distribution was identical for both domains. The interest of this analysis The regression model was trained on resting-state data, and predictions were made on the recordings of the somatosensory task. Re-centering the data led to slightly improved R 2 scores. Again, the re-scaling step did not lead to further improvements.

Correcting the rotation with the paired method contributed to improving 99% of the splits in comparison to only re-centering. (C) Here we used the data from the passive task as the source domain and the somatosensory task as the target domain. Re-centering and re-scaling steps did not affect the prediction performance. The paired rotation correction improved the scores in all splits.

is to assess what kind of shift is produced when only the task changes and if alignment methods can rectify this.

Results for each alignment method on age prediction for three source-target tasks associations are displayed in Fig- ure 4.

The z-score method led to scores similar or lower to what is obtained without alignment across all three Panels, as expected from the simulation results. For the two first Panels (A and B), the source domain contained the resting-state recordings, and the target tasks were, respectively, the passive and the somatosensory tasks. The R 2 scores we obtained after these two benchmarks were highly similar. When no alignment was done, the mean R 2 score was around 0.7. Re-centering the distributions led to a reduced standard deviation and an increased mean R 2 score in both situations, even though this was more pronounced in Panel (A). The re-scaling step had no obvious impact on performance. The paired rotation correction led to improved prediction scores on 87.4% of the bootstrap iterations in Panel (A) and on 99.1% of the iterations in Panel (B) compared to only re-centering. In Panel (C), the source domain was the passive task, and we made predictions on the somatosensory task, leading to quite different results. The performance reached with no alignment was already very high, with a mean R 2 of 0.9. Re-centering and re-scaling gave the same results as not performing any alignment. Then, the paired rotation correction step induced increased scores for all bootstrap iterations (R 2 = 0.951 ± 0.005). Going from rest to tasks affects the geometric mean of the covariance matrix distributions, but not when going between passive and smt tasks. In all three situations, the performance gain obtained with Procrustes paired implies the presence of a rotation of the tangent vector distribution. Panel, the passive task was the source domain, and the somatosensory task was the target. In this case, aligning was not helpful and led to the same performance as not performing any alignment.

Cam

tasks in non-overlapping samples of distinct subjects. As a consequence, the distributions of mixing matrices, necessarily, differed for the source and target domains. We performed a cross-validation in which 80% of the subjects were assigned to the source domain and 20% to the target domain. To keep relatively similar age distributions in the train and test splits, we did a stratification on the age decades of the subjects (cf StratifiedShuffleSplit of the Scikit-Learn software). Applying the paired rotation correction in this setup was no longer possible. Thus, it is impossible to analyze whether a rotation exists in the shift. We present the results of each alignment step in Figure 5.

When covariance matrices were not aligned, generalizing from rest to passive tasks led to a R 2 score of 0.55±0.05 (A), and when the target task was the smt task, we observed an R 2 = 0.54 ± 0.04 (B). The z-score method performed again similarly to the procedure without align-ment. The re-centering step led to comparable results across generalization scenarios involving resting state and any event-related task (Panels (A) and (B)). Again, matching the source and target dispersions with re-scaling was not helpful. Finally, all methods showed similar performance when the passive and the smt tasks were the source and the target tasks, respectively (C). Our observations for this benchmark match those we made when the subjects were the same for the source and target sets. The R 2 scores reached after alignment in Figure 5 11(as for all the results we have previously presented). Without alignment, the prediction made on the LEMON data led to R 2 scores far lower than what was reported in [START_REF] Engemann | A reusable benchmark of brain-age prediction from M/EEG restingstate signals[END_REF] (10-fold cross-validation on LEMON data only: 0.54 ± 0.13 represented by the dashed gray line). When both domains are re-centered to identity, we reached performances similar to when the model is trained on LEMON. Re-scaling did not visibly improve results. (B) Results when the regression model follows the SPoC approach. Not aligning led again to poor R 2 scores. Unlike the first panel, the z-score method improved the predictions similarly to re-centering. Re-scaling helped to reach performances on par with the Riemannian model trained on LEMON.

and the target dataset is the smaller LEMON dataset (n=213). TUAB also has a broader age range. This way, the regression model will be asked to predict age values that fall within the range observed during model training.

We performed a bootstrap with 2000 iterations on TUAB data. The results are reported in Figure 6. In addition to the Riemannian approach we focused on in this work (A), we were also interested in the impact of the alignment methods on a non-Riemannian model like SPoC (B).

Without alignment, the Riemannian model and SPoC led to poor results with mean R 2 scores around 0.26. On Panel (A), the z-score method performed at the level of the dummy model. Re-centering the data drastically improved the age prediction performances with R 2 scores of 0.44 ± 0.06 with a visibly reduced variance. Adding the re-scaling step on top of re-centering did not bring any improvement in performance. In [START_REF] Engemann | A reusable benchmark of brain-age prediction from M/EEG restingstate signals[END_REF]

Discussion

In this study, we thoroughly explored domain adaptation methods that align M/EEG covariance matrix distri-butions for regression problems on both simulations and large datasets. We considered methods from BCI applications [START_REF] Bleuzé | Transfer learning for the riemannian tangent space: Applications to brain-computer interfaces[END_REF][START_REF] Maman | Domain adaptation using riemannian geometry of SPD matrices[END_REF][START_REF] Rodrigues | Riemannian procrustes analysis: Transfer learning for brain-computer interfaces[END_REF] et al., 2022). The re-centering step is simple to implement and has already been very effective in BCI classification to deal with variability between sessions [START_REF] Barachant | Classification of covariance matrices using a Riemannian-based kernel for BCI applications[END_REF] but also between subjects [START_REF] Zanini | Transfer learning: A riemannian geometry framework with applications to brain-computer interfaces[END_REF]. Our results suggest it is also effective in a regression context with variability between populations, tasks, and recording devices.

We extended our evaluation of the impact of alignment methods on different EEG datasets to the SPoC model [START_REF] Dähne | SPoC: a novel framework for relating the amplitude of neuronal oscillations to behaviorally relevant parameters[END_REF]. In this setting, the z-score method and re-centering performed both equally better than no alignment. Interestingly, re-scaling was beneficial and helped to reach performance close to the Riemannian model trained on LEMON. By inspecting the SPoC patterns and the associated log powers, we demonstrated that the observed gain in the performance of re-centering was enabled by more similar log powers between source and target than without alignment. In other words, data alignment adapts the target features to the regression equations fitted on the source data, which explains generalization.

Unfortunately, the two last benchmarks are missing a rotation correction method. As Procrustes paired led to an apparent score increase in the first benchmark, we ex-pect a rotation correction to be beneficial in the other benchmarks. Yet, the simulation study showed that the unpaired Procrustes method failed to correctly estimate the rotation when there is noise or when the shift gets too large. The result suggests that this method would likely fail with M/EEG data. The condition of matching subjects between the source and target in Procrustes paired is too restrictive and is not applicable in many settings.

The supervised rotation correction methods developed for classification problems [START_REF] Bleuzé | Transfer learning for the riemannian tangent space: Applications to brain-computer interfaces[END_REF][START_REF] Maman | Domain adaptation using riemannian geometry of SPD matrices[END_REF][START_REF] Rodrigues | Riemannian procrustes analysis: Transfer learning for brain-computer interfaces[END_REF] are unsuitable for regression.

Further investigations are needed to fill the lack of rotation correction in regression contexts.

Another limitation of this work is that we only performed benchmark that involved source and target covariance matrices formed from the same set of sensors. The dimensionalities of the source and the target data must be equal to apply the predictive model to it. It has been proposed to deal with different dimensionalities of covariance matrices via zero-padding [START_REF] Rodrigues | Dimensional-ity transcending: a method for merging bci datasets with different dimensionalities[END_REF].

However, this method is not applicable if there is no rotation correction afterward, so we could not use it. We also observed that our framework is not robust to sensor permutation in the covariance matrices, even with the same sets of sensors. In our EEG benchmark, we had to select the common channels between the two datasets and sort them to reach acceptable performance even for the re-center step. In addition to leading to a gain in performance, having a proper rotation correction would help to deal with issues related to different numbers or types of sensors in the source and target datasets.

Besides the limitations linked to the rotation correction, further points would deserve future studies. First, we explored unsupervised alignment methods that do not explicitly share any information between the source and the target domains. Comparing our results with supervised methods could allow us to quantify the gain of supervision and to have additional insights into the trade-off between the approaches. A second element to consider is that our goal was to evaluate the alignment methods on a regression problem by minimizing the prediction error. We focused on brain age prediction as age is a label that is easy to collect. But other prediction targets should equally benefit from the methods presented in this work. Importantly, we conducted our benchmarks on healthy participants sampled from the general population. Yet, the biggest impact of our results may be seen when bridging datasets from heterogeneous clinical populations, which remains to be demonstrated. Finally, in this work we focus on linear regression model but recent work demonstrated that this Riemannian framework can also be applied with nonlinear models [START_REF] Bonet | Sliced-Wasserstein on symmetric positive definite matrices for M/EEG signals[END_REF]. Kernel-based models have been shown to perform well on brain age prediction but were not investigated in a domain adaptation context.

Figure 1 :

 1 Figure 1: Alignment steps illustrated on simulated data. The three alignment steps are applied to data simulated following the generative model, as detailed in Section 3.1. We set the size of the matrices to P = 2 and generated N = 300 matrices in each domain. Each new step is applied on top of the previous one. The plots correspond to the two first principal components of the tangent vectors. (A) The simulated data are plotted on the tangent space before any alignment steps. (B) The original simulated data are centered to a common point, (C), then their distributions are equalized, and (D) finally, a rotation correction is applied.

  The target covariance matrices correspond to the source covariance matrices transformed with the square matrix D. A square matrix can be interpreted as a linear transformation: such a matrix can be decomposed into the product of an orthogonal matrix with a positive semi-definite Hermitian matrix (polar decomposition, a.k.a. QR factorization). Thus, we can interpret this scenario as the translation scenario (SPD matrix of the polar decomposition) with an additional perturbation by an orthogonal matrix.

  -3 to 1. Here, the mean mixing matrices A(S)

  of 1385 subjects (female = 775 and male = 610) with ages between 0 and 95 years (mean = 44.4 years and std = 16.5 years). Preprocessing: Data were band-pass filtered between 0.1 and 49 Hz with a zero-phase finite impulse response (FIR) filter using the firwin with Hamming window featuring a 0.0194 passband ripple, 53 dB stopband attenuation, a 0.1 Hz lower transition bandwidth and a 12.25 Hz upper transition bandwidth, and a filter length of 6601 samples (33.005 s). Data were then resampled to 200 Hz. We selected a subset of 21 channels common to all recording devices used in this dataset. When several recordings were available for one patient, we picked the first to get only one recording per subject. LEMON EEG data. The Leipzig Mind-Brain-Body database provides multimodal data from healthy groups of young and elderly subjects (Babayan et al., 2019). In our benchmark, we only used EEG recordings from this dataset. They were recorded with a 62-channel ActiCAP device and sampled at 2500 Hz. Each subject did two recordings at rest with two conditions: eyes closed and eyes open. Sample description: We included 213 subjects from the LEMON database in our benchmark. No selection criteria were applied, and we kept the data for which the processing and the feature extraction were successful. This led to a cohort with 134 males and 79 females aged from 20 to 77 years. The age distribution of the LEMON presents a peculiarity: it is split into two separate age groups, one with individuals being between 20 and 35 years old and the second between 55 and 77 years old.

Figure 2 :

 2 Figure 2: Pipeline for regression modeling with M/EEG with different dataset harmonization steps. For every subject, we

  setup, in addition to evaluating the alignment methods on the Riemannian regression model, we also applied them with a regression model based on Source Power Comodulation (SPoC). SPoC is a supervised spatial filtering method in which the filters W SPoC maximize the covariance between the power of the filtered signals and the outcome y (Dähne et al., 2014). Denoting by C = 1 N N i=1 C i the Euclidean average covariance matrix and C y = 1 N N i=1 y i C i the weighted average covariance matrix, the first filter w SPoC is given by: w SPoC = arg max w w ⊤ Cyw w ⊤ Cw . The same idea was proposed by (de Cheveigné and Parra, 2014). As W SPoC recovers the inverse of the mixing matrix A (Sabbagh et al., 2020), the SPoC regression model is defined as:

  handled by Riemannian geometry, and the presented alignment methods. For this benchmark, we transformed the covariance matrices into correlation matrices by removing the variance of each observation/subject individually. This way, we ended up with symmetric matrices with ones on the diagonal and correlation coefficients between -1 and 1 elsewhere. This transformation corresponds to making the signals zero-mean and with a unit variance before computing the covariance matrices. The correlation matrices were then used as input for the same estimation pipeline as the covariances: projection in the tangent space, vectorization, and regression. The chance level was represented by a dummy model, which always predicts the mean of the source domain (training data). Performances of all methods were evaluated with the coefficient of determination score or R 2 score, so the higher the score is, the better the model. Scores below 0 indicate predictions performing (arbitrarily) worse than the training-outcome-mean regressor that marks chancelevel prediction. Negative R2 values are in this sense worse than chance and point at bias due to systematic shifts or other distribution mismatches.

Figure 3 :

 3 Figure 3: Alignment method comparison across simulated dataset shift scenarios (R 2 score). Alignment methods (indicated by color) were evaluated on four different scenarios with an increasing shift. generated N = 300 matrices per domain to have data sets of the magnitude of real EEG datasets that would be considered as medium to large in terms of operational costs and curation effort. Error bars show standard deviations of the metric obtained with 50 random repetitions. The dashed vertical gray lines on (B) and (D) indicate the fixed parameter's value of the source set. Panel (A) displays the performance achieved when the target covariance matrices were created by multiplying the source mixing matrix with an SPD matrix: A (T ) = B α A (S) with B ∈ S + P . All methods that included re-centering the distributions on the same reference point performed well. (B) displays the performance achieved when the dispersion of covariances differs

Figure 4 :

 4 Figure 4: Impact of data alignment on age prediction across different tasks on the same subjects from Cam-CAN dataset (R 2 score). Alignment methods comparison for three different source-target tasks using 2000 repeat-bootstrap to select the subjects. Both domains contained the same subjects, only their task was different. Models are depicted along the y-axis, and standard boxplots represent their associated R 2 score. The dashed black lines represent chance-level performance. (A) Generalization of age prediction regression model from resting state to the passive task. Re-centering and the paired rotation correction led to an increased R 2 score with no obvious benefits for additional re-scaling. (B) The regression model was trained on resting-state data, and predictions were made on the recordings of the

  -CAN (MEG): different subjects. In this second benchmark, we focused on domain-shift differences between MEG

Figure 5 :

 5 Figure 5: Impact of data alignment on age prediction across different tasks for different subjects from Cam-CAN dataset (R 2 score). Alignment methods comparison for three different source-target tasks using 100 stratified Monte Carlo cross-validation (shuffle split) iterations to determine which subjects form the source and the target sets. We depict the models along the y-axis and represent the R 2 scores with standard boxplots. The dashed black lines represent chance-level performance. (A) The model was trained on the rest task, and predictions were made on the passive task recordings. When re-centering source and target distributions, prediction performance substantially improved, whereas re-scaling did change performance. (B) The target set was composed of recordings from the somatosensory task. The improvement of the re-centering step was smaller but still present. Re-scaling, still, did not lead to obvious improvements. (C) In the last

Figure 6 :

 6 Figure6: Impact of data alignment on age prediction across different EEG datasets (R 2 score). Data from the TUAB dataset were used as the source domain, and from the LEMON data as the target domain. We compare the alignment methods across 2000 bootstrap iterations on the source data (n=1385). The target set was always the same (n=213). The methods are represented along the y-axis, and we depict their associated R 2 scores with standard boxplots. The dashed black lines represent chance-level performance. (A) Results of alignment methods combined with the Riemannian approach of Equation11(as for all the results we have previously presented). Without

  , the filterbank-riemann pipeline trained on LEMON data only with a 10-fold cross-validation led R 2 scores of R 2 = 0.54 ± 0.13. Here, the training dataset only consisted of data from the TUAB dataset. The Riemannian re-center step made it possible to reach performance comparable to a model trained within the same dataset. With SPoC (B), re-centering led to a reduced standard deviation compared to no alignment. The highest R 2 scores were achieved when the re-scaling step was added to the alignment procedure and almost reached the performance of the Riemannian model trained on LEMON. Aligning the covariances distribution helped improve prediction performance even with a regression model like SPoC that does not leverage the geometry of the covariance matrices. This observation motivated an examination of how alignment affects the SPoC patterns, the inverse of the SPoC spatial filters W SPoC , and the resulting powers.

Figure 7 :

 7 Figure 7: Impact of alignment of different EEG datasets on their SPoC patterns and source powers. TUAB data were used as the source domain, and LEMON data as the target domain. Alignment refers to re-centering the source and the target distribution by whitening them respectively by their geometric mean. To obtain these figures, data were filtered in the alpha band. We included 19 channels (15 commons and 4 with similar locations on the scalp) in both datasets. (A) Topographic maps of the five first SPoC source patterns without alignment (first row) and target patterns without alignment (second row). The third row corresponds to the aligned source patterns adjusted with the target whitening inverse filter. These are the patterns applied to unaligned target data to obtain the target powers with alignment. The color map is normalized across each row. (B) Scatter plot of the target log powers as a function of the source log powers without and with alignment averaged across subjects. The dashed black line is the identity line. Alignment makes target and source log-powers more comparable.

  articulated in three alignment steps: re-center the geometric means, equalizing dispersions, and rotation correction. These alignment steps are evaluated in the regression context of generalizing age prediction across different domains. We investigated how dataset shifts can occur by analyzing a statistical generative model of M/EEG data. We presented simulated dataset shift scenarios based on this model for which alignment steps can effectively compensate the shift, plus a noise scenario to get a sense of how the methods would perform with real data. The simulation results showed that Procrustes paired is the most effective method in all scenarios. It was expected as it includes the three alignment steps and a rotation correction informed by the pairing of source and target subjects. We then designed M/EEG benchmarks with different domain definitions to determine the alignment methods' efficiency in those various settings. Coherently with the simulation results, Procrustes paired achieved the best performance, but since it cannot be applied in all situations, re-centering is the best option.We compared the alignment steps leveraging Riemannian geometry with a z-score method that transforms the covariance matrices into correlation matrices. This method systematically performed worse than all the others. Taking into consideration the geometry of the data space is essential. Among the three Riemannian alignment steps, recentering and the paired rotation correction of the source and target distributions help to improve the prediction performance in the M/EEG benchmarks. Re-centering and the paired rotation correction were shown to compensate for changes in the mixing matrices of the generative model, so we expected these steps to reduce the shift in benchmarks where the target population is not the same as the source population. In the first benchmark on Cam-CAN data, the source and target subjects were the same, but we still observed that re-centering and Procrustes paired led to better scores. On the other hand, equalizing their dispersions did not bring clear gains in performance in any benchmark. In the Cam-CAN benchmarks, the scores reached when the subjects are different in the domains are distinctly lower than with the same subjects: the shift is bigger (lower no alignment baseline) and harder to recover. For the EEG benchmark, we used the TUAB dataset as the source and the LEMON dataset as the target. Here, all recordings were done at rest but with different recording devices and in different populations. Re-centering the distributions in the EEG benchmark exceeded our expectations. Re-centering was sufficient to recover performance close to what is reached when training the Riemannian model on LEMON (Engemann

Table 2 :

 2 Summary of the simulation scenarios.

Table 3 :

 3 Definition of frequency bands.

Acknowledgements

This work was supported by the grants ANR-20-CHIA-0016, ANR-20-IADJ-0002 and ANR-20-THIA-0013 from Agence nationale de la recherche (ANR).

Numerical computation was enabled by the scientific

Python ecosystem: Matplotlib (Hunter, 2007), Scikitlearn (Pedregosa et al., 2011), Numpy (Harris et al., 2020), Scipy (Virtanen et al., 2020) and PyRiemann (Barachant et al., 2022), MNE (Gramfort et al., 2013).

Data and code Availability

All data used are publicly available. The Cam-CAN data repository is accessible via http://www.mrc-cbu. cam.ac.uk/datasets/camcan/, the TUH EEG data corpus via https://isip.piconepress.com/projects/tuh_ eeg/ and the Leipzig Mind-Brain-Body database via https: //fcon_1000.projects.nitrc.org/indi/retro/MPI_LEMON/ downloads/. The scripts and code for the alignment methods and the results presented in this paper are available on GitHub (https://github.com/apmellot/harmonizing_ aligning_meeg).

Author contributions

Conceptualization

×N is at most P . This implies that, computing the left singular vectors U ∈ R

×N that has at maximum P rows with non zero variances.

Thus, the other rows must be discarded using, for example, the class VarianceThreshold from the scikit-learn library [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF]. Otherwise, numerical issues can be encountered using functions like StandardScaler from the scikit-learn library. To prove these assertions, we begin by recalling the mixing model of C i with no noise in the next hypothesis.

Assumption 1 (Mixing model). We have a set of covariance matrices {C i } N i=1 that respect the following mixing model

with A ∈ R P ×P an invertible mixing matrix and with D i ∈ R P ×P diagonal with strictly positive elements.

This assumption induces that the Riemannian mean C of {C i } N i=1 defined in Equation ( 10) and the associated Riemannian logarithmic mappings {log(

have particular structures. These structures are computed in the following lemma.

Lemma 1. Knowing A and {D i } N i=1 , the Riemannian mean C has a closed form expression which is

Furthermore, the Riemannian logarithmic mappings of any

3)

The log function on the right-hand side of the equation applies the scalar logarithm on the diagonal elements. such that log(B) and log(EBE -1 ) exist, we have

These structures induce that the vectorized logarithm maps {z i } n i=1 only span a subspace of R P (P +1) 2

.

Proposition 1. The matrix Z = {z i , i = 1 . . . N } ∈ R P (P +1) 2

×N has a maximum rank of P .

Proof. We begin by defining the full vectorization counterpart of ( 11)

where vec is the operator that concatenates the columns of a given matrix. Then, by denoting s = vec(S), we get

where ⊗ is the Kronecker product. Denoting To conclude, the rows of Z are included in those of Z, hence rank(Z) ≤ P.