Bernard Ng
email: bernardyng@gmail.com

Martin Dressler

Gael Varoquaux

Jean Baptiste Poline

Michael Greicius

Bertrand Thirion

Transport on Riemannian Manifold for Functional Connectivity-based Classification

Keywords: Classification, Connectivity, fMRI, Riemannian manifold

We present a Riemannian approach for classifying fMRI connectivity patterns before and after intervention in longitudinal studies. A fundamental difficulty with using connectivity as features is that covariance matrices live on the positive semi-definite cone, which renders their elements inter-related. The implicit independent feature assumption in most classifier learning algorithms is thus violated. In this paper, we propose a matrix whitening transport for projecting the covariance estimates onto a common tangent space to reduce the statistical dependencies between their elements. We show on real data that our approach provides significantly higher classification accuracy than directly using Pearson's correlation. We further propose a non-parametric scheme for identifying significantly discriminative connections from classifier weights. Using this scheme, a number of neuroanatomically meaningful connections are found, whereas no significant connections are detected with pure permutation testing.

Introduction

Recent years have witnessed significant growth in the application of resting state functional magnetic resonance imaging (RS-fMRI) for neuroscience research [START_REF] Castellanos | Clinical Applications of the Functional Connectome[END_REF].

Besides addressing the fundamental question of how the brain is wired, substantial effort has been placed on translating RS-fMRI for classification tasks, such as separating diseased and healthy subjects [START_REF] Castellanos | Clinical Applications of the Functional Connectome[END_REF]. Recently, there is a growing interest in using RS-fMRI to examine the effects of intervention on intrinsic connectivity [START_REF] Castellanos | Clinical Applications of the Functional Connectome[END_REF]. For example, comparing the functional connectome before and after motor learning can provide insights into neuroplasticity [START_REF] Castellanos | Clinical Applications of the Functional Connectome[END_REF]. In this work, we focus on the task of classifying connectivity patterns before and after intervention in longitudinal studies. Directly using conventional Pearson's correlation, i.e. normalized covariance, as features for classifier learning has a fundamental limitation. Since covariance matrices, Σ, live on the space of positive semidefinite cone, i.e. v T Σv ≥ 0 for all vectors v, elements of Σ are inherently inter-related [START_REF] Varoquaux | Detection of Brain Functional-connectivity Difference in Post-stroke Patients Using Group-level Covariance Modeling[END_REF][START_REF] Pennec | A Riemannian Framework for Tensor Computing[END_REF], which violates the independent feature assumption implicit in most classifier learning algorithms. In particular, many classifier learning algorithms become unstable in the face of correlated features, i.e. small perturbations on the training data can alter the relative weighting of the features [START_REF] Tolosi | Classification with Correlated Features: Unreliability of Feature Ranking and Solutions[END_REF]. Worsening the situation is the lack of training samples. Connectivitybased classifier learning is usually performed on data from multiple subjects to keep scan time reasonable for each subject, i.e. reliable covariance estimation requires at least several minutes of RS-fMRI data. For typical studies, the number of subjects is no more than thirty or so. Hence, the sample size is rather limited. As a result, the estimated classifier weights would have high variance. Combined with the effects of correlated features, the generalizability of the classifiers is limited, and identification of significantly discriminative connections from these classifiers is nontrivial.

In this paper, we present a Riemannian approach for classifying connectivity patterns before and after intervention in longitudinal studies. The heart of our approach is to treat the positive definite cone as a Riemannian manifold and project the covariance estimates onto a common tangent space of this manifold. On the tangent space, elements of the covariance matrices are no longer linked by the positive definite constraint [START_REF] Pennec | A Riemannian Framework for Tensor Computing[END_REF]. Thus, the impact of correlated connectivity features on classifier learning is alleviated. To bring the covariance estimates of all subjects to a common tangent space, we propose a matrix whitening transport. The underlying idea is to find a subject-specific base covariance estimate that is close to all covariance estimates of each subject and use it for matrix whitening. The resulting covariance estimates of all subjects would be close to the identity matrix. Thus, we can use the tangent space at the identity matrix as the common space for projection. For comparison, we examine the concept of parallel transport from differential geometry [START_REF] Hauberg | Unscented Kalman Filtering on Riemannian Manifolds[END_REF]. To perform parallel transport, we use the Schild's ladder algorithm [START_REF] Hauberg | Unscented Kalman Filtering on Riemannian Manifolds[END_REF][START_REF] Lorenzi | Schild's Ladder for the Parallel Transport of Deformations in Time Series of Images[END_REF]. For validation, we apply our approach to RS scans of fifty one subjects before and after a memory task, and compare it against directly using Pearson's correlation and its regularized variants as features. Furthermore, we propose a non-parametric scheme that combines bootstrapping and permutation testing for identifying significantly discriminative connections from the classifier weights. Detection of neuroanatomically relevant connections is shown.

Methods

Treating the positive definite cone as a Riemannian manifold and using the associated manifold operations (Section 2.1), we propose a matrix whitening transport for covariance projection (Section 2.2), and compare it against parallel transport with Schild's ladder (Section 2.3). Various methods for estimating base covariance matrices are examined (Section 2.4). A non-parametric scheme is devised for identifying significantly discriminative brain connections from classifier weights (Section 2.5).

Manifold Operations for the Space of Positive Definite Matrices

Let + + ∈ S A be a d×d matrix that lives on the positive definite cone, S ++ . Due to the positive definite constraint: v T Av > 0 for all vectors v, elements of A are inter-related.

One way to decouple the elements of A is to treat S ++ as a Riemannian manifold and project A onto the tangent space at a d×d base point,

+ + ∈ S B
, using the Log map [START_REF] Pennec | A Riemannian Framework for Tensor Computing[END_REF]:

2 / 1 2 / 1 2 / 1 2 / 1 ) ( ) ( B AB B B A B - - = logm Log , (1)
where logm(•) denotes matrix logarithm and Log B (A) is the tangent vector at B "pointing towards" A with B assumed to be close to A [START_REF] Hauberg | Unscented Kalman Filtering on Riemannian Manifolds[END_REF]. The key property of Log B (A) for our purposes is that its elements are not linked by the positive definite constraint [START_REF] Pennec | A Riemannian Framework for Tensor Computing[END_REF], which helps alleviate impact of correlated features on classifier learning. More generally, a major difficulty with working in S ++ is that the resultant from even standard operations, such as subtraction, may not reside in S ++ . An elegant solution to this problem is to operate in the tangent space and project the resultant back onto S ++ using the inverse mapping, i.e. the Exp map (2), which guarantees positive definiteness [START_REF] Pennec | A Riemannian Framework for Tensor Computing[END_REF].

2 / 1 2 / 1 2 / 1 2 / 1 ) ( ) ( B AB B B A B - - = expm Exp , (2)
where expm(•) denotes matrix exponential. By combining ( 1) and ( 2), one can compute the geodesic, i.e. local shortest path, between two positive definite matrices [START_REF] Pennec | A Riemannian Framework for Tensor Computing[END_REF]:

] 1 , 0 [ )), ( ( ) ( ∈ ⋅ = t Log t Exp t A γ B B . (3)
This concept of geodesic will be important for parallel transport (Section 2.3).

Proposed Matrix Whitening Transport

Let

+ + ∈ S s c ) , ( C
be the d×d covariance matrix of subject s associated with scan c, where d is the number of brain regions. Also, let

+ + ∈ S s) ( C
be a d×d base covariance matrix of each subject s that is close to C(c, s) for all c. We note here again that our goal is to classify connectivity patterns derived from scans with interventions interleaved in between. Now, if we simply apply (1) to project C(c, s) onto the tangent space at ) (s C , all subjects' projected covariance matrices will lie in different tangent spaces (i.e.

) (s C

are different across subjects), hence not comparable with each other. Instead, under the assumption that

) (s C is close to all C(c, s) of subject s, we can use ) (s C
to approximately whiten C(c, s) for each subject, so that the result-

ing covariance matrices, 2 / 1 2 / 1 ) ( ) , ( ) ( - - s s c s C C C
, would be close to the identity matrix, I d×d , for all subjects. Thus, the tangent space at I d×d can serve as the common space across subjects for projection, which reduces to taking the matrix logarithm:

) ) ( ) , ( ) ( ( ) , ( 2 / 1 2 / 1 - - = s s c s logm s c d C C C C , (4) since
) ( ) ( ) ( 2 / 1 2 / 1 2 / 1 2 / 1 A I AI I I A I logm logm Log d d d d d d d d d d = = × - × - × × ×
. The proposed matrix whitening transport is summarized in Fig. 1(a). In effect, we are isolating the distinctive connectivity attributes by removing commonalities between scans in a nonlinear fashion, and using the residual for classification. It turns out that ( 4) is analogous to the manifold operation deployed in [START_REF] Varoquaux | Detection of Brain Functional-connectivity Difference in Post-stroke Patients Using Group-level Covariance Modeling[END_REF] for building one-class classifier to discriminate subject types, but we are proposing a new perspective on the entailed matrix multiplication as a transport to the neighborhood of I d×d , which justifies why projection onto the tangent space at I d×d is legitimate. To generate well-conditioned, positive definite covariance estimates, we employ oracle approximating shrinkage (OAS), which has a closed-form solution for determining the optimal level of l 2 regularization [START_REF] Chen | Shrinkage Algorithms for MMSE Covariance Estimation[END_REF].

To examine the necessity of matrix whitening, we compare (4) with the case of setting dC(c, s) to logm(C(c, s)), i.e. no transport before projection, which is exactly the Log-Euclidean approach employed by the diffusion MRI community [START_REF] Arsigny | Fast and Simple Calculus on Tensors in the Log-Euclidean Framework[END_REF]. Another simplification to (4), which we refer to as the Euclidean approximation, is to remove the commonalities between the C(c, s)'s in a linear fashion and apply a first order approximation to the matrix logarithm, i.e. logm(A) ≈ A -I d×d , resulting in dC(s, c) = C(c, s)

) (s C -
for all off-diagonal elements. Since linear subtraction can result in non-positive definite matrices, the first order logm approximation is necessary. Also, for connectivity-based classification, only the off-diagonal elements are of interest.

Parallel Transport with Schild's Ladder

Another way of transporting the covariance matrices of all subjects to a common tangent space is to employ parallel transport, which provides the least deforming way of moving geometric objects along a curve on a manifold [START_REF] Hauberg | Unscented Kalman Filtering on Riemannian Manifolds[END_REF]. One way of performing parallel transport is to use the Schild's ladder algorithm [START_REF] Hauberg | Unscented Kalman Filtering on Riemannian Manifolds[END_REF][START_REF] Lorenzi | Schild's Ladder for the Parallel Transport of Deformations in Time Series of Images[END_REF], as summarized in Fig.

C PT1 (c, s) = ))) , ( ( 2 ( 1
) ( ) ( s c Log Exp s s M C C ⋅ .
This procedure is repeated for all G i (s), i = {1, …, g} until we reach R, where we perform projection to find the parallel transported T(c, s): dC(c, s) = Log R (C PTg (c, s)).

In effect, the Schild's ladder algorithm is parallel transporting T(c, s) by forming parallelograms on S ++ but with geodesics in place of straight lines. To enable direct comparisons with the proposed transport, we set R to I d×d . As for g, we have tried values from 1 to 10 with close to exactly the same classification results obtained.

Base Covariance Estimation

Let X(c, s) be a t×d matrix containing d brain region time courses of subject s associated with scan c = 1 to N. For estimating a base covariance matrix that is close to all covariance estimates of each subject s, we examine three mean covariance estimation methods: Euclidean mean, Log-Euclidean mean, and time course concatenation. We exclude the Frechet mean [START_REF] Pennec | A Riemannian Framework for Tensor Computing[END_REF] due to the observed numerical instability, e.g. using Euclidean mean versus Log-Euclidean mean of the covariance estimates as initialization result in different Frechet mean estimates. The Euclidean mean is simply given by: Σ c C(c, s)/N, but does not preserve the spectral properties of the covariance estimates [START_REF] Arsigny | Fast and Simple Calculus on Tensors in the Log-Euclidean Framework[END_REF]. One way to mitigate this problem is to first apply matrix logarithm, take the mean, and apply matrix exponential to bring the mean back to S ++ , i.e. expm(Σ c logm(C(c, s))/N) [START_REF] Arsigny | Fast and Simple Calculus on Tensors in the Log-Euclidean Framework[END_REF]. Yet, another way for estimating a mean covariance matrix is to concatenate X(c, s) across scans into a ct×d matrix and apply OAS [START_REF] Chen | Shrinkage Algorithms for MMSE Covariance Estimation[END_REF].

Discriminative Brain Connection Identification

Critical to neuroscience studies is result interpretability. For identifying significantly discriminative connections from classifier weights, we propose here a non-parametric scheme that combines bootstrapping with permutation testing. Bootstrapping enables identification of the more stable discriminative features, while permutation testing facilitates generation of a null distribution. Importantly, the chance of assigning large weights to the same brain connections for different bootstrap samples would presumably be lower with labels permuted. This intuition is exploited in our proposed scheme, which proceeds as follows. Let w pq be the classifier weights for scan p versus scan q learned from all subjects' dC(c, s) for c = p and q. We first randomly permute the scan labels p and q 10,000 times. For each permutation, we perform classifier learning on each of B = 500 bootstrap samples (with replacement). Denoting the classifier weights for each bootstrap sample b as w b pq , we compute the normalized mean over bootstrap samples: pq w = 1/B • Σ b w b pq /std(w b pq ), and store the maximum element of pq w for each permutation. We then compute pq w without label permutation, and declare its elements as statistically significant if they are greater than the 95th percentile of the maximum normalized mean distribution, corresponding to a p-value threshold of 0.05. We highlight that using maximum statistics implicitly accounts for multiple comparisons [START_REF] Nichols | Controlling the Familywise Error Rate in Functional Neuroimaging: a Comparative Review[END_REF]. The same procedure is applied for finding significantly negative elements of w pq but with maximum replaced by minimum.

Materials

RS-fMRI data were collected from 51 healthy subjects. All subjects underwent two scans of 8 min each. Between the scans, the subjects were asked to memorize 72 words in 12 min and were tested immediately after encoding in a recognition task of 12 min. Acquisition was performed a 3.0T GE scanner with TR = 2.5 s, TE = 30 ms, and flip angle = 90 o . The RS-fMRI data were motion corrected using FSL and spatially normalized using SPM8. Motion artifacts, white matter and cerebrospinal fluid confounds, and principal components from high variance voxels extracted using CompCor [START_REF] Behzadi | A Component based Noise Correction Method (CompCor) for BOLD and Perfusion based fMRI[END_REF] were regressed out from the voxel time courses. A bandpass filter with cutoff frequencies at 0.01 and 0.1 Hz was subsequently applied. For defining brain regions, we employed the atlas in [START_REF] Shirer | Decoding Subjectdriven Cognitive States with Whole-brain Connectivity Patterns[END_REF], which comprises 90 functional regions that span 14 large-scale networks. We restricted our analysis to 78 of these regions due to incomplete cerebellum coverage. Gray matter voxel time courses within each region were averaged to generate brain region time courses. These time courses were normalized by subtracting the mean and dividing by the standard deviation.

Results and Discussion

To evaluate the gain of decoupling connectivity features, we compared Pearson's correlation and its OAS regularized variant against using projected covariance estimates for classifying the connectivity patterns before and after the described memory task. We used l 2 -reguralized linear SVM for classifier learning with the soft margin parameter left at its default value of 1. We opted to use SVM since l1-regularized classifiers tend to be unstable with limited samples [START_REF] Tolosi | Classification with Correlated Features: Unreliability of Feature Ranking and Solutions[END_REF]. For estimating classification accuracy, we performed repeated subsampling over 10,000 random splits: 34 subjects used for training and 17 subjects used for testing. Subsampling on subjects prevents dC(c, s) of different c from the same subject to be used for training and testing, which avoids introducing correlations between the training and test samples. Both Pearson's correlation and its OAS regularized variant resulted in an accuracy of 76%. Projection using the Log Euclidean approach [START_REF] Arsigny | Fast and Simple Calculus on Tensors in the Log-Euclidean Framework[END_REF] without transport obtained an accuracy of 80%. Using the Euclidean approximation to the matrix whitening transport obtained accuracy of 87%-92%. Accounting for the manifold structure of S ++ using the matrix whitening transport and parallel transport achieved accuracy of 98%, which is a 22% increase in accuracy compared to the standard approach of directly using Pearson's correlation as features. Also, the substantial improvement compared to the Log Euclidean approach suggests the necessity for transporting prior to projection. Further, lower variability in accuracy was observed compared to the Euclidean approximation. Overall, our results show that treating S ++ as a Riemannian manifold and using the associated operations to decouple connectivity features is highly beneficial for connectivity-based classification. Since the matrix whitening and parallel transport provided similar performance, from a computational standpoint, the proposed transport is preferred considering the increasing temporal resolution offered by recent acquisition protocols that enable larger connectivity matrices to be reliably estimated. We note that using sparse Gaussian graphical model for covariance estimation resulted in close to exactly the same classification accuracy.

To identify significantly discriminative brain connections, we applied our proposed non-parametric scheme (Section 2.5) on classifier weights learned with all subjects' dC(s, c) as features. We present here only results for dC(s, c) estimated using the matrix whitening transport with mean covariance computed from concatenated time courses. Statistical significance is declared at a p-value threshold of 0.05 with multiple comparisons implicitly corrected for using maximum statistics [START_REF] Nichols | Controlling the Familywise Error Rate in Functional Neuroimaging: a Comparative Review[END_REF]. We note that no significant connections were found with pure permutation testing on classifier weights derived from the same dC(s, c), hence illustrating superior detection sensitivity with our proposed scheme for discriminative connection identification.

Using our non-parametric scheme, a number of significant connections between memory-relevant regions were found (Fig. 2(b)). The posterior cingulate cortex (PCC), precuneus, angular gyrus, and medial prefrontal cortex (mPFC) are key regions of the default-mode network (DMN) [START_REF] Greicius | Functional Connectivity in the Resting Brain: A Network Analysis of the Default Mode Hypothesis[END_REF], which are closely associated with memory processes [START_REF] Kim | Dissociating the roles of the default-mode, dorsal, and ventral networks in episodic memory retrieval[END_REF] and targeted early in Alzheimer's disease, i.e. the quintessential disorder of memory impairment. The dorsolateral PFC is implicated in episodic memory encoding and retrieval [START_REF] Sandrini | The Role of Prefrontal Cortex in Verbal Episodic Memory: rTMS Evidence[END_REF], and the hippocampus, which is often found to be coupled with the DMN and nearby temporal regions, is critical for encoding and retrieval of recent information [START_REF] Squire | The neuropsychology of memory[END_REF]. The detected connections thus conform to what one would expect for the memory task employed. An interesting implication from our results is that even within a short duration of 30 minutes, detectable functional rewiring appears to be present, which provides further evidence for the shorter temporal scale of functional plasticity compared to its anatomical counterpart.

Fig. 1 .

 1 Fig. 1. Proposed matrix whitening transport versus parallel transport.

=

  be a d×d projected covariance matrix, i.e. a tangent vector, at ) (s C . To parallel transport T(c, s) from the tangent space at ) s ladder, we first determine using (3) the point on S ++ that is unit distance from ) (s C along the geodesic uniquely defined by T(c, s)[START_REF] Pennec | A Riemannian Framework for Tensor Computing[END_REF], which is exactly C(c, s) since )to R: G i (s), i = {1, …, g}, where g is the number of points along the geodesic, and find the midpoint of the geodesic joining G 1 (s) and C(c, s): M 1 (c, s) = Exp C(c, s) (0.5•Log C(c, s) (G 1 (s))). We then construct a geodesic from ) (s C to M 1 (c, s) and move twice the distance to find the one-step parallel transported covariance matrix:

Fig. 2 .

 2 Fig. 2. Classification accuracy and significantly discriminative connections.

Conclusions

We presented a Riemannian approach for classifying functional connectivity patterns before and after interventions. By using manifold operations to decouple connectivity features, significantly higher classification accuracy was obtained compared to directly using connectivity estimates as features, which are inherently inter-related. Also, higher detection sensitivity was shown with our proposed discriminative connection identification scheme compared to conventional permutation testing. The overall classification framework thus provides both accuracy and interpretability.