Deep Geodesic Canonical Correlation Analysis for Covariance-Based Neuroimaging Data
Authors: Ce Ju
Reinmar J Kobler
Liyao Tang
Cuntai Guan
Motoaki Kawanabe
Venue: ICLR 2024 spotlight (2024)
Keywords: Geometric Deep Learning
Self-Supervised Learning
Brain-Computer Interfaces
Neuroimaging
Neuroscience
Primary area: applications to neuroscience & cognitive science
TL;DR: A geometric deep learning-based approach to learn the SPD matrix-valued latent representation for paired covariance-based neuroimaging modalities under the self-supervised learning framework.

Abstract: In human neuroimaging, multi-modal imaging techniques are frequently combined to enhance our comprehension of whole-brain dynamics and improve diagnosis in clinical practice. Modalities like electroencephalography and functional magnetic resonance imaging provide distinct views to the brain dynamics due to diametral spatiotemporal sensitivities and underlying neurophysiological coupling mechanisms. These distinct views pose a considerable challenge to learning a shared representation space, especially when dealing with covariance-based data characterized by their geometric structure. To capitalize on the geometric structure, we introduce a measure called geodesic correlation which expands traditional correlation consistency to covariance-based data on the symmetric positive definite (SPD) manifold. This measure is derived from classical canonical correlation analysis and serves to evaluate the consistency of latent representations obtained from paired views. For multi-view, self-supervised learning where one or both latent views are SPD we propose an innovative geometric deep learning framework termed DeepGeoCCA. Its primary objective is to enhance the geodesic correlation of unlabeled, paired data, thereby generating novel representations while retaining the geometric structures. In simulations and experiments with multi-view and multi-modal human neuroimaging data, we find that DeepGeoCCA learns latent representations with high geodesic correlation for unseen data while retaining relevant information for downstream tasks.

Published as a conference paper at ICLR 2024
DEEP GEODESIC CANONICAL CORRELATION ANALY- SIS FOR COVARIANCE-BASER NEUROIMAGING DATA
Ce Ju Reinmar J. 2 Liyao Tang
Cuntai Guan' Motoaki Kawanabe
1 School of Computer Science and Engineering, Nanyang Technological University, Singapore
2 Advanced Telecommunications Research Institute International, Japan
3 RIKEN Artificial Intelligence Project, Japan
4 School of Computer Science, University of Sydney, Australia
Equal contribution
ABSTRACT
In human neuroimaging, multi-modal imaging techniques frequently com- bined to enhance our comprehension of whole-brain dynamics and improve diag- nosis in clinical practice. Modalities like electroencephalography and functional magnetic resonance imaging provide distinct views to the brain dynamics due diametral spatiotemporal sensitivities and underlying neurophysiological cou- pling mechanisms. These distinct views pose considerable challenge to learning a shared representation space, especially when dealing with covariance- based data characterized by their geometric structure. To capitalize on the geometric struc- ture, we introduce measure called geodesic correlation which expands traditional correlation consistency to covariance-based data the symmetric positive defi- nite (SPD) manifold. This measure is derived from classical canonical correlation analysis and serves to evaluate the consistency of latent representations obtained from paired views. For multi-view, self-supervised learning where one or both la- tent views are SPD propose innovative geometric deep learning framework termed DeepGeoCCA, Its primary objective is to enhance the geodesic correlation of unlabeled, paired data, thereby generating novel representations while retain- ing the geometric structures. In simulations and experiments with multi-view and multi-modal human neuroimaging data, we find that DeepGeoCCA learns latent representations with high geodesic correlation for unseen data while retaining rel- evant information for downstream tasks.
INTRODUCTION
In human ncuroimaging, is common to use multi-modal imaging techniques to improve our un- derstanding of the dynamic relationships among brain networks and their alterations in pathol-
ogy (Debener et al., [006]] Ou et al., 2010: Huster et al.] 2012] /aadahhaani & Wirsich] (2000). For instance, combining non-invasive electreencephalography (EEG) measurements with functional magnetic resonance imaging (IMRI) can offer a high-resolution view of neurophysiological events relevant to treat mental disorders like epilepsy (Moeller al.,
2009). However, achieving sense
of consistency between EEG and fMRI frequently presents significant challenges due to a range of complications stemming from their disparities (Daunizeau al. 2009 Warbrick 2022).
A classical statistical learning approach to identify paired subspaces for consistency through corre- lation Canonical Correlation Analysis (CCA) (Hotelling 1992). Its versatility is evident in a range of scenarios even when the sample size is limited relative to data dimensionality or when data di- mensionality exceeds human interpretability (Uurtio et al. (0017). Beyond classical CCA, extensions coupling proved useful in uni-modal Dahne et al., (214)) and multi modil neuroimaging (Dähne et al. 2015) Deligianni et al. 2014). For instance, Deligianni et al. (001 2016) used inverse modeling to align paired, resting-statt EEG and fMRI signals spatially. After alignment they used a CCA variant to link the covariance-batt EEG and fMRI data, and assessed the goodness-of-ft with geonetry-aware evaluation criterion.
to study
coupling and amplitude-amplinne

---

Published as a conference paper at ICLR 2024
expa(ff) expo (hi)
Dimension Adaptation
Tangent Space Projection
W
(te)
S
EEG
exp. (tar)
Transformation
Zi
S
T,ST S S
fMRI
S
Figure I: Diagram of DeepGeoCCA Framework. We consider paired views of multivariate time- series data (e-g.. simultaneously recorded EEG and fMRI data). Individual views are transformed (e.g., covariance matrices) to form N paired views on SPD manifolds with different dimensions and rıg, i.c., {(S,, | E 4 St E 8T2 +++ A Dimension Adaptation module transforms the views so that they lie on + a cone manifold R"x*, On this cone, a Transformation module learns to center the data around the identify matrix and project them onto geodesics expe (Iv) and exp. (tw) to form latent SPD representations {(Z, Z() The aim of this problem is to maximize geodesic correlation which measures correlation along the tangent speeds v and w, after the Tangent Space Projection at e.
When dealing with paired covariance- based neuroimaging data, using conventional CCA methods may lead to a loss of certain mathematical structures that are important for practical results, such as properties of symmetry and positive definiteness (SPD). Modeling techniques based on Rieman- nian geometry have become the primary approach to address such constraints in covariance based
neuroimaging data (Fillard et al., 2007: Barachant et al., 2011 Qiu et al.) 2015: You & Park [021]]
Dan et al. 2022). In the case of CCA, extensions generalized the concept of subspace projection in
Euclidean vector space non-linear projections on geodesic submanifolds for general Riemannian
manifolds (Kim et al. 2014) and SPD manifolds (Kim et al. 2014 Fallah & Yang 2020).
In the context of multi-view deep representation learning (Wang al. 2015), CCA inspired fam-
ily of self-supervised learning (SSL) approaches (Caron et 2020 Zbontar et al. 2021 Ermolov
202T Bardes al., 2022), as categorized in Balestriero al. (2023). Methods, such Deep- CCA (Andrew et al. 2013) and VICReg (Bardes et al.) (022), proposed to learn the correlation relationship between two latent representations by analyzing their covariance and cross-covariance structure. However, for covariance-based neuroimaging data with specific geometric characteris- ties, these methods fail maintain these particular geometric properties when learning correlation relationships.
al.
Here, we propose a geometric deep learning framework, named Deep Geodesic Canonical Corre- lation Analysis, or DeepGeoCCA. It belongs to the SSL-based deep representation learning family (Chen et al. 2020) as well as geometric modeling techniques to learn latent SPD representations with high correlation consistency among paired views (Figure[). I]. To characterize this correlation consistency, we present a new correlation measure called geodesic correlation, which is applicd to SPD) manifolds equipped with the affine-invariant Riemannian metrie (AIRM) (Pennec al. 2006). Intuitively, geodesic correlation essentially measures correlation of paired data along geodesics pass- ing through the identity matrix e. To maximize this measure, DeepGeoCCA introduces novel geodesic loss function and SPD matrix-valued neural network models to learn latent SPD represen- tations. Specifically, we devise relaxation method, called e- geodesic constraint. that enables deep neural networks on SPD manifolds (Huang & Gool [017]] Ju & Guan 2022] 2023 Kobler et al.) 2022ah to map covariance-basee neuroimaging data onto or near geodesics within double-cone region around e.
Using theoretical analysis and experiments with simulations and real data, we systematically study our framework, detailed in section and show that DeepGeoCCA:
generates latent representations which maintain the properties of symmetry and positive definiteness.
supports a controllable projection region. Specifically, we prove that adjusting in our relaxation approach control the deviation of the mapped points from the geodesics.
2

---

Published as a conference paper at ICLR 2024
allows unit speeds for predefined geodesics, in the sense that they can be pre- pre-defined by characteristics of downstream tasks or driven by data.
effectively learns latent representations exhibiting significant geodesic correlation for pre- viously unseen data while preserving pertinent information (section
2 PRELIMINARY
We first briefly recap standard CCA, a deep learning-based CCA extension, and a CCA method on general Riemannian manifolds. Please refer to Appendix A for introducing SPD manifolds equipped with the affine-invariant Riemannian metric, Appendix B for Fréchet means and vari- ances, Appendix foror geodesics. We extensively examine and contrast the following CCA variants in Appendix
2.1 NOTATIONS AND CONVENTIONS
0, Vr R" /0) to be the space of real n X n symmetric and positive definite matrices. The set S"h := {S E R"*" ST) } is the space of n X n real symmetrie matrices. We always use e to represent the n x identity matrix In E S (ST (ST ) to refer to the Riemannian manifold formed by equipping + with the AIRM. 11-11. is the Riemannian norm at identity The notation J denotes the Frobenius norm. (-. ) represents the -n-irner product. exp(.) and log(-) are the matrix exponential and the matrix logarithm. Tr(-) is the trace of a matrix. The operator vec(-) represents the vectorization of a matrix.
Sx In the following context, employ the set S4+ := IS : S= and TSr >
2.2 CANONICAL CORRELATION ANALYSIS (HOTELLING 1992)
Consider two paired matrices of data distributions X E N and X RN9, where N represents the number of observations, while p and y denote the dimensions of two variables. Paired means that each 7th row in both X and X are coupled (e.g., recorded at the same time). We assume that the columns in both X and X are zero-mean. The aim of CCA is to find two linear transformations, de- noted as w E Rº and w E R2, in such a way that the cosine of the angle between the projected vec- tors 2t and gl onto the unit ball is maximized, as described below, cos θ max,t, ERR 21, ,
where z! (X / will E R and gl == (X w 1| E R
2.3 DEEP CANONICAL CORRELATION ANALYSIS (ANDREW ET AL., 2013
Deep canonical correlation analysis (DeepCCA) employs neural networks to learn latent represen- tations with the goal maximize the CCA objective between latent representations of two modali- ties, as follows, max,,0 ā corr (<p(1;; 0), ((;; 0)), where p((:00) and (f.;; 0) are neural networks with parameters and 0, respectively. Inspired by CCA, Andrew et al. (2033) proposed to compute and maximize Σー Σ ppp"" Σ=1/2 1/2 for each batch, where Enp and Σg represent within-
modality covariances of latent representation, and Enp represents cross-modality covariances of latent representations. Their preprocessing procedure includes centering and regularization.
2.4
RIEMANNIAN CANONICAL CORRELATION ANALYSIS (IM ET AL. 2014)
Given a set of paired siiifolddvalued data {ax) and { (M,g). Riemannian canonical correlation analysis (RieCCA) is the correlation between orthogonal projections of {ar): and J onto two geodesic submanifolds exp,, (U), exp,, C M, where two open neighborhoods U T and Ü Ty.M are spanned by tangent vectors Fréchet means μ and jú, respectively, and two orthogonal projections r/. == minggexn,,(U) d2 (z;, x) and कt == minggexp,, (0) d (T; F) the geodesic submanifolds, for all i N). where d..-.)) is the Riemannian distance with respect to the Riemannian metric g.
is the symbol that is typically utilized the identity element in Lie group and Lie algebra contexts.
3

---

Published as a conference paper at ICLR 2024
3 DEEP GEODESIC CCA (DEEPGEOCCA) FRAMEWORK
Consider paired views of novariance-based neuroimaging data, whether from the same or different modalities, using established criteria. For example, EEG and IMRI views are paired based on their recording timestamp. In the proposed framework, the need for labels is eliminated, offering conve- nience for neuroimaging data derived from complex brain dynamics. Deep neural networks on SPD manifolds are utilized to learn latent SPD matrix-valued representations by maximizing geodesic correlation. As such, DeepGeoCCA can be considered as a specific instance within the SSL family, designed particularly for handling paired SPD matrix-valued representations.
Without loss of generality, suppose N paired SPD matrix- representations { and (5:) are on ($"+,gAARM, ALRMM, ) each with Fréchet mean e, and two arc-Iength geodesics expe(tv) and exp, (tw) S" - S with time t E R and tangent unit speed v, w E T.S _ S", respectively. Geodesic correlation is defined follows with an illustration in Figure 2
Definition (Geodesic Correlation). Geodesic correlation with respect to the given langent speeds is defined as the correlation between their Riemannian logarithm of the orthogonal projections on the langent space at the identity matrix E S4 as follows,
e and w on
Σ 9c
AIRM (loge (S)), loge (S:))
1.
VE N log. (s))l . VE log. (S!)
where S and S represent the orthogonal projections on the given geodesics.
Note that all the orthogonal projections are unique, which is a result of the affine-invarian Riemannian metric that leads the SPD space to be a Hadamard manifold (Bacák) 2014 Pennec 2020). Hence, the geodesic correlation is well- defined. Please refer to Appendix C for orthogonal projec- tion and geodesic constraints, Appendix E for additional mo- tivations and main difference with Riemannian Correlation. Besides, incorporating geodesic correlation and introducing geodesic constraints can be reformulated as constrained non- linear optimization problem, and its Lagrangian duality nat- urally serves as a loss function for the neural network-based solutions, as explained in Appendix F
exp, (te) exp, (tre
T.S2 S
Figure 2: Geodesic Correlation.
3.1 NEURAL NETWORK- -ASSE SOLUTION
As outlined in Figure] transform paired views of covariance-based neuroimaging data into latent SPD matrix-valued representations centered around the identity matrix e. To that end, we em- ploy neural network architectures SPD manifolds. These specialized architectures can maintain the properties of symmetry and positive definiteness throughout their computations.
The proposed framework comprises three modules. The first module, known as Dimension Adapta- tion network, is responsible for converting input data varying dimensions to a consistent dimen- sion. The second module, referred to Transformation network, learns non-linear transformations to center the data around the identity matrix e and project them close to geodesics passing through Lastly, Tangent Space Projection maps the observations to the tangent space to form inputs for our proposed loss terms. In particular, to enable computations to be performed at e, we employ Riemannian Batch Normalization on (S g. g ATRM) (Brooks al. 2019) with specially tuned step momentum parameter to parallel transport Fréchet means of SPD data distribution to (Kobler et al. 2022ah. For a more detailed overview of the employed SPD layers, please refer to Appendix
3.2 RELAXED ORTHOGONAL PROJECTION
In this section, we propose relaxation method for the geodesic constraints of geodesic correlation, which is particularly suited for the neural network- based solution. In intuitive terms, this relaxation assesses whether the projection Z = St tangent space T.S after being normalized, falls
4

---

Published as a conference paper at ICLR 2024
within the double cone centered at identity e, with the unit geodesie speed the symmetric line, by calculating the cosine value between and u, defined as follows.
9- g.TRM (loge (Z). u)
cos ())
||og [Z)le ||u].
relaxation parameter e E [0, 1] is introduced to control how much the projections can deviate from the geodesic by S [cos(E)]. We refer to this new constraint the geodesiic constraint. Fur- thermore, we provide an estimate of the discrepancies between projections precisely on the geodesic and those within the tubular neighborhood of the geodesic exp, (tu) as follows,
Theorem. Suppose { are a set of latent representations on (5 tubular neighborhond of are-length geodesic exp, (fu) with unit speed u is defined (expe (tu: + w)le (T,S++) ||e|le I 6, t E RR, where (T.S4+)+ represents the orthogonal complement of TeS+.. Then, the width of this rubular neighborhood is upper bounded by the following inequality,
and a 8-width
62
85
max [)...NN
ge gALRa (loge( (Zi),
E2
where relaxation parameter & E (0, regulates the deviation of the projections from the geodesics by 8 < [cos()) In particular, when e + 1, we have 6 + O.
For comprehensive proof, please consult Appendix H
3.3 Loss FUNCTION
Within the proposed framework, SPD neural nets are trained to minimize an overall loss function L comprising correlation loss term Lp. relaxed e geodesic constraints Lg, and variance-pressrvvnng loss L
= a1 Lp + 02 L. + ag Cg. with term coefficients 01 02, and 03 2 0. (Tips of choice for term coefficients refers to Remark
(3)
Specifically, the correlation loss Lp is utilized maximize geometric correlation among paired views as follows,
EN Tr (log. (7a) log, (2:)
EE Tr(log (1)) EM= Tr( (lo)
where Z, and Zi denote the latent SPD matrix-valuee representations for two paired views, respec- tively, and symmetric matrices log, (2,) and log, (Z:) are their projections on the tangent space at e; The relaxed E -geddesic constraints L. are utilized to penalize projections outside the double-cone regions as follows.
Tr log, (Zx) bb)
N
Tr log, (2,) ww
Le = min (e, [log. (z)) 1l. ||0l|."
Σ i=1
min (c, 1
logg (Z, ||. ||wlle
i=1
where geodesic speeds v, w E S" are predefined or learnable parameters, and E is a preset relaxation parameter; Furthermore, variance- preeeevvngeeeegnggggge loss L. proposed in (Bardes et al. (022) is utilized to mitigate potential dimensional collapse, which frequently encountered in the SSI framework (Jing
al. 2022) as follows,
C, Σ max {0, -
10.1-V
--a--VarrSs Var(S(Z))) + +
Σ max (0, Var(S(())) + +c),
where d is the dimension of vee(log, (Z)), S(Z)² is the vector composed of values the gth di- mension of (vec(log, (Zi)), vec(log, (2ologg..()))) (Zx))) E Rxxxd and is scalar to prevent numerical instabilities, and
10 in practice.
Remark 1. Choice of Initial Geodesic Speeds: The selection of the initial geodesic speeds for two covariance based modalities can be influenced by the specific requirements of downstream tasks or solely dictated by the data. For example, a choice can be the largest canonical correlation obtained from the conventional CCA vec(u), vec(w) CCA([vec(S)) (. (vec(S.) (.)
5

---

Published as a conference paper at ICLR 2024
Remark 2. Choice of Loss Term Coefficients: Depending the specific scenario and architecture, the coefficients (01, Q22, az)) can be determined to give more importance to geodesic correlation or the constraints. By default. we set a 1,02 0.25, 03 to give more weight Er If specific architecture tends to dimensional collapse, we enable the variance-preserse
loss (og > 0);
Remark 3. Extending Overall Loss function S Since all loss terms in the overall Loss (3) operate with tangent space projections loge (Za), loge (Z.) E T corresponding to the space of symmetric matrices S", the proposed loss also be computed for any neural network generating latent representations in S For example, the output x E [2((++))22 of any standard neural net-
work layer can be transformed to 8 via applying upper where the norm-preserving map upper vectorizes a matrix's 's elements along the upper triangular part.
EXPERIMENTS
We conducted simulations and two experiments with public human datasets to empirically evaluate our proposed framework in challenging neuroimaging tasks. The nature of these tasks poses sev- eral unique challenges, including relatively small datasets, multivariate non-stationary signals with complex dynamics, distribution shifts (severe across subjects; also indicated in Figure low signal- to-noise ratios, and non-linear coupling mechanisms between EEG and fMRL. See Appendix for more details.
4.1 SIMULATIONS
To study the effect of the E e-geodesic constraint and demonstrate that the proposed framework can recover the activity of a latent, shared source, generated paired SPD matrix-valued data. That is, we generated a shared noise source with fixed, view-specific projections to S and paired it with a view-specific noise source. Parallel transport and exponential mapping were used to project the data to S with a random, view-specific Fréchet mean. Figure 3 summarizes the generative process and the considered architecture. For both modalities, we used SPDNet (Huang & Van Gool (017) as dimension adaptation and transformation module (Figure 3b). As tangent space projection module, we combine SPD momentum batch norm (SPDMBN) layer (Kobler et al., 2022a) and LOG layer. For thorough introductions of neural networks on SPD manifolds, please refer to Appendix G
To quantify performance, we report conventional Euclidean vector space scores. For geometric approaches, we converted the latent SPD matrix-valued representations to Euclidean vector space features via tangent space projection. Conventional CCA was then applied to the training data to estimate one component. We report Pearson's 's correlation coefficients and coefficients of determina- tion R22 on the held out test data.
Figure 3 summarizes cross validation (CV) results of our framework for various e compared to sev- eral baseline methods. The baseline methods were either fitted the generated SPD data (CCA and RicCCA) or used as alternative loss terms (DeepCCA, VICReg) for the considered architccture. For E (no geodesic constraint), DeepGeoCCA shared the highest correlation score with DeepCCA among the deep learning methods (left panel) at the cost of small explained variance (i.e., small R²) in the latent representations (right panel). This result can be expected because both objectives only aim to maximize correlation in latent space. For able reduction in correlation for minimizing deviations from the geodesic constraint. Moreover, for some CV splits, the network got stuck in poor local minima (noticeable lower correlation and R²). We found that = 0.75 resulted in good trade-off between both objectives while avoiding poor local minima. Scatterplots in Figure3p visualize the learned representations (s 0.75) for rep- resentative CV split. Comparing the correlation scores of CCA and RieCCA, clearly demonstrate the importance of utilizing the mathematical structure of SPD matrices. Overall, RieCCA obtained the highest correlation score. This is likely because RieCCA had access to the entire training data to find the optimal solution for this toy problem, while the deep learning approaches were fitted to randomly sampled mini-batches
the trained network traded notice-
2 The coefficient of determination H² is used to quantify how much variance in the latent representations can be explained by the first CCA component, i.c., R² = 1 4 data distributed along line or geodesic for geometric approaches.

---

Published as a conference paper at ICLR 2024
generative process and generated SPD data
b latent SPD representations
C quantitative results correlation
Z
R2
inpul CCA
shared source
Randern Bios
1.00
R ReOGA
modality
EXP
Sp SPDNOR
LOG
Geodesice Correlation
Intent UneeODA VICRlg
0.8
noise
proposed
c11
(-0.72
Próchet moan #
0.90
COO.
modality 2
Randem Bles
S
noise
sspo SENON EEDNGR
EXP
LOG
0.4
shared source
0.60
Figure 3: Simulations on S a. Generative process and visualization of generated observations (N 500) with shared latent source. b, Visualization of the latent representations learned with DeepGeoCCA (e = 0.75). c, Boxplots summarize the simulation results for test data (10-fold CV). DeepGeoCCA obtain high geodesic correlation (left panel; higher is better) and the same time approx. constrain the representations to geodesics (right panel; coefficient of determination R2; mean across modalities; higher is better).
4.2 SIMULTANEOUS EEG-FMRI DATA
The aim of this experiment is to exhibit the feasibility of DeepGeoCCA, Since EEG and fMRI provide distinet views to brain dynamics due to their diametral spatiotemporal sensitivities and un- derlying neurophysiological coupling mechanisms, the proposed approach is fundamentally differ- ent from previous data-driven approaches which used extensive
-nwaalnn-kwwwddee to extract low- dimensional representations (Deligianni et al.] 2014 Dähne et al. [2015] Keynan et al. 2019 Wirsich
et al. 2020 Abreu et al., 2021 Ogawa et al. 2021 Philiastides et al. 2021) or tried to transcode
them (Liu al. (2022). Using simultaneous EEG and fMRI measurements,
aim learning la-
tent representations where EEG and IMRI dynamics are highly correlated, and remain so for held data. A publicly available simultaneous EEG-IMRI dataset (van der Meer et al.! 2016 used to demonstrate feasibility. It contains data of 8 subjects while they were resting in two different conditions (eyes closed or open). Both conditions lead changes in whole-brain dynamics, which are captured by EEG (Barry et al.) 2007) and fMRI (Costumero et al. 2020), rendering the dataset an ideal candidate to evaluate our proposed framework against baseline methods. For detailed description of the dataset, EEG and fMRI preprocessing and implementation details, please refer to Appendix[K]
4 summarizes our experimental approach. We used t-SNE to visualize the preprocessed EEG and (MRI data ((S,, S as well as the learned representations {(z, Z for representa- tive CV split. For EEG data, tested two architectures, both extracting latent representations in $10.. The first one was similar to the architecture used in the simulations. The second one (TSM- Net) is conceptually similar except that it applies convolutional feature extractors EEG time-series data before covariance pooling (Kobler et al.) 2022a Wilson et al.) 2022). To accommodate large distribution shifts (e.g., g.. in top-left (-SNE plot Figure 4h across subjects and runs often encountered
in EEG (Jayaram et al. 2016 Rodrigues et al. 2019 Wu et al., 2020: Défossez et al. 2023). we 2022a). To l'aiilitate fair comparison with RieCCA and CCA, we re- re-centered the inputs for these methods per run according to (Zanini et al. 2018). 3. For fMRI data, we used a standard
used subject and run-specific batch normalization layers (Kobler
cewasformer-enccodee layer (Vaswani al.) proved useful in learning latent representations from large-scale, diverse fMRI datasets (Thomas et al. 2022p. This approach considered each fMRT volume inside sliding window as token. At the end of each window, concatenated Iearnable summary (SUM) token. Finally, linear layer projected the transformed SUM token into a symmetric matrix which we used according to Remark as log, (Z;) in the overall loss
Table Ssummarizes the results for a 10-fold cross- validation scenario. As in section 4.1 we use Pear- son's 's correlation to assess the consistency between the learned EEG and fMRI representations and R² to quantify the geodesic constraints. Utilizing TSMNet in DeepGeoCCA obtained significantly higher correlation than all other candidates. Learning the considered architectures with alternative
7

---

Published as a conference paper at ICLR 2024
subjoct run run
subject resting eyes dosed
EEG
74 F S"
SP
eyee open
(saunnox Ir
Cav Pool
EEG
ISNE dimz
Goedesic Correlation
1SNE dim
CR
S. R
dim1
time
wincowg
fMRI
ISNE dimz
vokuma
&
poc. ww..embeaa
Transformn erzzddr
1MRI
ISIEE
S E k2
18NEE dim1
Z4 E 5
Figure 4: Diagram of DeepGeoCCA for Simultancous EEG-IMRI Data. Using simultaneously recorded observations (S,, S.). we aim to learn a latent space where brain dynamics (74, Z,). shared between EEG and IMRI covary with high-congruence. After preprocessing and extracting sliding windows, we use established decoder models to convert IMRI activity (Thomas et al. (022) and oscillatory EEG activity (Kobler et al., 2022a) to latent representations. Visualizations (t-SNE, perplexity=30) summarize input and latent data distributions after fitting.
Table Results for Simultaneous EEG-IMRI Data. Test-set model performance (higher is better) is evaluated with correlation and R² metrics across 10-fold CV with stratification across subjects and runs. Exhaustive permutation t-tests (df=9, tests with t-max adjustment) were used identify significant differences between TSMNet=DDeepGeCCCA and baseline methods.
IMRI model
EEG model
Method
Correlation 1
R²
mean (std)
t-val (p-val)
mean (std)
t-val (p-val)
CCA
0.02 (0.15)
-9.9 (0.002)
0,07 (0.08)
-19.8 (0.002)
RieCCA
0.38 (0.10)
-5.4 (0.002)
0.01 (0.01)
-88.4 (0.002)
DeepCCA
0.42 (0.11)
-3.5 (0.031)
0.30 (0.07)
-12.3 (0.002)
SPDNet
VICReg
0.29 (0.19)
-4.0 (0.018)
0.28 (0.10)
-92 (0.002)
Transformer
DeepGeoCCA
0.44 (0.11)
-3.0 (0.049)
0.58 (0.02)
3.8 (0.023)
Encoder
DeepCCA
0.29 (0.10)
(0.002)
0.18 (0.09)
-15.2 (0.002)
TSMNet
VICReg
0.38 (0.10)
(0.002)
0.24 (0.03)
40.5 (0.002)
DeepGeoCCA
0.58 (0.08)
0.55 (0.02)
loss terms (VICReg, DeepCCA) resulted in obvious performance drops for either or both EEG chitectures. Moreover, compared to the performance of RieCCA, the results clearly demonstrate the utility of DeepGeoCCA extract shared, latent brain dynamics whose coupling generalizes to unseen data. Note that cross subject generalization results, listed in Table of Appendix[[] show similar effects, confirming the robustness DeepGeoCCA. Taking a look the EEG architectures. the results suggest advantage of combining learnable convolutional layers with latent covariance pooling over broad-band EEG spatial covariance matrices as inputs.
4.3 MULTI-VIEW EEG DATA
The aim of this experiment is to evaluate the performance of DeepGeoCCA in a downstream task, namely baG-based motor imagery classification. During planning, executing but also imagining limb movements, EEG captures spatio-spectrally localized fluctuations in sensorimotor rhythms (Pfurtscheller Da Silva 1999). These fluctuations give rise to identifiable patterns that consis- tently enable classification of EEGs.
Motivated by the emergence of mobile, few-channel EEG devices, we generate pairs two-channel EEG and multi-channel (full sensorimotor coverage) EEG views. We use DeepGeoCCA for model pre-training with the goal to improve downstream classification performance of the 2-channel EEG views. Here, kept geodesic speeds component directions of the EEG spatial covariance matrices of the full-channel EEG views. For detailed settings of this experiment, please refer to Appendix
fixed and initialized them with tangent space principal
8

---

Published as a conference paper at ICLR 2024
Using DeepGeoCCA to pre-train Tensor-CSPNet (Ju & Guan 2022) resulted in statistically signif- icant group-level improvements in both considered datasets (Table and Supplementary Table (5) compared to the lower bound (i.e.. two EEG channels without SSL pre-training). Although our pro- posed loss consistently yielded significant improvements, the effect remained small across scenarios compared other SSL. pre-training losses. Nonetheless, the results clearly indicate that our frame- work for pre-training SPD neural nets like Tensor-CSPNG with paired views can yield performance increases downstream tasks. For detailed results the subject-level. refer Figures and 10 in [
Table 2: Downstream Results for Multi-view EEG Data. Average (std across subjects in brack- ets) classification accuracy for the motor imagery task (higher better). We considered two public datasets, and two scenarios, namely, within-session 10-fold CV (sessions Si and S2) and a hold-out scenario (train S test S2). For the KU and BNC1201500) datasets the all channel EEG views com- prise 20 and 13 channels above the sensorimotor cortex. Pensor-CSPNe is utilized as the network architecture in DeepGeoCCA.
Channel
SSL pre- pre-training
KU (54 subjects, classes)
BNC12015001 (12 subjects, classes)
81
S2
S1 + 522
S1
S2
S1
82
Two
x
62.97 (15.72)
65.85 (16.03)
62.89 (12.18)
74.12 (18.20)
75.62 (14.22)
75,00 (14.10)
DeepCCA
66.24 (15.73)
68.11 (15.67)
65.20 (13.24)
76.46 (16.83)
78.33 (11.33)
76.58 (13:23)
VICReg
66.12 (15.80)
67.35 (16.80)
65.20 (13.51)
77:21 (15.96)
77.25 (13.46)
77.00 (13.97)
DecpGeoCCA
66.26 66.26 (17.05)
68.27 (17.23)
68.93 (13.68)
77.58 (16.26)
78.46 (12.82)
77.83 (13.11)
All
X
70.56 (16.26)
72.86 (15.86)
68.85 (14.18)
84.33 (10.10)
86.46 (14.19)
85.17 (11.1)
5 DISCUSSION
Motivated by wide-spread applications of covariance- neuroimaging data, capture the cor- relation consistency of paired views on SPD manifolds using a novel measure called geodesic cor- relation. This novel correlation is generalized concept of correlation for SPD matrix-valued data. is well suited for deep neural networks on SPD manifolds compared to conventional correlation measures. To maximize geodesic correlation in latent space between paired views of SPD matrix- valued representations, we propose a novel loss function using the relaxed geodesic constraints, and present novel geometric deep learning-based SSL framework, referred to as DeepGeoCCA, Deep- GeoCCA is more generally applicable than classical non-deep learning approaches, because, in a general sense, deep learning is capable of handling datasets with intricate high-dimensional data distributions. and align latent representations of neural data collected from diverse population. In particular, our framework respects the geometric structure of SPD matrix-valued data. In simula- tions and experiments with EEG and fMRI measurements, DeepGeoCCA can learn representations whose dynamics generalize favorably to held-out data, while preserving task-relevant information.
Overall, this work can be regarded as proof-of-concep that geometric deep learning- dearning-based SSL can be used to extract meaningful representations from simultaneous EEG-IMRI measurements solely from paired observations. Although experiments with EEG-IMRI data is limited to a small dataset with task-related effects in either modality, our work is the first completely data-driven SSL frame- work in this context. In practice, one may expect further gains by utilizing large-scale datasets and uni-moda SSL pre- pre-training objectives proposed in carlier studics, for example Thomas et al. (2022);
Banville et al. (2020).
Besides, for applications beyond the considered tasks, it is necessary to identify interesting factors of variation to generate paired views. In clinical setting, task could be disease diagnosis. Naturally, views associated with the same subject could be considered as paired. For example, our framework could contribute relating functional, derived from fMRI, and anatomical, derived from diffusion tensor imaging, connectomes with the goal to identify latent, potentially clinically relevant, modes of variation across subjects.
9

---

Published as a conference paper at ICLR 2024
ACKNOWLEDGMENTS
This work was supported by Innovative Science and Technology Initiative for Security Grant Num- ber JPJ004596, ATLA and KAKENHI (Grants-in-Ains for Scientific Research) under Grant Numbers 21H03516, 216.,21111 21K12055, JSPS, Japan.
This work was supported under the RIE2020 Industry Alignment Fund Fund-Industry Collaboration Projects (IAF-ICP) Funding Initiative, as well as cash and in-kind contributions from industry part- ner(s); This work was also supported by the RIE2020 AME Programmatic Fund, Singapore (No. A20G8b0102).
REFERENCES
Pierre Ablin, Gean-Francois Cardoso, and Alexandre Gramfort. Faster Independent Component Analysis by Preconditioning With Hessian Approximations. IEEE Transactions on Signal Pro- cessing, 66(15):4000-4449, 2018,
doi: 10. 10.110999PP.0088844400000
Rodolfo Abreu, João Jorge, Alberto Leal, Thomas Koenig, and Patrícia Figueiredo. EEG Mi- crostates Predict Concurrent fMRI Dynamic Functional Connectivity States. Brain Topography, 34(1):41-55, 2021. doi: 10.1007/s 10001000000000000--1111111
P-A Absil, Robert Mahony, and Rodolphe Sepulchre. Optimization algorithms on matrix manifolds. Princeton University Press, 2009.
Philip J. Allen, Oliver Josephs, and Robert Turner. A Method for Removing Imaging Artifact from Continuous EEG Recorded during Functional MRL. Neurolmage, 12(2):230-239, 2000. doi: 10. (0000///mm..0000000..00
Galen Andrew, Raman Arora, Jeff Bilmes, and Karen Livescu. Deep canonical correlation analysis. In International conference on machine learning, pp. 1247-1255, PMLR, 2013.
Miroslav Bacák. Convex analysis and optimization in Hadamard spaces, volume 22. Walter de Gruyter GmbH & Co KG, 2014.
Randall Balestriero, Mark Ibrahim, Vlad Sohal, Ari Morcos, Shashank Shekhar, Tom Goldstein, Florian Bordes, Adrien Bardes, Gregoire Mialon, Yuandong Tian, et al. A cookbook of self- supervised learning, arXiv preprint arXiv:2304. 12210, 2023.
Hubert Banville, Omar Chehab, Aapo Hyvärinen, DeeianAlexander Engemann, and Alexandre Gramfort. Uncovering the structure of clinical EEG signals with self-supervised learning. Journal of Neural Engineering, 2020. doi: 1088/1741-2552/abbe
Alexandre Barachant, Stéphane Bonnet, Marco Congedo, and Christian Jutten. Multiclass brain- computer interface classification by riemannian geometry. IEEE Transactions Biomedical Engineering, 5944)9920-928, 2011.
Adrien Bardes, Jean Ponce, and Yann LeCun. VICReg: Varianceiinvarancc--vvaii regulariza- tion for self-supervised learning. In International Conference on Learning Representations, 2022.
URL /tR//tttp////perevieee wet.iewww.eeeeeee..eee..ee.eeeee forum 1dt.mmmmmm0xxmmmmm0000000
Robert J. Barry, Adam R. Clarke, Stuart J. Johnstone, Christopher A. Magee, and Jacqueline A. Rushby, EEG differences between eyes-closed and eyes-open resting conditions. Clinical Neuro- physiology, 118(12):275-2777) 2007. doi: 10. h001//..iinhh.200..00
Gary Becigneul and Octavian-Eugen Ganea. Riemannian Adaptive Optimization Methods. In In- lernational Conference Learning Representations, 2019, URL https::/ //tt:::///ppereiiew net paet/forune htt.:/oooumn11x=leeiii035.
Daniel Brooks, Olivier Schwander, Frederic Barbaresco, Jean- Yves Schneider, and Matthieu Cord. Riemannian batch normalization for spd neural networks. Advances Neural Information Pro- cessing Systems, 32, 2019.
10

---

Published as a conference paper at ICLR 2024
Madeleine Bullock, Graeme D. Jackson, and David F. Abbott. Artifact Reduction in Simultaneous EEG-IMRI: A Systematic Review of Methods and Contemporary Usage. Frontiers in Neurology, 12:622719, 2021.
doi: 10.3389/feeu:22211.22117
Mathilde Caron, Ishan Misra, Julien Mairal, Priya Goyal, Piotr Bojanowski, and Armand Joulin. Unsupervised learning of visual features by contrasting cluster assignments. Advances in neural information processing systems, 33:9912-999 2020.
Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey Hinton. A simple framework for contrastive learning of visual representations. In International conference machine learning, pp. 1597-1607. PMLR, 2020.
Víctor Costumero, Elisenda Bueichekú, Jesús Adrián- Ventura, and César Ávila. Opening or clos- ing eyes rest modulates the functional connectivity of V1 with default and salience networks. Scientific Reports, 10(1):9137, 2020. doi: 10.1038/s41198- 11188000--6610---- 66100-y.
Frederick Emory Croxton. Applied general statistics. In Applied general statistics, pp. 754-754. 1967.
Kamalaker Dadi, Gaël Varoquaux, Antonia Machlouzarides-Shallee Krzysztof Gorgolewski, Demian Wassermann, Bertrand Thirion, and Arthur Mensch. Fine-grain atlases of functional modes for fMRI analysis. Neurolmage, 221:117126,, 2020. 117126.
doi: nour01/..
Tingting Dan, Zhuobin Huang, Hongmin Cai, Rohert G Lyday, Paul J Laurienti, and Guorong Wu. Uncovering shape signatures of resting-state functional connectivity by geometric deep learning on riemannian manifold. Human Brain Mapping, 43(13):3970-3988 2022.
Jean Daunizeau, Helmut Laufs, and Karl Friston. EEG-IMRI Information Fusion: Biophysics and Data Analysis. In Christoph Mulert and Louis Lemieux (eds.), EEG fMRI, pp. 511-526. Springer Berlin Heidelberg, 2009, ISBN 978-3-540-87918-3 978-3-540-8791990, doi: 10.1007/
Alain de Cheveigné. Sparse time artifact removal. Journal of Neuroscience Methods, 262:14-20 2016. doi: jneumehh.n..
Alain De Cheveigné. ZapLine: A simple and effective method to remove power line artifacts. Neurolmage, 207:1 16356, 2020. doi: ner.000ee. 16356.
Alain de Cheveigné and Lucas C. Parra. Joint decorrelation, versatile tool for multichannel data analysis. Neurolmage, 98:487-505, 2014. doi: 10.1016/j neroommg...o00.0...0..5.
Stefan Debener, Markus Ullsperger, Markus Siegel, and Andreas K Engel. Single-trial eeg-fmri reveals the dynamics of cognitive function. Trends in cognitive sciences, 10(12):558-56) 2006.
Fani Deligianni, Maria Centeno, David Carmichael, and Jonathan Clayden. Relating resting- state fmri and eeg whole-brain connectomes across frequency bands. Frontiers in neuroscience, 8:258, 2014.
Fani Deligianni, David W Carmichael, Gary H Zhang, Chris Clark, and Jonathan D Clayden. Noddi and bassor-based microstructural indices as predictors of functional connectivity. Plos one, 11(4):e115340, 2016.
Sven Dähne, Vadim V. Nikulin, David Ramírez, Peter J. Schreier, Klaus-Robert Müller, and Stefan Haufe. Finding brain oscillations with power dependencies in neuroimaging data. Neurolmage, 96:334-348, 2014. doi: 10. 1100ff..eerrimgee.004.0000000000
Sven Dähne, Felix Bießman, Wojciech Samek, Stefan Haufe, Dominique Goltz, Christopher Gund- lach, Arno Villringer, Siamac Fazli, and Klaus-Robert Müller. Multivariate Machine Learning Methods for Fusing Functional Multimodal Neuroimaging Data. Proceedings of the IEEE, 103 (9):1-22, 2015. doi: 10.1 10//PP00...005..2225077

---

Published as a conference paper at ICLR 2024
Alexandre Défossez, Charlotte Caucheteux, Jérémy Rapin, Ori Kabeli, and Jean-Rémi King. De- coding speech perception from non -on-iivvsivv brain recordings. Nature Machine Intellignnce, 2023. doi: 10.1038//42256-022-007744,
Alcksandr Ermolov, Aliaksandr Siarohin, Enver Sangineto, and Nicu Sebe. Whitening for self- supervised representation learning. In International Conference on Machine Learning, pp. 3015- 3024. PMLR, 2021.
Oscar Esteban, Christopher Markiewicz, Ross W. Blair, Craig A. Moodie, A. Ilkay Isik, Asier Erramuzpe, James D. Kent, Mathias Goncalves, Elizabeth DuPre, Madeleine Snyder, Hiroyuki Oya, Satrajit S. Ghosh, Jessey Wright, Joke Dumez, Russell A. Poldrack, and Krzysztof J. Gor- golewski. IMRIPrep: a robust preprocessing pipeline for functional MRI. Nature Methods, 16 (1):111-116, 2019. doi: 10.1038644152--08--027
Stephen H. Fairclough and Fabien Lotte. Grand Challenges in Neurotechnology and System Neu- roergonomics. Frontiers in Neuroergonomics, 1:602504, 2020.
Faczeh Fallah and Bin Yang. Supervised canonical correlation analysis of data on symmetric pos- itive definite manifolds by riemannian dimensionality reduction. In ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 0369-8373. IEEE, 2020.
Pierre Fillard, Xavier Pennec, Vincent Arsigny, and Nicholas Ayache. Clinical dt-mri estimation, smoothing, and fiber tracking with log-euclidean metries. IEÉE Transactions on Medical Imag- ing. 26(1 1):1422-1482, 2007.
Thomas Fletcher and Sarang Joshi. Riemannian geometry for the statistical analysis of diffusion tensor data. Signal Processing, 87(2):250-262, 2007.
Thomas Fletcher, Conglin Lu, Stephen M Pizer, and Sarang Joshi. Principal geodesic analysis for the study of nonlinear statistics of shape. IEEE transactions medical imaging, 2358) 2004.
Thomas Fletcher. Geodesic regression riemannian manifolds. In Proceedings of the Third Inter- national Workshop on Mathematical Foundations of Computational Anatomy-Geometee and Statistical Methods for Modelling Biological Shape Variability, pp. 75-86, 2011.
Maurice Fréchet. Les éléments aléatoires de nature quelconque dans un espace distancié. In Annales de I l'institut Henri Poincaré, volume 10, pp. 215-310, 1948.
Lukas AW Gemein, Robin T Schirrmeister, Patryk Chrabkaszez, Daniel Wilson, Joschka Boedecker, Andreas Schulze-Bonbage, Frank Hutter, and Tonio Ball. Machine-leaming--aate based diagnostics of eeg pathology. Neurolmage, 220:1 17021, 2020.
Karsten Grove and Hermann Karcher. How to conjugate 1-close group actions. Mathematische Zeitschrift, 132(1):11-22) 1973.
Nicolas Honnorat and Mohamad Habes. Covariance shrinkage can assess and improve functional connectomes. Neurolmage, 256: 255:119228 19229, 2022.
Harold Hotelling. Relations between two sets of variates. In Breakthroughs in statistics: methodol- ogy and distribution, pp. 162-190. Springer, 1992.
Zhiwu Huang and Lue Van Gool. A Riemannian Network for SPD Matrix Learning. In Proceedings of the Thirty- Marrt-FFirttaa AAAI Conference on Artificial Intelligence, AAAI'17, pp. 2036-2042. AAAI Press, 2017. event-place: San Francisco, California, USA.
Zhiwu Huang and Lue Van Gool. riemannian network for spd matrix learning. In Thirty-Fiie Firstt AAAI Conference Artificial Intelligence, 2017.
René Huster, Stefan Debener, Tom Eichele, and Christoph Hermmann. Methods for simultancous eeg-fmri: an introductory review. Journal of Neuroscience, 32(1 (2(18)600536060, 2012.
12

---

Published as a conference paper at ICLR 2024
Catalin lonescu, Orestis Vantzos, and Cristian Sminchisescu. Matrix Backpropagation for Deep Networks with Structured Layers. In 2015 IEEE International Conference on Computer Vision (ICCV), pp.
IEEE, 2015. doi: 10.11099/CCV.2015.3399
Mainak Jas, Denis A. Engemann, Yousra Bekhti, Federico Raimondo, and Alexandre Gramfort. Autoreject: Automated artifact rejection for MEG and EEG data. Neurolmage, 159:417-429, 2017. doi: 10. new.01i..eeroimmgeee...... 201..0...00..0000000
Vinay Jayaram, Morteza Alamgir, Yasemin Altun, Bernhard Scholkopf, and Moritz Grosse- Wentrup. Transfer Learning in Brain- Computer Interfaces. IEEE Computational Intelligence Magazine, 11(1):20-31, 2016. doi: 1090//..001..2500000
Li Jing, Pascal Vincent, Yann LeCun, and Yuandong Tian. Understanding dimensional collapse in contrastive self-supervised learning. In International Conference Learning Representations, net/ rew..//orrrnn.eeeeessee....
2022. URL
Ce Ju and Cuntai Guan. Tensor-cspnet: A novel geometric deep learning framework for motor imagery classification. IEEE Transactions on Neural Networks and Learning Systems, 2022.
Ce Ju and Cuntai Guan. Graph neural networks spd manifolds for motor imagery classification:
A perspective from the fen-frreyenccy analysis. IEEE Transactions on Neural Networks and Learning Systems, 2023.
Jackob N. Keynan, Avihay Cohen, Gilan Jackont, Nili Green, Noam Goldway, Alexander Davi- dov, Ychudit Meir-Hasson, Gal Raz, Nathan Intrator, Eyal Fruchter, Keren Ginat, Eugene Laska, Mare Cavazza, and Talma Hendler. Electrical fingerprint of the amygdala guides neu- rofeedback training for stress resilience. Nature Human Behaviour, 3(1):63-73, 2019. doi: 10. 00.10088441562-018-444833
Hyunwoo J. Kim, Nagesh Adluru, Barbara B. Bendlin, Sterling C. Johnson, Baba C. Vemuri, and Vikas Singh. Canonical correlation analysis on riemannian manifolds and its applications. Com- pul Vis ECCV, (8690):251-667, 2014.
Diederik P. Kingma and Jimmy Ba. Adam: A Method for Stochastic Optimization, 2017.
Reinmar Kobler, Jun-ichiro Hirayama, Qibin Zhao, and Motoaki Kawanabe. SPD domain batch normalization to crack interpretable unsupervised domain adaptation EEG. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (eds.), Advances in Neural Information Processing Systems, volume 35, pp. 6219-6235. Curran Associates, Inc., 2022a.
Reinmar J. Kobler, Jun-ichiro Hirayama, and Motoaki Kawanabe. Controlling The Fréchet Variance Improves Batch Normalization on the Symmetric Positive Definite Manifold. ICASSP 2022 2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp- 3863-3867, Singapore, Singapore, 2022b. IEEE, doi: 10. 101111CASSP449922222227777
Max Kochurov, Rasul Karimov, and Serge Kozlukov. Geoopt: Riemannian Optimization in PyTorch, 2020.
Serge Lang. Fundamentals of differential geometry, volume 191. Springer Science & Business Media, 2012.
Xueqing Liu, Tao Tu, and Paul Sajda. Inferring latent neural sources via deep transcoding of simul- ott..r/gaabbb......
taneously acquired EEG and fMRL. arXiv, 2022. URL
02226
Maher Moakher. A differential geometric approach to the geometric mean of symmetric positive- definite matrices. SIAM journal on matrix analysis and applications, 26(3):735-74) -4777 2005.
Moeller, L. Tyvaert, D. K. Nguyen, LeVan, A. Bouthillier, E. Kobayashi, D. Tampieri, F. Dubeau, and J. Gotman. EEG-IMRI: Adding to standard evaluations of patients with non- lesional frontal lobe epilepsy. Neurology, 73(23):2023-25 2030, 2009. doi: 10.1212/WNL. 0b013e31 1155d117.
13

---

Published as a conference paper at ICLR 2024
R.K. Niazy, C.F. Beckmann, G.D. lannetti, J.M. Brady, and S.M. Smith. Removal of FMRI environ- ment artifacts from EEG data using optimal basis sets. Neurolmage. 28(3):720-737, 2005. doi: 10.1016/j 10.1016/j.meurommaag wtwp://aa......0..00000
Takeshi Ogawa, Hiroki Moriya, Nobuo Hiroe, Motoaki Kawanabe, and Jun ichiro Hirayama. Eeg- based neurofeedback with network components extraction: data-driven approach by multilayer ica extension and simultaneous eeg-fmri measurements. bioRxiv, 2021. doi: 10.1101/2021.0..2 449196.
Wanmei Ou, Aapo Nummenmaa, Jyrki Ahveninen, John W Belliveau, Matti S Hämäläinen, and Polina Golland. Multimodal functional imaging using fmri-informed regional eeg/meg source estimation. Neuroimage, 52(1)997-108, 2010.
Hae-Jeong Park and Karl Friston. Structural and functional brain networks: from connections to 2013.
cognition. Science, 342(6158):
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming I Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang. Zachary Veto,o Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. PyTorch: An Imperative Style, Hirf--Peforraace Deep Learning Library. In H. Wallach, H. Larochelle, A. Beygelzimer, d' Alché-Bus E. Fox, and R. Garnett (eds.), Advances in Neural Information Processing Systems 32, pp. 8024 8035. Curran Associates, Inc., 2019.
Pedregosa, G. Varoquaux, Gramfort, V. Michel, Thirion, O. Grisel, M. Blondel, Pretten- hofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine Learning in Python. Journal of Machine Learning Research, 12:2825-2830, 2011.
Xavier Pennec. Manifold-vad -aauue image processing with spd matrices. In Riemannian geometric statistics in medical image analysis, pp. 75-134. Elsevier, 2020.
Xavier Pennec, Pierre Fillard, and Nicholas Ayache. A riemannian framework for tensor computing. International Journal of computer vision, 66:41-66, 2006.
Gert Pfurtscheller and FH Lopes Da Silva. Event-related eeg/meg synchronization and desynchro- nization: basic principles. Clinical newrophysiology, 110(11)1142--857) 1999.
Marios G. Philiastides, Tao Tu, and Paul Sajda. Inferring Macroscale Brain Dynamics via Fusion of Simultaneous EEG-IMRI. Annual Review of Neuroscience, (4(1):3153334, 2021.
Luca Pion-Tonachins Ken Kreutz-Delgado, and Scott Makeig. ICLabel: An automated electroen- cephalographic independent component classifier, dataset, and website. Neurolmage, 198:181- 197. 2019. doi: www0l//j.nuroimmae..oo 2019.05.026.
Anqi Qiu, Annie Lee, Mingzhen Tan, and Moo K Chung. Manifold learning on brain functional networks aging. Medical image analysis, 20(1):52-60, 2015.
Pedro Luiz Coelho Rodrigues, Christian Jutten, and Marco Congedo. Riemannian Procrustes Anal- ysis: Transfer Learning for Cwain-Computee Interfaces. IEEE Transactions Biomedical Engi- neering, 6638):2390- 2019. doi: 28.1009//BEE..11.28888888
Raphaälle Roy, Marcel Hinss, Ludovic Darmet, Simon Ladouce, Emilie S Jahanpour, Bertille Somon, Xiaoqi Xu, Nicolas Drougard, Frédéric Dehais, and Fabien Lotte. Retrospective the first passive brain-computer interface competition on cross-session workload estimation. Fron- fiers in Neuroergonomiis, 3:838342, 2022.
David Sabbagh, Pierre Ablin, Gaël Varoquaux, Alexandre Gramfort, and Denis A Engemann. Pre- dictive regression modeling with meg/cog: from source power to signals and cognitive states. Neurolmage, 222:1 16893, 2020.
Sepideh Sadaghiani and Jonathan Wirsich. Intrinsic connectome organization across temporal scales: New insights from cross-modal approaches. Network Neuroscience, 4(1):1-29, 2020. doi: 10.1162/netn 001 14.
14

---

Published as a conference paper at ICLR 2024
Armin Thomas, Christopher Ré, and Russell Poldrack. Self- Self-Super Learning of Brain Dynamics from Broad Neuroimaging Data. In S. Koyejo. S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (eds.), Advances in Neural Information Processing Systems, volume 35, pp. 21255- 21 21269. Curran Associates, Inc., 2022.
Thomas Fletcher. Geodesic regression and the theory of least squares on riemannian manifolds. International journal of computer vision, 105:171-1811 2013.
Viivi Uurtio, João M Monteiro, Jaz Kandola, John Shawe-Taylor Delmiro Fernandez-Reyee and Juho Rousu. A tutorial on canonical correlation methods. ACM Computing Surveys (CSUR), 50 (6):1-33, 2017.
Johan van der Meer, André Pampel, Eus van Someren, Jennifer Ramautar, Ysbrand van der Werf, German Gomez-Herrero, Jöran Lepsien, Lydia Hellrung. Hermann Hinrichs, Harald Möller, and Martin Walter. "Eyes Open Eyes Closed" EEG/MMI data set including dedicated "Carbon Wire Loop" motion detection channels. Data in Brief. 7:990-999 994, 2016a. doi: 11.0016...ibb2016.000 001.
Johan N. van der Meer, André Pampel, Eus J.W. Van Someren, Jennifer R. Ramautar, Ysbrand D. van der Werf, German Gomez-Herrero, Jöran Lepsien, Lydia Hellrung, Hermann Hinrichs, Har- ald E. Möller, and Martin Walter. Carbon-wint loop based artifact correction outperforms post- processing EEG/MMRI corrections AA validation of a real-time simultaneous EEG/MMI correc- tion method. Neurolmage. 125:880-894 2016b. doi: 10.10166/. www.:/maae...21.........
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, łukasz. Kaiser, and Illia Polosukhin. Attention U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (eds.), Advances in Neural Information Processing Systems. volume 30. Curran Asso- ciates, Inc., 2017. URL https::/pproceeddiiiii heweeeiiiisnnnnninneeeeeeee www.inuu.i.... wttpuuupps..cpppppreeeeepepe cttp://ppee//217//1eeeee.
you Need. In Guyon,
Weiran Wang, Raman Arora, Karen Livescu, and Jeff Bilmes. On deep multi-view representation learning. In International conference machine learning, pp. 1083-10922 PMLR, 2015.
Tracy Warbrick. Simultaneous EEG-IMRI: What Have We Learned and What Does the Future Hold? Sensors, 22(6):2262, 2022. doi: 10.3390/s2206226..
Daniel Wilson, Robin Tibor Schirmeister, Lukas Alexander Wilhelm Gemein, and Tonio Ball. Deep Riemannian Networks for EEG Decoding, 2022. URL http:://arx.v.org/aa/////2
10426
Jonathan Wirsich, Anne-Lise Giraud, and Sepideh Sadaghiani. Concurrent EEG- and IMRI-derived functional connectomes exhibit linked dynamics. Neurolmage, 219:116998. 2020. doi: 10.1016/ j.neuroimage.2000 2020.. 116998.
Dongrui Wu, Yifan Xu, and Bao-Liang Lu. Transfer Learning for EEG-Based Brain-Computer Interfaces: A Review of Progress Made Since 2016. IEEE Transactions on Cognitive and Devel- opmental Systems, pp. 1-1, 2020. doi: 2.019//TTCDDS2220000744..3.
Or Yair, Mirela Ben-Chen, and Ronen Talmon. Parallel Transport the Cone Manifold of SPD Matrices for Domain Adaptation. IEEE Transactions on Signal Processing, 6777)1197-1811, 2019. doi: 10.1109/SSP20192299499
Ayumu Yamashita, Yuki Sakai, Takashi Yamada, Noriaki Yahata, Akira Kunimatsu, Naohiro Okada, Takashi Itahashi, Ryuichiro Hashimoto, Hiroto Mizuta, Naho Ichikawa, et al. Generalizable brain network markers of major depressive disorder across multiple imaging sites. PLoS biology, 18 (12):e3000966, 2020.
Kisung You and Hae-Jeong Park. Re-visiting riemannian geometry of symmetric positive definite matrices for the analysis of functional connectivity. Neurolmage, 225:117464 2021.
15

---

Published as a conference paper at ICLR 2024
Paolo Zanini, Marco Congedo, Christian Jutten, Salem Said, and Yannick Berthoumieu. Transfer Learning: A Riemannian Geometry Framework With Applications to Brain Computer Interfaces. IEEE Transactions on Biomedical Engineering, 65(5):100711116, 2018. doi: 10.1109/TBME.
2017.2742541.
Jure Zbontar, Li Jing, Ishan Misra, Yann LeCun, and Stéphane Deny. Barlow twins: Self-supervised learning via redundancy reduction, In International Conference on Machine Learning, pp. 12310- 12320. PMLR, 2021.
16

---

Published as a conference paper at ICLR 2024
A SPD MANIFOLDS WITH AFFINE INVARIANT RIEMANNIAN METRIC A
In the section, we will elucidate covariance-based neuroimaging data on the SPD manifold. First, we present its definition follows,
Definition (Covariance-Basee Neuroimaging Data). The covariance matrix of neuroimaging data segment, represented by X I xee, is denoted as S = XX /nr S44,, where ne is the number of spatial dimensions and wer the number of timesteps.
Covariance-basse data is inherently symmetric and positive definite. Even for degenerate cases (e.g., co-linear spatial dimension), the positive definiteness can be ensured via shrinkage regulariza- tion Honnorat Habes (2022).
Note that this definition holds for any multi-channel neuroimaging time series. However, in this study, we focus on two neural signal modalities, namely EEG spatial covariance matrices and fMRI functional connectivities. EEG spatial covariance matrices capture relevant aspects of neural ac- tivity within short to long-time scales. Riemannian geometry-aware methods operating with these covariance-bases data are considered state-of-the-an in various EEG application domains, including clinical applications (Sabbagh et al. 2020) Gemein et al.) 2020) and recent brain-computer interface competitions (Roy et al., [2022). Concerning fMRI, covariance (or correlation) matrices between regions of interest reflect static and dynamic communication between brain networks and their al- terations in tasks (Park & Friston 2013) and pathological conditions (Yamashita et al. 2020).
Equipped with the affine invariant Riemannian metric g ATRM, as follows,
gp ATRM (v, w) = (P- PP-pPPP ) TT(P-'uP--u0
for any langent vectors v and in the tangent space TpS++ + at the base point P S++, the space of covariance-bast neuroimaging data becomes an SPD manifold, denoted as (S++. X AiRMM, ), which Hadamadd manifold. A Bartan-Hadamart manifold is a Riemannian manifold that is complete, connected, and possesses everywhere non-positive sectional curvature. In particular, Riemannian distance between two matrices S₁ and S2 defined as follows,
dyarnsr (S1, S2) = ( (ST -S2)
B FRÉCHET MEAN AND VARIANCE
In this section, we provide a brief introduction to the Fréchet mean and variance (Frechet 1948), also known as Karcher mean (Grove & Karcher 1973). The two names are often used interchangeably in literature.
Given set of points (PI....,P PN E (M,g) g). the Fréchet mean μ is defined as the minimizer of the sum-of-squared Riemannian distances as follows.
arg min PEM
(5)
And, the corresponding Fréchet variance is defined as follows,
N 0² = E = N Σ (..
(6)
Remark. We consistently employ parallel transport to relocate the Riemannian means of each data cluster to the identity matrix e. Consequently, the Fréchet variance can be simplified as follows:
N
1
2
llog, (5:)
N
which is given by dystrer S) = |loog. (S) where Riemannian norm ||p||__ := g.. AARR (v for any E T.S++
17

---

Published as a conference paper at ICLR 2024
C ORTHOGONAL PROJECTION AND GEODESIC CONSTRAINTS
This section first provides an clementary introduction to the geodesic on Riemannian manifolds.
Let (M" .g) be Riemannian manifold with dimension equipped with the Levi-Civita connection
V. and y + .M be a smooth curve on M. We say a vector field X is parallel along y V-rt = for all & R. A map Pla.. ( + T M is parallel transport of Xo M along y if Pi (X((()) is the parallel vector field along (t) such that ((>(0)) Xo-
We say (t) is geodesic if (t) is parallel along 7(),, ie.,
7 0.
Note that ( (t) = 0 enforces ((t) to be curved with zero acceleration. The concept geodesic extends the idea straight lines from Buclidean spaces. is worth mentioning that every geodesic is parameterized in proportion to its are length.
In local chart ()) the geodesic curve 7(1) satisfies the following second-order ordinary differ- ential equations:
(t) + Σ 1 (t))(()) 0, for k 1,,.7 R, 19
where each (1) = & o (t). For local frames {0: 1 denotes the Christoffel symbols of V given by Va, (a,) Or. Vi,j, k- = 1,..., n. By the Picard-Lindelöf theorem, it yiclds the local existence and uniqueness of the solution for geodesic: For any given p E M and E Tp.M, there exists an interval (-c,c) c) and unique geodesic (-c,c) b M satisfying the initial conditions y(0) = and 7((0) := V.
Given that the SPD manifold (S++ g A/RM) is geodesically complete, meaning that all geodesics are defined for all time / E R. any geodesic exp, (lu) on S++ is solely determined by u E T and varies smoothly with respect u.
To clarify the geodesic constraints, we will begin by introducing the idea orthogonal projection: given unit speed E S" the orthogonal projection of E S on arc-length geodesic exp, (fu) 1. ALRM) is St exp. (t' with the following optimality condition:
arg min d ArRR (S, (t))
BER
Plugging into Equation
4 Equation becomes into the following expression:
Tr 1r(1og' loge ( exp(
arg min (E(-6,))
(8)
Take the derivative of Equation with respect to time I, set to zero, and we obtain:
Tr loge (S-1. exp)(
Ht
2 Ht exp, (tu)
Tr log. (S-1 (tp. exp.(-tu).
which yields an equality as follows,
Tr TT(1og. logg (
cxp. (tu)) expe(-tu) expe 0.
(9)
Equality (*) is valid by virtue of the following proposition:
Proposition 1 (Moaheer (2005)). Consider real matrix-valued function S(t) defined on time t IR. Assume that S(t) is an invertible matrix and does not have eigenvalues on the closed negative real line. Then, have the following equality:
at T1
i Tr log2 S(t) 2- Tr log, S(t) S(t)-
s(t))
at
3 A covariance matrix S marked with t symbol in the upper right corner denoted as St, signifies that it represents the orthogonal projection of S on the geodesic.
18

---

Published as a conference paper at ICLR 2024
Suppose := tf is optimal and thus St exp. (tpx(tlu). Hence, Equality 9 Tr (log, (S-1. s!) St u St = 0.
becomes to be as follows,
(10)
According to the problem setting given in Definition Equation 10 yields the following two geodesic constraints of paired covariance based neuroimaging data for all E NI,
Tr Tr(loge logg (5 .s) S 5 0;
Tr logg (5.1 5!) S w 5 0.
Remark. In this remark, we briefly discuss the reason for introducing geodesic and the orthogonal projection on the geodesic. The concept of geodesic constraint originates from the field of Geomet- ric Statistics, referred to as the geodesic submanifolds, initially introduced in principal geodesic
analysis (PGA) and related works Fletccher et al. 2004 Fletcher & Joshi 2007 Thomas Fletcher 2013).
In the context of Riemannian geometry, geodesic is a curve that locally represents the shortest path between points, acting as a generalization of a straight line. Moreover, the geodesic subspace can be viewed as an expansion of linear spaces from Euclidean to Riemannian geometry. The concept of the geodesic subspace is pivotal in PGA as the projection on this subspace preserves the Riemannian distance, representing the munifold-valued data's variance.
In a formal sense, a submanifold H of manifold M is considered geodesic at the point p E H if all geodesics of N that pass through p remain geodesic of M. A geodesic submanifold H at the Fréchet mean J. is represented as the exponential mapping of the linear span of tangent vectors {w; TTuM as N e (ppan(ww:)==
In PGA, tangent vectors (w=)._1 are computed by sequentially maximizing the Riemannian distance between each data point and the Fréchet mean. This process captures the variance of the data at each step while progressively removing the inftuence of previously derived directional projections. This completes the manifold valuee data's 's dimensionality reduction. In this context, w1 can be analo- gous to the first principal component in principal component analysis (PCA). This work goes beyond finding the geodesic submanifold's 's tangent vectors { K. obtaining novel low oow-dimesiionl data representation is crucial.
CCA aims to identify the maximum correlation between two variables. Specifically, as outlined in Section 2.2, CCA aims to find linear transformations for paired matrices E NRN and X E RNN denoted as w E RR" and w E R9, respectively, to maximize the correlation between these two linear combinations Xww and X w. CCA is typically transformed into an eigenvalue de- composition problem using the Lagrange multipliers method. Consequently, the objective often involves identifying a set of r pairs of canonical variables along with corresponding linear trans- formations W R.x and W Raxr, where each column in these matrices signifies mapping to the respective pair of the canonical covariate. A constraint is imposed to render them uncorrelated one another ensure that each pair of canonical variables captures distinct phenomena.
When extending the CCA method to Riemannian manifolds, Kim et al. (2014) was the first to use geodesic submanifolds as a replacement for linear subspaces. They did not underscore the or- thogonal constraints for each pair of canonical covariates and obtained them through a nonlinear optimization approach. Each resulting pair of canonical covariates corresponds to tangent vec- tors [wi)-= and [w; Mm for two geodesic submanifolds. Therefore, orthogonal projection on the geodesic can be seen as projecting the data onto a lower- dowver-dmensional space by analogy with a linear transformation in classical CCA. Each langent vector in the langent space corresponds to column vector in the linear transformation (canonical covariates).
The initial target of this study was to enhance the temporal dynamic consistency of simultaneously recorded EEG-fMRI signals. We opt for only one geodesic path for each modality to streamline the method. This design simplification indeed leads to well-defined convex problem and computational advantages in our deep learning architecture. In particular, opting for only one geodesic path for each modality is equivalent utilizing the first pair of canonical covariates in CCA. In other words, the tangent vector Wn a symmetric matrix (tangent vector for SPD manifolds), precisely serves as the first canonical covariate. Thus, while the mechanism for obtaining wi differs from that obtained
19

---

Published as a conference paper at ICLR 2024
through eigenvalue decomposition in CCA DeepCCA, our proposed component analysis method inherently possesses a directional aspect.
D COMPARISONS BETWEEN DIFFERENT CCAS
In Table 3
we analyze the similarities and differences between several related variants of the CCA method and our proposed method from three perspectives, including:
Nonlinear Projection: Projections on the nonlinear subspace not.
Eigendecomposition: Utilizing the eigendecomposition method not, i.e. eigendecompo- sition of
Σ>. 1/2
>
1/2
Σ Exp
Eppe1 in Subsection 2.2 and
1/2
2.3.
Neural Networks-Bas -Bsse Approach: the neural networks-based approach or not.
Geometry: SPD matrix-valued outputs or not.
Table 3: Comparisons between different CCAs.
Methodology
Nonlinear Projection
Eigendecomposition
Neural Networks
Geometry
CCA
DeepCCA
V
RieCCA
X
DeepGeoCCA
V
X
Please note that RicCCA is solved analytically through approximation approach within the log-
Euclidean framework using CCA (Kim et al. 2014). Therefore, we say that RicCCA involves
eigendecomposition.
CORRELATION, RIEMANNIAN CORRELATION, AND GEODESIC CORRELATION
In this chapter, we will provide detailed explanation of the distinctions between correlation, Rie- mannian correlation, and geodesic correlation. Jillutrtelss the distinctions between the three terms. All of them are utilized to elucidate statistical relationships between two random variables, whether causal or not. However, their definitions encompass different scopes, as outlined below:
Correlation: it is a general statistical term for any random variable (Croxton 1967).
Riemannian Correlation: it is a generalized term of correlation for general Riemannian manifolds. quantified by computing the correlation between tangent vectors, which represent manifold-valus v-valeed data projected onto the respective geodesic submanifold of the tangent spaces the Fréchet mean of each modality (Kim et al. 2014).
Geodesic Correlation: is a generalized term of correlation for covariance-bari neu- roimaging data. It is quantified by calculating the correlation between tangent vectors, which represent manifold-valued data projected onto the respective geodesics of the tan- gent spaces of each modality at the identity matrix
The main motivation for introducing geodesic correlation is due to the fact that (S++, AIRM) is Cartan-Hadamard manifold that is complete, connected, and has an everywhere non-positive sectional curvature (Pennec, 2020). Consequently, (S++. A.Rt. ALRAS) geodesically complete, and the geodesics can be extended to the whole IR. These excellent geometric properties enable the simplification of Riemannian correlation on the SPD cone. We utilize this fact as follows:
Geodesics: We simplify the constraint geodesic submanifold Riemannian correlation to geodesic.
20

---

Published as a conference paper at ICLR 2024
expiir)) esptred (aptrng) casuiey)
esp roperst
exp,(r))
exp(rw)
(eg))
M3
(xp
Mx
Figure 5: Illustrations of Correlation, Riemannian Correlation, and Geodesic Correlation: From left to right, have Correlation, Riemannian Correlation, and Geodesic Correlation with the spe- cific descriptions. Correlation: Correlation is only required to be demonstrated between paired data; Between each pair of paired data, there is no need for proximity; Riemannian Correlation: It computes the correlation between paired projection data on the tangent spaces of each modal- ity, considering their respective Fréchet means. The projection data is obtained by projecting the raw data onto geodesic submanifolds; Geodesic Correlation: It calculates the correlation between paired projection data in the langent space the identity point The projection data obtained by projecting the raw data onto or near geodesics.
Base Point of Tangent Space: We exclusively compute the correlation the tangent space at the identity point e, rather than the Fréchet mean for each modality. In particular, defining geodesic correlation at equivalent defining it at any point s E S++. This equivalence is established through parallel transportation (see Appendix G.2).
Correlation Measurement: We use Riemannian metric to measure the correlation between two tangent vectors on the tangent space the identity point In contrast, Riemannian Correlation relies on the inner product of scalar values of tangent vectors based on different Fréchet means.
Given these simplifications, geodesic correlation fulfills the requirements of neural network- based solution. We summarize the primary distinctions from the four aspects discussed earlier Table as follows:
Table 4: Difference between Correlation, Riemannian Correlation and Geodesic Correlation
Correlation
Region
Base Point
Correlation Measurement
Solver
Correlation
Lines
Dot product
CCA
Riemannian Correlation
Geodesic submanifolds Geodesics
Two Fréchet means
Scalar dot product
RieCCA
Geodesic Correlation
Riemannian metric
DeepGeoCCA
Geodesic correlation has relationship with correlation according to the following proposition:
Proposition 2. Suppose and w are tangent vectors on T.S of ( then we have
AIRM (v, w) = (vec(u), vec(w))) ge
Proof. This proposition arises from the following derivation,
AIRM (v, w) == Tr(v- - w) = Σ oii, w[: o = (vec(v)) )) vec(w) e (+) ge
(vec(v), vec(w)) )
AIRM) is a symmetric where (*) is derived from the property that the tangent space of ASR.. space.
This proposition demonstrates that maximizing the geodesic correlation on g AIRM) at the base point is equivalent maximizing the correlation if view the covariance-basse data as
21

---

Published as a conference paper at ICLR 2024
flattened tensor. Furthermore, this proposition also implies that the coefficient of determination
R in the Riemannian geometry- based sense, as discussed in Fletcher 2011). 1) is equivalent to the classical R².
F CONSTRAINED NONLINEAR OPTIMIZATION PROBLEM
This section will show how to convert geodesic correlation into a constrained nonlinear optimization problem in matrix and Euclidean space.
In matrix space, the constrained nonlinear optimization problem is depicted as follows,
Σn ge
AIRM
(loge (S.), log. (S!))
minimize U, w. S..
E logo (S
|1: 2
2
(Ei=
|ogg (S (5))
(11)
subject to
Tr log, (5 .S! S vv. -s
=0
Tr logg (5,1 5/) 5t w 5t =0
with variables v, E S" and st. E S++ for i I,,...N N. The consistency targeted by this constrained nonlinear optimization problem using correlation on the tangent space is referred to as geodesic correlation.
The constrained nonlinear optimization problem can also be written in Euclidean space. However, when the results obtained through constrained nonlinear optimization algorithm in Euclidean space are converted back into matrices, they cannot he guaranteed to be positive definite and symmetric.
Next, we will derive the constrained nonlinear optimization problem in Euclidean space, equivalent to Problem II in the matrix space.
We substitute the correlation (84 AIRM) with conventional correlation denoted follows, vee( log, (S.)) and tu vee( log, (S.)). for i I...,N.. N. Then, we obtain the constrained nonlinear optimization problem in Euclidean space with total of 2N constraints, depicted as follows,
(Po, ( EM 11al) (E 11a)
minimize U, w. by tu subject
Tr (log, (ST (t. expp(-t))) exp. (t.u))
(12)
Tr log. (5. cxp, (f w) exp, (-th - w) w (tp
with variables v, w E S and H, t E RRn for = 1, N.
By solving Constrained Nonlinear Optimization Problem (2) through constrained nonlinear opti- mization algorithms, even when the resulting t, and tw are converted into a matrix, it no longer maintains the positive definite symmetric properties.
In the following, we will derive the Lagrangian for Constrained Nonlinear Optimization Prob- lem (11) in the matrix space. We arrive the expression for Lagrangian L(ST, SS ,, u,w, A,7, : ST ST x S" x St x EN x R + R, as follow,
En AIRM (log, (S!), log, (S!))
L(St, St, U. U.W.W A. X, 7)
N D
log. (S)) [3) -EN |ogg (p) (5!)
N A Tr (oog (5 ST) st - ss)
N
7s Tr loge (5.1 -T) 5. w St).
...s))
i=1
22

---

Published as a conference paper at ICLR 2024
where {Ax R) and {% R) Lagrange multipliers.
In the context that follows, we refer to the orthogonal projection St or St as Z or Z, respectively. These orthogonal projections are obtained by neural network-basee training using gradient descent. In practical implementation, these values precede the input to the LOG layer. The results produced by the LOG layer are denoted as log, (Z:) and log, (Zi).
Hence, the Lagrangian duality is naturally a loss function for DeepGeoCCA, depicted as follows,
Σ Tr( (looge (z:) log, (2.))
C(Z,2):
V Zi
N
(oo0 V V Tr( log (2i))
N
+
Σ
A;- Tr TD((10gg) (log. (5.² Zi) Z v Zi)
2=1
N
+
Σ
Ni Tr? log, (5,1 Zi) Z w 2.).
8=1
where ds, E R are preset coefficients, for i = I...., N, and v, E S" are preset tangent vectors.
G NEURAL NETWORKS ON SPD MANIFOLDS
The SPD matrix-valued network layers within the self-supervised learning framework in the pro- posed architecture are outlined in the following three studies:
G.1 SPDNET
Huang & Van Gool (2017) introduced the SPDNet network architectures which specifically de- signed to operate with SPD matrices to ensure that the essential properties of symmetry and positive definiteness are maintained throughout the learning process. It incorporates the following layers:
BiMap: This layer performs the bi-map transformation W on S. The transformation matrix W is typically required to have a full-row rank.
ReEig: This layer performs max(c e, Σ) UT.
LOG: This layer maps to its tangent space at the identity matrix using U Ulogg (2) U.
where signal covariance matrix = -..
G.2 BATCH NORMALIZATION
Brooks et al. 2019 conceived a layer, termed Riemannian BN, form of SPD matrix-valueds network architecture to center and re-bias batches of latent SPD data B. Formally, the computation of the batch Fréchet mean, denoted B, involves solving (5) and utilizing parallel transport (Yair let al. (019)) to establish connections between each element E B and the identity matrix e. This is achieved through the following operations:
T BB [B = B Z B
T To-o(Z) G. 7 G!
where G ST is learnable biasing parameter that is optimized through matrix backpropagation during training. To manage computational overhead per minibatch, Brooks et al. (2019) performed only step of the iterative update algorithm to solve (5). During testing, the layer uses the dataset's estimated Fréchet mean, computed via exponential smoothing of the batch means.
To mitigate potential limitations due to erude Fréchet mean estimates Kobler et al. (202aa) proposed an iterative update scheme, denoted SPDMBN, that can control the error to the latent data's 8 changing
Letter Z commonly recognized as symbols representing latent representations in deep learning.
23

---

Published as a conference paper at ICLR 2024
Fréchet mean during matrix backpropagation training and, upon convergence of preceding feature extraction layers. converge to the latent data's true Fréchet mean. They also showed that if the biasing parameter is kept fixed to the identify matrix (G e), the combination of SPDMBN and LOG layers can learn to transport the SPD data to vary around e and project the data to the tangent space T
Different from Riemannian BN, SPDMBN can also learn to re-scale the SPD data via keeping track of the dataset's Fréchet variance, defined in (6, which aids the matrix backpropagation learning process (Kobler et al. 2022b).
H RELAXATION DEVIATION
Proof. The 5-width tubular neighborhood of a geodesic (lu) on g AARM always exists
inside its normal bundle according to general tubular neighborhood theorem Lang 2012).
For any predefined unit speed TeST+. let the subspace of speed vectors U {tul IR} C
T +- We define the orthogonal complement of li as u- = (u-|gAIRM (u u) 0) C
TaST+ Note that the orthogonal complement 4 not an empty set, as the space of skew-
symmetric matrices is the orthogonal complement of the space of symmetric matrices.
For each i [I. suppose log, are orthogonal components of log, St with respect as follows,
(S.)) = log, (S:) ge AIRM (l.g.(S.): u) E
Assuming E (0, represents the deviation magnitude, is defined as magnitude with the direction that perpendicular to the vector that passes through (S)) in the tangent space. The proposed approach requires that its projections fall within 5-width tubular neighborhood around u. This is expressed as:
log. (8.)² Ug = 6. [loge (1:
EU Vi=1,.,,
1, measue the cosine value cos(e) between u + ge ATRM (loge (S), u and u, and get the following inequality:
Since the relaxation parameter 0
g AIRM (u)) + ge ATRM (loge (S:). u) u,u)) u + gATRM (log, (S;), u) u/ll.
26, Vi 1i.1,.,,N
which implies a formula for 8 as follows,
E-
5 $
max (1.....
De
ALRM (log,
E2
Since max;;(1,....N) 96 have 5 > 0.
AIRM
(aaaa(oo88 (log. (S:), 2)) is upper bounded as it is finite set, when > 1,
I COMMON EXPERIMENTAL SETTINGS
In all experiments, learn the considered architectures from random initializations and optimize the parameters with ADAM (Kingma & Ba 2017). To propagate gradients through the network standard minibatch-basee backprop learning with extension for structure matrices (Ionescu et al.)
2013) and manifold-constranned gradients (Absil et al. 2009 Becigneul & Ganea 2019) is used.
To implement the considered architectures, either used publicly available python code or imple- mented the methods in python using the packages torch (Paszke al. 2019). (Perregee et al. 2011)
and geoopt (Kochurov et al. 2020
24

---

Published as a conference paper at ICLR 2024
expx(Iu)
(S)
"
0
1
T TOSH
Figure 6: Illustration Relaxation Deviation: The yellow region represents the double cone, where the maximum cosine value of the included angles is The region within the red dashed lines denotes the tubular neighborhood of width d along the speed of geodesic u on the tangent space.
J EXPERIMENTAL SETTINGS FOR MOTOR IMAGERY CLASSIFICATION
J.1 EXPERIMENTAL DATASET
J.1.1 KOREAN UNIVERSITY DATASET
The Korean University (KU) dataset, also known as Lee2019M1 in MOABB contains EEG signals collected from 54 subjects who participated in a binary-class EEG motor imagery task. The EEG signals were recorded at a sampling rate of 1,000 Hz using 62 electrodes. The EEG signals were segmented from the first second to 3.5 seconds relative to the stimulus onset, resulting a total duration of 2.5 seconds for each trial. The KU dataset is the largest dataset for the two-movement publicly available motor imagery dataset.
The dataset was partitioned into two sessions, denoted S₁ and S2, each encompassing training and testing phases. In each phase, 100 trials were evenly allocated between right and left-hand imagery tasks. As result, there were total of 21,600 trials available for evaluation, where 21,600 trials can be expressed as the product of 2 sessions, 54 subjects, and 200 trials per subject.
J.1.2 BNC12015001 DATASET
The BNC12015001 dataset from BNCI Horizon 2020 involved participants tasked with per- forming imagery of right-hand movement versus imagery of movement in both feet. Data were recorded at a sampling rate of 512 Hz and filtered using bandpass filter ranging from 0.5 to 100 Hz, along with notch filter 50 Hz. Recording started at 3.0 seconds following the prompt and continued until the cross period's conclusion at 8.0 seconds, resulting in a total duration of 5 seconds for each trial.
Most of the studies (Subjects 8) were conducted in two sessions, denoted as SA and Sg, which occurred on consecutive days. However, for Subjects 10, 11, and 12, three sessions were denoted as SA. Sg. and Sc. In each session, every participant completed 100 trials for each class, totaling
5 The MOABB package includes benchmark dataset designed for motor imagery decoding algorithms, encompassing open-access datasets and involving over 250 subjects, the following address: https:
The datasets from the Horizon 2020 project available for access at the following website http:
25

---

Published as a conference paper at ICLR 2024
200 trials per participant per session. Therefore, for self- supervuue learning, this dataset provide a total 5,600 trials for training. (5,600 trials = 200 trials/session x 28 sessions).
J.2 EXPERIMENTAL SCENARIOS
There are two primary evaluation scenarios utilized in this study.
The first scenario involves a 10-fold validation approach. This setup divides each subject's data into 10 equally sized and class-balanced folds. Nine folds are used for training, while one is reserved for testing. This process is repeated ten times using a different fold for testing.
The second scenario is referred to as the holdout scenario. We employ the data from the first session for training, the first 100 trials of the second session for validation, and the remaining 100 trials for testing. It's important to note that the EEG data for these two sessions is typically collected on different days, introducing potential variability.
J.3
IMPLEMENT DETAILS
J.3.1 KOREAN UNIVERSITY DATASET
The following two-view DeepGeoSSL model is employed:
1d View: Two electrodes (C3 and C4);
2d View: Twenty electrodes (FC-5/3/1/2/4/6, C-5/3/1///21455, and CP-5/3/1//2/116)
All 21,600 trials are utilized in training the self-supervised learning framework.
For the downstream classification task,
tensor-based approach for segmentation in time and frequency domains (Ju & Guan 2022). In the time domain, there are three segments as follows:
0.0 1.58, 0.5 28, and 1.0 2.58.
The frequency band has been partitioned into 9 segments, each spanning 4Hz bandwidth. These segments cover the range from 4Hz to 40Hz and are evenly distributed within this range. To process the digital signals, Chebyshev Type filters with Hz intervals were employed. These filters were meticulously designed to ensure a maximum passband loss of 3 dB and a minimum stopband attenuation of 30 dB.
It results in 27 channels (3 temporal segments x 9 frequency bands). Additionally, the BiMap layer in the architecture takes the input dimension of 20, transforms 30, and then returns to an output dimension of 20. The architecture details can be found in Figure
J.3.2 BNC12015001 DATASET
The following two-view DeepGeoSSL model is employed:
Vie- Two electrodes (FCz and C3);
View: Twenty electrodes (FC-3/z/4, C-5/3/1/2/2/416, and CP-3/z/4).
All 5,600 trials are utilized in training the self-supervised learning framework.
In the time domain, there are three segments follows:
0.0 ~ 1.0s, 1.0 ~ 2.0s, 2.0 ~ 3.0s, 3.0s ~ 4.0s, and 4.0 ~ 5.0s.
The frequency division was carried out at the same intervals and filter settings as in the KU dataset, resulting in 45 channels (5 temporal segments 9 frequency bands). Additionally, the BiMap layer in the architecture takes the input dimension of 13, transforms to 30, and then returns it to an output dimension of 13. The architecture details can be found in Figure
26

---

Published as a conference paper at ICLR 2024
J.3.3 GEODESIC SPEED
Since the downstream motor imagery tasks emphasize spectral information for classification, uti- lize the largest principal component from Tangent Principle Component Analysis (Tangent PCA) as the initial speed for both v and follows,
V = w = TPCA (PCA(fvec(log)) (uec(log, (ee)l (S.))
where {S:) represents set of 20-channel EEG signals. In particular, we compute v = w using Tangent PCA for each channel of 27 channels for the KU dataset and 45 channels for the BNCI2015001 dataset.
The choice of this speed be understood as both modalities projecting signals onto high- information--ontent spectral principal component directions and keeping the two signals geodesi- cally consistent. This is akin neural networks augmenting lov-informationn-ooeenn signals with sufficient information.
J.3.4 ESTIMATION OF COVARIANCE MATRICES
is well-known that the estimation of the sample covariance matrix can be inaccurate, particularly in nigh-dimenssonal settings with limited observations. To address this challenge, we employ reg- ularization technique known as shrinkage, which involves adding a scaled identity matrix (typically multiplied by small value, such as le-2 in this experiment) to the sample covariance matrix.
J.4 AVERAGE ACCURACY CHANGES
In this experiment, we have two figures that present the average accuracy changes for each par- ticipant across six scenarios in two datasets. The bascline is the result obtained by running the 2-channel EEGs using Tensor-CSPNet. The changes are calculated by subtracting the baseline re- sults from DeepGeoCCA (s 0.05, an 1, Q2 Supplementary Table additionally summarizes the significance of these differences. In these fig- ures, green areas indicate the subject who has improvements, while blue areas represent those with declines. Our method is effective for a large proportion of subjects in all scenarios.
a3 0.1 ). shown in Figures 9 and I0
Table Simultaneous EEG-IMRI Dataset Results. Statistical test results for downstream test set accuracy differences in the multi-view EEG experiment, extending Table We used permutation, paired t-tests to identify significant differences between the lower bound (i.e., no pre-training: channel EEG) with DeepGeoCCA for three scenarios (S1. S2. S₁ + Sz; t-max adjustment for multiple comparisons) and the KU (df=53, le4 permutations) and BNCI datasets (df=11, exhaustive permutations). Significant differences highlighted in bold (p-val < 0.05) and trends in italic (p-val S 0.1).
KU dataset
BNCI dataset
scenario
SSL pretraining
mean (std)
t-val (p-val)
mean (std)
I-val (p-val)
62.97 (12.46)
-6.1 (0.0001)
74.12 (16.50)
-3.1 (0.0391)
DeepGeoCCA
65.90 (12.95)
77.58 (14.56)
S2
65.85 (12.95)
-4.0 (0.0006)
75.62 (14.58)
-2.9 (0.0469)
DeepGeaCCA
68.11 (13.64)
78.46 (13.95)
Si + S2
62,89 (12.29)
-4.4 (0.0003)
75.00 (14.73)
-2.6 (0.0825)
DeepGeoCCA
65.43 (13.42)
77.83 (13.70)
J.5 ABLATION STUDY
The loss function of DeepGeoCCA has total of four hyperparameters: for the e-geodesic con- straints, and three coefficients dr a2 and 0g for the three loss functions. the following two ablation studies, we analyze these parameters:
7 We implemented the initialization of and w using the Tangent PCA approach according to the Geomstats Package in the following address: htdpss//cctsss/eeeetttttt....ri http://wo...ooooooooooooooo////a ia met hothoss ri emanni an_ oaaaaannnn__eeeeeeeeoonnnnneeeen
html
27

---

Published as a conference paper at ICLR 2024
In the first experiment, the overall performance of the model does not vary significantly with differ- ent values of E. Smaller values of E lead to a slight improvement classification performance. The results reported in the main text are based on the setting where 0.05. In the second experiment, we observed that using the variance loss slightly improved the performance compared to not using it. Therefore, in the results reported in the main text, we used the variance loss with coefficient of = 0.1.
De
S
S2
C
S
S2
Lower Bound
62.97 (15.72)
65.85 (16.03)
Lower Bound
62.97 (15.72)
65.85 (16.03)
8 0.9
65.94 (16.10)
68.19 (16.02)
0g =
65.74 (16.59)
67.33 (16.34)
8 0.7
65.82 (15.35)
67.77 (16.04)
0g = 0.1
66.26 (17.05)
68.27 (17.23)
E = 0.5
66.23 (16.10)
68.06 (16.66)
= 0.3
66.14 (15.93)
68.09 (16.07)
0.1
66.03 (15.56)
68.13 (15.87)
0.05
66.26 (17.05)
68.27 (17.23)
66.32 (16.40)
68.13 (16.55)
(a) Le in L (a1, 02 1, 03 0.1).
(b) L. in L (a1, 02 1, 0.05).
Table 6: Ablation Studies: Average (std across subjects in brackets) classification accuracy for the motor imagery task (higher is better). We considered the KU dataset with scenario 10-fold CV (sessions S1 and S2).
Se 52,
S 20
Downstream
Classifier Network
BiMap (2->30)
BiMap (20->30)
BiMap (20 220->30
Dimension Adaptation Network
BIMap (30 -0->20
BiMap (80~>20)
BiMap (30->2)) ~20)
ReEig
ReEig
ReEig
Riemannian BN
Riemannian BN
Riemannian BN
BiMap (20->30)
LOG
BiMap -30->20 220)
Transformation Network
log, (Z) E 200
ReEig
Geodesic Consistency
RieBN
Z E 820
Riemannian BN
Projection Head
BiMap (20 (20=>20)
LOG
Figure 7: Illustration of the network architecture for the KU dataset. The spatial covariance matrices from both the 2-channel EEG setup (i.c., C3 and C4) and the 20-channel EEG setup (i.c., FC- 5/3/1/2/4/6, C-5/3/1////14/5, and CP-5/3/1/z/2/4/6 are transformed into common dimension using the Dimension Adaptation network. During the training of the self-supervised learning (indicated by the solid black lines), the Transformation network and Projection Head retain the parameters in the neural networks, and geodesic correlation is calculated. The dashed red lines represent the downstream task classifier. Specifically, in the upper right corner of Riemannian BN the asterisk denotes that we use specially designed Riemannian BN to perform center- martarl translation of outputs to point Please refer to Section for information on the layers within the architectures.
28

---

Published as a conference paper at ICLR 2024
SE 8++
S E 8
Downstream Classifier Network
BiMap (2->30)
BiMap (13->30)
BiMap (13->30)
Dimension
BiMap (30 13)
BiMap (30->1) 100~>13)
Adaptation Network
BiMap (30->13)
ReEig
ReEig
ReEig
Rlemannian BN
Riemannian BN
Riemannian
BiMap (13->30)
LOG
BiMap (30 101~>13
Transformation Network
log, (Z) E 8
ReEig
Geodesic Consistency
RieBN
Z E 8/
Riemannian BN
Projection Head
BiMap (13 (33~)33)
LOG
Figure 8: Illustration of the network architecture for the BNC12015001 dataset. The spatial covari- ance matrices for the 2-channel EEG setup are derived from FCz and C3, while for the 20-channel EEG setup, they are derived from FC-3/z/4, C-5/3/1/2/214/6, and CP-3/z/4. The description this architecture is consistent with it provided in Figure 7 For detailed explanation, please refer to it.
29

---

Published as a conference paper at ICLR 2024
KU S
KU S2
KU S
Sub.54
+11.5
Sub.20
+14.0
Sub.28
+16.0
Sub.06
+9.0
Sub.46
+14.0
Sub.06
+11.0
Sub.19
+9.0
Sub.02
+8.55
Sub.09
+10.0
Sub.45
+7.5
Sub.35
+8.0
Sub.32
+9.0
+7.5
Sub.18
+7.5
Sub.27
+8,0
Sub.21
+6.0
Sub.33
+7.5
Sub.24
+7.0
+6.0
Sub.51
+7.5
Sub.21
+7.0
Sub.16
+6.0
Sub.26
+6.0
Sub.54
+7700
Sub.25
+5.5
Sub.04
+5.5
Sub.44
+6.0
Sub.05
+5.5
Sub.39
+5.0
Sub.04
+6.0
Sub.44
+5.5
Sub.42
+5.0
Sub.19
+6.0
Sub.28
Sub.44
+5.0
Sub.10
+5.0
Sub.33
+5.0
Sub.37
+4.5
Sub.50
+5.0
Sub.10
+5.0
Sub.01
Sub.22
+5.0
Sub.37
+4.5
Sub.10
+4.5
Sub.07
+4.0
Sub.22
+4.5
Sub.32
+4,0
Sub
+4.0
Sub.43
+4.5
Sub.52
+3.5
Sub.29
+4.0
Sub.09
Sub.13
+3.5
Sub.31
Sub.04
+4.55
Sub.14
+3.5
Sub.17
+4.0
Sub.18
+4.0
Sub.15
Sub.01
Sub.13
+4.0
Sub.16
+3.0
Sub.37
+3.0
Sub.03
+4.0
Sub.22
Sub.52
Sub.32
+4,0
Sub.40
+3.0
Sub.46
+3.0
Sub.31
Sub.28
Sub.
Sub.14
+3.5
Sub.47
+2.5
Sub.39
+3,0
Sub.20
Sub.08
Sub.34
Sub.36
+3.5
Sub.49
+2.5
Sub.11 +2.0
Sub.52
Sub.48
+2.0
Sub.12
+2.0
Sub.30
+3.5
Sub.05
+1.5
Sub.48
+2.0
Sub.49
Sub.25
Sub.13
Sub.01
+3.0
Sub.23
+1.5
Sub.30
+2.0
Sub.27
Sub.43
Sub.42
Sub.07
+3,0
Sub.30 +1.5
Sub.35 +2.0
Sub.12
Sub.36
Sub.38 +2.0
Sub.50
+2.5
Sub.38 +1.0
Sub.45 +100
Sub.29
Sub.45
Sub.41
Sub.51
+2.0
Sub.27 +1.0
Sub.36 +1.0
Sub.24
Sub.29
Sub.18
Sub.08
+2.0
Sub.03 +1.0
Sub.33 +1.0
Sub.15
Sub.21
Sub.51
Sub.47 +1.0
Sub.33 +0.5
Sub.25
Sub.35
Sub.53
Sub.17 +0.5
0.55 Sub.54
Sub.05
Sub.34
Sub.06
-1.0
Sub.08
Sub.02
0.51 Sub.24
-1.0
Sub.26
-1.0
Sub.38
-0.5 Sub.50
-2.0
Sub.20
-1.5
Sub.41
11.0 Sub.19
2..0
Sub.49
-1.5
Sub.23
-1.0
Sub.17
33.0
Sub.12
-1.5
Sub.46
11.5
Sub.31
33.0
Sub.02
2.5
Sub.48
Sub.11
Sub.23
2.5
Sub.40
-2.5
Sub.07
-3.0
Sub.40
-2.5
-3.5
-4.0
-4.0
Sub.42
-8.5
Sub.53
-4.0
Sub.43
-9.5
-9.0
Sub.12
-8.0
Figure 9: This figure displays the average accuracy changes (%) for each subject in the KU dataset
across two 10-fold cross-validation scenarios S1 and S2, and holdout scenario S + Sg. Subjects with larger improvements are positioned at the top, arranged in descending order from top to bottom. Green represents increase, while blue represents decrease. Corresponding increases labeled with the respective subject number.
30

---

Published as a conference paper at ICLR 2024
BNCI2015001 St
BNCI2015001 Sp
BNCI2015001 St
Sub.05
+30.0 Sub.20
SUb.09
949.11
+7.5
Sub.04
6.5
Sub 05
7.0
500.30
1.0
Sub.05
10.0
Sue 10
16.0
940.08
+5.5
9ib.00
+5.
966 12
+5.0
940.09
+5.5
+4.
+10
14.0
43.5
1..0
+2.0
Sub 02
Sub.12
Sub.01
Sub.01
2.5 Sulu2
5,0 3ub.1
20
-2.5
0.0
25
5.0
7.5
Figure 10: This figure displays the average accuracy changes (%) for each subject in BNC12015001 dataset across two 10-fold cross-validation scenarios SA and Sg. and holdout scenario SA + Sg. Subjects with larger improvements are positioned at the top, arranged in descending order from top to bottom. Green represents increase, while blue represents de- crease. Corresponding increases are labeled with the respective subject number.
31

---

Published as a conference paper at ICLR 2024
K SIMULTANEOUS EEG-FMRI DATA
K.1 EXPERIMENTAL DATASET
We used publicly available dataset Dattaserocontaininng simultaneous EEG-AMRI recordings collected from 8 subjects while they were resting in two different conditions (eyes closed or open). This dataset was originally recorded to evaluate EEG artifact correction methods specific EEG data recorded inside MR scanners (van der Meer et al., 2016b). As a result, it contains several experimental settings. For this analysis, we only consider the simultaneous EEG-IMRI setting where EEG and fMRI data were recorded synchronously. In this setting, the dataset contains 2 runs (i.e., continuous recordings lasting 4.5 minutes) per subject. Within each run, subjects were instructed with visual cues (every 30 seconds) to alternate between the conditions.
For details about the experimental protocol, recording devices and settings, please refer to the ac- companying dataset paper (van der Meer et al. (0T66)). In brief, EEG data was recorded with 30 MR compatible EEG electrodes placed at standardized locations to cover the whole head. Their electri- cal activity was sampled at 5 kHz with MR compatible EEG amplifiers. To facilitate EEG artifact correction (Bullock et al. 202T) the EEG amplifiers' clocks were synchronized with the MR scan- ner. The dataset was recorded inside two different MR scanners (4 subjects inside a Siemens TIM Trio and the rest inside a Siemens Verio). In both MR scanners standard echo-planar imaging (EPI) sequences were used record fMRI activity. The acquisition parameters were identical except for a slight difference in the repetition time (TR) 1.95 for TIM Trio and 2.0 for Verio. As a result, the fMRI recordings contain volume (i.c., discrete grid with 3 mm isotropic resolution inside the volume covering the subject's 's head) activity with a sample interval of 2.0s (or 1.95 for TIM Trio).
K.2 DATA PREPROCESSING
The data of each run were preprocessed separately.
EEG data were preprocessed with custom pipeline. The dataset provides raw and MR gradient- artifact (Allen et al., 20000. Níazy et al. 2005 corrected EEG data; we used the latter. The activity of each EEG channel was bandpass filtered (0.5 to 125 Hz), before attenuating pulse artifacts (van der Meer et al., 2016b). After these EEG- GEG-IMRI specific preprocessing steps, we resampled the data at 250 (70 Hz cut-off frequency) and applied further commonly used EEG preprocessing steps. In nutshell, we detected and interpolated bad EEG channels (Jas et al. 2017). rejected transient single- channel artifacts (de Cheveigné 20166), re-refereneed the channels to their common average, attenu- ated line-noise at HZ (De Cheveigné 2020), and rejected transient muscle artifacts (de Cheveigné & Parral 2014). The residual EEG activity was decomposed into independent components (ICs) (Ablin et al. 2018). These ICs were submitted to an automatized model (Pion- Tonachini et al.) 2019) to classify ICs into various artifacts or brain activity. We used tolerant thresholds to maintain on average. 85% (9% std) of the ICs. After rejecting artifact ICs, the residuals were projected back to the original EEG channel space. We excluded 6 temporal and pre-frontal channels, known to cap- ture mostly artifacts, resulting in total of 24 EEG channels. As last step, their data was bandpass filtered from to 36 Hz. and resampled the data at 125 Hz to extract broad-band brain-activity.
fMRI data preprocessing was similar TThmas et al. 2022). First, used fmriprep (Esteban et al. 2019) to perform standard anatomical and functional preprocessing steps on the raw anatomi- cal and functional scans. In a nutshell, fmriprep performed spatial normalization of fMRI data to applied slice-time correction and estimated various con- founds. As confound regressors, we selected global signal, motion parameters basic parameters) and the first 5 principal components in white matter (WM) and cerebrospinal fluid (SCF). Using the nilearn package ? we regressed out these confounds from the fMRI data, applied spatial smooth- ing (Gaussian kernel; FWHM at 3 mm), extracted activity of 128 regions of interest (DiFuMo atlas (Dadi et al. (020)), applied a bandpass-fibl (0.01 Hz to 0.2 Hz), and resampled the output at Hz (i.e.. 1 volume per second).
standard space
http://wwlllgggeleee ferr f
httpa://nileamm.giin
32

---

Published as a conference paper at ICLR 2024
EEG and [MRI data were segmented into epochs (20s duration) using triggers marking the onsets of eyes open and closed conditions. As suggested by Défossez et al. (202)), we used robust z-scoring to standardize the EEG channel and fMRI ROI data per and clip resulting activity beyond the interval [-20, 20) to reduce effects of gross outliers. The paired data within these epochs were then distributed among train, validation and test sets. We considered 10-fold stratified, leave- one-run-out (LORO) and leave-one-subjeettour (LOSO) cross-validation scenarios, as detailed in Figure To boost the number of observations in each set, we extracted sliding windows (10s duration, 9 overlap) per epoch.
10-fold spilt the 128 i- oudjects x runs x epochas epocns randomiy into train and leot swt with stroniication for subject run.
Scenarios (example splits)
Extraction of sliding windows
(10s duration)
LOSO
Dataset Structure
Wifrodowe
Supject
Validation
Truin
Thain
Fi
IIfS vildution
Epoch 1
(duration 20s)
Bubject
a
Figure 11: EEG-IMRI dataset. Visualization of the cross- vals-vaaii (CV) scenarios and sliding windows extraction. (left) The considered dataset contains runs for each subject. Each run contains 8 trials with 30 duration (without breaks in-between) during which the subjects with eyes open (EO) or eyes closed (EC). We extract epoch per trial (20s duration, 2 offset). (middle) The epochs fold stratified (10-fold), leave-one-run-out (LORO), and leave-one-subject-ouu CV. (right). After forming the CV splits, sliding windows were extracted from each epoch.
either resting
split across CV folds depending on the scenario. We considered 10-
K.3 IMPLEMENTATION DETAILS
The implementation details provided in this section apply to the simulation and EEG-MMRI experi- ments.
Architecture Figures and visualize the architecture used in the simulation and EEG-IMRI experiments.
Geodesic Speeds We fitted the geodesic speeds in data-driven fashion. Specifically, we treated them as learnable parameters on the Stiefel manifold (random initializations). In addition to L. we additionally employed them to project the latent tangent space representations (log, (Z:), log, (2)) scalars and quantify correlation among these in Lp-
Parameter Optimization We used a minibatch based training scheme to fit the learnable param- eters for 250 epochs. Batches contained 128 paired views (stratified for subject and run). As op- timizer, we used ADAM with learning rate of 1e-3 (5e-3 for the simulations experiment) and applied a weight decay 1e-3 to all learnable parameters except the manifold-conssrrined ones. To put more emphasis
maximizing correlation than fulfilling the geodesic projection constraint, we set a1 = 1, 02 0.25 and ag = O. To discourage the architectures to collapse dimensions, fixed the learnable rescaling parameter of SPDMBN to 1 so that the Fréchet variance (6) of the latent representations remains approx. constant.
During training, used dedicated validation set to monitor the total loss held out data after ery epoch. The parameters that resulted in the lowest validation set loss, were then used to transform the test data.
33

---

Published as a conference paper at ICLR 2024
S 8
5 E SES+
Dimension Adaptation
BiMap (2->2)
BiMap (2 (2-22)
Transformation Network
ReEig(1.0)
ReEig(1.0)
Z 8
Ž E 52
Tangent Space Projection
SPDMBN(e.1.0)
SPDMBN(e.1.0)
LOG
LOG
Geodesic
(Z) E 82
log, (2) 82
Correlation
Figure 12: Illustration of Network Architecture for the simulations with paired SPD manifold-valued data. The BiMap layers used unconstrained parameters. For the ReiEig layers, we used a threshold of 1.0. The tangent space projection module used a SPDMBN layer without learnable parameters, i.e. the rebias parameter was fixed at the identity matrix and the scaling parameter was fixed to 1..
$ 24
OR
S 1250x24
Š R 10x128
Dimension Adaptation & Transformation Network
BiMap (24->10)
Transfomer-Encodan layer
TSMNet temp filt=4 40
p_drop=0.5
ReEig(0.55)
lat dim=64
Z 810
2 E R 55
Tangent Space Projection
S.DDSMBNic.1.0)
Batch Norm(0.1100
LOG
upper
Geodesic Correlation
(Zg E 100
log,(2) 10
Figure 13: Illustration of Network Architecture for the EEG-MMI experiment. For EEG views S. we used either SPDNet or TSMNet. For the ReiEig layer within SPDNet, we used a threshold of 0.5 and unconstrained parameters in BiMap. For TSMNet. used standard parameters except for 10 latent SPD dimensions instead of 20. The tangent space projection module used a domain- specific SPDMBN layer without learnable parameters. For fMRI views S. use standard single Sunsformer-encodee layer as implemented in torch.
34

---

Published as a conference paper at ICLR 2024
K.4 SUPPLEMENTARY RESULTS
In addition to the 10-fold CV results reported in the main text, this section presents additional results that test generalization across runs (LORO CV scenario) and subjects (LOSO CV scenario). The results are summarized in Table[]and and Figure along with the 10-fold CV results provided in the main text. Comparing the results, we noticed large variability for the LORO and LOSO evaluation scenarios. Particularly for subjects 1 and 5 (highlighted with red dots in Figure 14, all models failed to generalize. To understand why this is the case, we consulted the accompanying publication van der Meer et al. 6b) and identified that for these subjects, there was task effect (eyes open vs. closed) in the EEG data. Since the effect is prominent in the other subjects, naturally, models that utilize the "majority" effect will not generalize to these particular subjects. That is why report the results in Table with (w/) and without (w/o) test-set folds associated with subjects 1 and S. Apart from the large variability (mostly due to subjects I and 5), we found that the overall effects identified in the 10-fold CV scenario are preserved in the other evaluation scenarios. If we exclude the results for subjects 1 and 5, the average correlation score for all methods increases drastically (approx. 0.2). Particularly, the combination of DeepGeoCCA and TSMNet yields correlations close to 0.7 on held-out data. These additional results clearly demonstrate that our proposed framework can extract representations that generalize to new subjects, and pave the way for potential clinical applications with sufficiently large and diverse datasets.
Table 7: Simultaneous EEG-MMRI Dataset Results. Additional test-set results for the EEG-MMR dataset that extend Table[] with regard to the considered evaluation scenario. In addition to the 10-fold CV scenario results, reported in Table[]] this table also summarizes results for leave-one- run-out (LORO) and saxe--one-uujecc--on (LOSO) CV scenarios. For LORO and LOSO, we report summary statistics (mean and std) with (w/) and without (w/o) outlier subjects (i.e., subjects 1 and 5).
Evaluation Outlier
10-Fold CV
LORO w/o
LOSO w/o
EEG model
w/
Method
Correlation
X
CCA
0.02 (0.15)
0.18 (0.35)
0.34 (0.21)
0.13 (0.30)
0.27 (0.17)
RieCCA
0.38 (0.10)
0.37 (0.44)
0.59 (0.13)
0.24 (0.39)
0.45 (0.10)
Tensor- CSPNet
DeepCCA
0.42 (0.11)
0.36 (0.43)
0.56 (0.23)
0.30 (0.32)
0.45 (0.16)
VICReg
0.29 (0.19)
0.24 (0.43)
0.46 (0.21)
0.31 (0.28)
0.43 (0.19)
DeepGeeCCA
0.44 (0.11)
0.39 (0.45)
0.58 (0.27)
0.24 (0.47)
0.44 (0.35)
TSMNet
DeepCCA
0.29 (0.10)
0.17 (0.34)
0.25 (0.32)
0.12 (0.46)
0.30 (0.35)
VICReg
0.38 (0.10)
0,42 (0.36)
0,58 (0.22)
0.25 (0.37)
0.42 (0.19)
DeepGeoCCA
0.58 (0.08)
0.46 (0.45)
0.68 (0.18)
0.45 (0.43)
0.68 (0.14)
R2 f
X
CCA
0.07 (0.08)
-2.22 (1.29)
-2.26 (1.31)
-1.05 (0.70)
-1.19 (0.76)
RieCCA
0.01 (0.01)
0.00 (0.02)
-0.00 (0.02)
0.01 (0.01)
0.01 (0.01)
Tensor- CSPNet
DeepCCA
0.30 (0.07)
0.27 (0.13)
0.28 (0.13)
0.23 (0.10)
0.24 (0.12)
VICReg
0.28 (0.10)
0.18 (0.12)
0.20 (0.14)
0.23 (0.07)
0.24 (0.07)
DeepGeoCCA
0.58 (0.02)
0.56 (0.05)
0.57 (0.05)
0.54 (0.07)
0.54 (0.08)
TSMNet
DeepCCA
0.18 (0.09)
0.10 (0.06)
0.11 (0.05)
0.13 (0.04)
0.13 (0.05)
VICReg
0.24 (0.03)
0.22 (0.11)
0.23 (0.10)
0.22 (0.09)
0.22 (0.09)
DeepGeoCCA
0.55 (0.02)
0.53 (0.06)
0.54 (0.06)
0.54 (0.02)
0.54 (0.02)
L SPECIFIC CHALLENGES ASSOCIATED WITH NEUROIMAGING DATA
For the neuroimaging tasks that we study, we have to deal with relatively small datasets (recording human neuroimaging data is expensive and time- consuming,, non-stationary signals and distribution shifts that limit generalization across days and subjects (also indicated in Figure 4), inherently low signal-to-noke -o-ie ratios, transient high-variance artifacts (i.e., outliers), and non-line lineir, poorly under- stood coupling mechanisms that link views of different modalities (Philiastides al. 2021). For
35

---

Published as a conference paper at ICLR 2024
LORO
Figure 14: Graphical representation of the results summarized in Table 7 for the correlation metric on test data. Each dot summarizes a CV split result. (left) 10-fold CV scenario. (middle) LORO scenario (16 runs=folds). Red dots highlight fold for which the test set was sampled from outlier subjects (subjects and 5). (right) LOSO scenario (8 subjects folds).
more details, let refer to recent reviews elaborating not only challenges but also opportunities of EEG BCI (Fairclough & Lotte 2020) and simultaneous EEG-MMRI (Warbrick) 2022).
36


Reviews:
[{"review_id": "lTT9mjScRC", "reviewer": "Reviewer_Mrh1", "summary": "The manuscript presents a CCA method, called DeepGeoCCA, that can be applied to SPD-valued data.  The target application is to relate the covariance matrices (encoding brain connectome) from multiple modalities of neuroimaging data. The method first maps the covariance matrices from two-domains (may be of different sizes, e.g., from EEG and fMRI) to a common SPD cone, where two geodesics are found such that the projections onto the geodesics achieve maximum correlation in the tangent space.\n\nThe method was evaluated on a synthetic dataset and on a real dataset of paired EEG-fMRI data. Results show that the method achieved the highest correlation in a cross-validation setting compared to similar non-linear CCA methods.", "questions": "1. Is $e$ the identify or the Frechet mean? The text suggests both.\n\n2. There are two \"Appendix ??\"\n\n3. What is the definition of \"double cone\"?\n\n4. How do we determine the $\\alpha$ weighting parameters?\n\n5. \"DeepGeoCCA shared the highest correlation score with Deep-CCA among the deep learning methods (left panel) at the cost of small R2 (right panel)\" I thought high correlation automatically means high R2?\n\n6. 10s sliding window is really short for a 2s TR. Are there references for this choice?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. CCA analysis on connectome data, brain connectivity matrices, has a great potential in advancing neuroimaging research. The target application can be very broad, e.g., to relate fMRI with structural connectome derived from Diffusion Tensor Imaging.\n\n2. The method can relate covariance matrices of different sizes, which increases the flexibility and can be applied to even broader scenarios, e.g., relating connectome data from different brain parcellations.", "weaknesses": "1. The motivation of using the geodesic relaxation (section 3.2) in intuitive a CCA setting. Why do the data points have to lie near the geodesics? In normal CCA or even deep CCA, we never have such constraints to encourage data points to align with the components. Component analysis is about finding directions that explain the data, not the other way around.\n\n2. My biggest concern is the cross-validation setting. These are small datasets (8, 54, and 12 subjects in three datasets). How do you perform 10-fold CV on 8 subjects? I would then assume that the data split is not on the subject level, but on the matrix level (each subject has many covariance matrices generated by a sliding window). These matrices are derived from overlapping time series and have high dependency. Having subject-level split is essential to make the results clinically relevant, otherwise one would just need to find a method that can overfit to the training data.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "cjEEzyUPWj", "reviewer": "Reviewer_7bLk", "summary": "This paper introduces a novel measure called geodesic correlation and proposes an innovative geometric deep learning framework called DeepGeodesic Canonical Correlation Analysis (DeepGeoCCA) to enhance the consistency of latent representations obtained from paired views in covariance-based neuroimaging data. The authors demonstrate the effectiveness of DeepGeoCCA on both EEG-fMRI and fMRI-fMRI datasets, achieving state-of-the-art performance on several benchmark tasks. The paper's contributions include the introduction of geodesic correlation, the development of DeepGeoCCA, and the demonstration of its effectiveness on various neuroimaging datasets.", "questions": "The authors demonstrate the effectiveness of DeepGeoCCA on both EEG-fMRI and multiview EEG datasets, but it is not clear how the proposed framework would perform on other types of neuroimaging data, such as diffusion tensor imaging (DTI) or positron emission tomography (PET) data. Can the authors comment on the potential applicability of DeepGeoCCA to these types of data?", "limitations": "", "rating": 6, "confidence": 1, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper presents a novel measure called geodesic correlation and proposes an innovative geometric deep learning framework called DeepGeoCCA to enhance the consistency of latent representations obtained from paired views in covariance-based neuroimaging data. The paper is well-written and clearly explains the concepts and methods used in the study. \n\nOriginality: The introduction of geodesic correlation is a novel measure that expands traditional correlation consistency to covariance-based data on the symmetric positive definite (SPD) manifold. The proposed DeepGeoCCA framework is also highly original, as it is the first completely data-driven approach to learn latent representations from paired EEG-fMRI data. The effectiveness of DeepGeoCCA is demonstrated on both EEG-fMRI and fMRI-fMRI datasets, achieving state-of-the-art performance on several benchmark tasks.\n\nQuality: The paper is clearly presented and well-structured. With detailed descriptions of the datasets, experimental setups, and evaluation metrics used in the study, it is easy for other researchers to reproduce their results.\n\nClarity: The paper is well-written and the figures are clearly plotted to illustrate the concepts and methods used in the study. \n\nSignificance: The paper is highly significant, as it presents a novel measure and framework for enhancing the consistency of latent representations obtained from paired views in covariance-based neuroimaging data. The proposed DeepGeoCCA framework has the potential to improve our understanding of the dynamic relationships among brain networks and their alterations in pathology, which could lead to improved diagnosis and treatment of neurological disorders. The framework is also not limited to neuroimaging data and can be applied to other types of data with geometric structures, making it a valuable contribution to the field of deep learning.", "weaknesses": "The clarity of presentation needs to be improved, with typos scattered throughout the paper:\nIn Section 2.1, I would suggest using the notation $S^n$ (instead of Sym_n) to denote the space of n × n real symmetric matrices, to be consistent with the the notation S^n_++ for the space of real n × n symmetric and positive definite matrices. \n\n\"For a more comprehensive proof, please consult Appendix ??.\"\n\"In addition, we have use(d) the relaxed\"\nSection 3.3, in the definition of $L_epsilon$ formula, there is a redundant right parenthesis \")\".", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "gkFFp2WqEU", "reviewer": "Reviewer_2nXA", "summary": "While researchers began to utilize multi-model data for neuroimaging analysis with the hope to combine distinct information from multiple modalities, the challenge in learning a shared representation space has arisen. The authors pointed out that this is especially nontrivial when dealing with covariance-based data characterized by their geometric structure.\n\nTo that end, the authors proposed a metric termed “geodesic correlation” that expands traditional correlation consistency to covariance-based data on symmetric positive definite manifolds and measures the consistency of latent representations obtained from paired views/modalities. They further proposed a deep learning framework termed “DeepGeoCCA” for data on symmetric positive definite manifolds that optimizes the geodesic correlation of unlabeled paired data.", "questions": "1. Excuse my poor knowledge, but it seems to me that the term “covariance-based data” is not a standard or widely recognized term in the field of statistics or data analysis. Would it be better to define the term somewhere in the paper?\n\n2. On page 4 and 5, you have LaTeX failures (Appendix ??).\n\n3. Possible typo on page 6, “the activity of a latent, shared source, we generated paired views in $S_{++}^2$ $^2$”.", "limitations": "", "rating": 8, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The motivation to learn a joint latent space with coherence between two latent representations from different modalities is well grounded. \n\n2. Compared to existing methods like DeepCCA and VICReg, both of which learn the correlation relationship between two latent representations, the proposed DeepGeoCCA better retain the geometric relations as analyzed theoretically and empirically.\n\n3. Figures 3 and 4 are well designed.", "weaknesses": "1. The authors proposed a metric and a learning framework for \"covariance-based data\" on symmetric positive definite manifolds, but they have not yet motivated why the EEG/fMRI data are “covariance-based data” and/or why they lie on symmetric positive definite manifolds. While this might be obvious to some people, it is a bit confusing to me.\n\n2. The experiments are a bit weak. While the proposed DeeGeoCCA managed to achieve a higher geodesic correlation, empirical evidence on how this may benefit the learned representation is not clear. Table 2 only showed marginal improvement on the downstream task on relatively small datasets.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (spotlight)