Parallel Transport on the Cone Manifold of SPD Matrices for Domain Adaptation
Or Yair Mirela Ben-Chen Ronen Talmon












Parallel Transport on the Cone Manifold of SPD Matrices for Domain Adaptation
    Or Yair, Student Member, IEEE, Mirela Ben-Chen,
and Ronen Talmon, Member, IEEEThe authors are with the Viterbi Faculty of Electrical Engineering, Technion-Israel Institute of Technology, Haifa 32000, Israel (e-mail: http://oryair@campus.technion.ac.iloryair@campus.technion.ac.il; http://ronen@ee.technion.ac.ilronen@ee.technion.ac.il).



     December 30, 2023
=======================================================================================================================================================================================================================================================================================================================================================


In this paper, we consider the problem of domain adaptation.
We propose to view the data through the lens of covariance matrices and present a method for domain adaptation using parallel transport on the cone manifold of symmetric positive-definite matrices. We provide rigorous analysis using Riemanninan geometry, illuminating the theoretical guarantees and benefits of the presented method. In addition, we demonstrate these benefits using experimental results on simulations and real-measured data.




positive definite matrices, domain adaptation, transfer learning, parallel transport.







§ INTRODUCTION


The increasing technological sophistication of current data acquisition systems gives rise to complex, multimodal datasets in high-dimension.
As a result, the acquired data do not live in a Euclidean space, and applying analysis and learning algorithms directly to the data often leads to subpar performance.



To facilitate the analysis and processing of such data, one approach is to observe complex high-dimensional data through the lens of objects with a known non-Euclidean geometry. Notable examples of such objects are SPD matrices, which live on a cone manifold with a Riemannian metric. One of the most common forms of SPD matrices is a covariance matrix, which captures the linear relations between the different data coordinates. These relations are typically simple to compute, and therefore, recently, have become popular features in many applications in computer vision, medical imaging, and machine learning <cit.>. In particular, in <cit.>, the Riemannian geometry of covariance matrices was studied and exploited for medical imaging and physiological signal analysis.

Typically, Riemannian geometry is used to map objects from the non-Euclidean manifold to a linear Euclidean space by projection onto a tangent plane of the manifold. In existing work, the use of Riemannian geometry is usually limited to a single tangent plane. This indicates a hidden assumption that the SPD matrices corresponding to the data are confined to a local region of the manifold. However, the SPD matrices of the data often do not live in a small neighborhood on the manifold, and thus the resulting calculations may be inaccurate. 

One such particular scenario, in which SPD matrices span a large portion of the cone manifold, occurs when the data comprise multiple domains corresponding to multiple sessions, subjects, batches, etc.

For example, we will show that in a Brain-Computer-Interface (BCI) experiment, the covariance matrices of data acquired from a single subject in a specific session capture well the overall geometric structure of the data. Conversely, when the data consist of measurements from several subjects or several sessions, then the covariance matrices do not live in the same region of the manifold.

Often, multi-domain data pose significant challenges to learning approaches. For example, in the BCI experiment, it is challenging to train a classifier based on data from one subject (session) and apply it to data from another subject (session). This problem is largely referred to as domain adaptation or transfer learning, and it has attracted a significant research effort in recent years <cit.>. 

Broadly, in domain adaptation, the main idea is to adapt a given model that is well performing on a particular domain, to a different yet related domain <cit.>. Specifically, in the context of the cone manifold of SPD matrices, previous work proposed (geometric) mean subtraction as a simple method for domain adaption of BCI data <cit.>. Although this approach provided reasonable results for overcoming the differences between multiple sessions of a single subject, we show here that it fails to overcome the differences between multiple subjects. In
<cit.>, a PT approach was proposed, which can be applied either directly to the data, or to a generative model of the data to reduce the computational load for large datasets. However, their approach considers a general Riemannian manifold. Since there is no closed-form expression of PT on Riemannian manifolds, besides the sphere manifold and the manifold of all SPD matrices <cit.>, no specific scheme or algorithm was provided. We note that PT can be approximated using Schild's Ladder <cit.>, an approach that has been used extensively on the manifold of imaging data <cit.>.  

In this paper, we propose a domain adaptation method using the analytic expression of PT on the cone manifold of SPD matrices. 
We claim that this is a natural and efficient solution for domain adaptation, which enjoys several important benefits. First, the solution is especially designed for SPD matrices, which have proven to be good features of data in a gamut of previous work <cit.>. Second, the analytic form of PT on the cone manifold circumvents approximations. Third, PT can be efficiently implemented, in contrast to the computationally demanding Schild's Ladder approximation. 
We establish the mathematical foundation of the proposed domain adaptation method. To this end, we provide new results in the geometry of SPD matrices. In addition, we show applications to both simulation and real recorded data, obtaining improved performance compared to the competing methods.

In parallel to our study, recent work <cit.> has proposed a scheme for transfer learning using the Riemannian geometry of SPD matrices, with a tight connection to the present work. We will show that the affine transformation proposed in <cit.> can be recast as PT. In this paper we provide the mathematical foundation to analyze this transport, we discuss the advantage of our solution compared to  <cit.>, and we point out the special case in which the two methods coincide.

This paper is organized as follows. In Section <ref>, we present preliminaries on the Riemannian geometry of SPD matrices. In Section <ref>, we formulate the problem, present the proposed domain adaptation method, and provide mathematical analysis and justification. Section <ref>
shows experimental results on both simulation and real data. Finally, we conclude the paper in Section <ref>.











§ PRELIMINARIES ON RIEMANNIAN GEOMETRY OF SPD MATRICES



In this section we provide the preliminaries regarding SPD matrices, and we refer the reader to the book <cit.> for a detailed exposition of this topic. We note that in this paper we focus on covariance matrices, however the statements also hold for general SPD matrices. By definition, an SPD  matrix P∈ℝ^n× n has only strictly positive eigenvalues. An alternative definition is that for any vector v≠0 the quadratic form is strictly positive, i.e., v^TPv>0.  



 §.§ Metric and distance

The definition of an SPD matrix entails that the collection of all SPD matrices constitutes a convex half-cone in the vector space of real n×n symmetric matrices.
This cone forms a differentiable Riemannian manifold ℳ equipped with the following inner product 

    ⟨S_1,S_2⟩ _𝒯_Pℳ=⟨P^-1/2S_1P^-1/2,P^-1/2S_2P^-1/2⟩,

where  𝒯_Pℳ is the tangent space at the point P∈ℳ, S_1,S_2∈𝒯_Pℳ, and ⟨·,·⟩ is the standard Euclidean inner product operation. The symmetric matrices S∈𝒯_Pℳ in the tangent plane live in a linear space, and therefore, we can view them as vectors (with a proper representation). Throughout this paper, we interchangeably use the terms vectors and symmetric matrices when referring to S∈𝒯_Pℳ. 

This Riemannian manifold is a Hadamard manifold, namely, it is simply connected and it is a complete Riemannian manifold of non-positive sectional curvature. Manifolds with non-positive curvature have a unique geodesic curve between any two points, a property that will later be exploited.
Specifically, the unique geodesic curve between any two SPD matrices P_1,P_2∈ℳ  is given by <cit.>

    φ(t)     = P_1^1/2(P_1^-1/2P_2P_1^-1/2)^tP_1^1/2,    0≤ t≤1.

 The arc-length of the geodesic curve defines the following Riemannian distance on the manifold <cit.>:

    d_R^2(P_1,P_2)     = ‖log(P_2^-1/2P_1P_2^-1/2)‖_F^2
        = ∑_i=1^nlog^2(λ_i(P_2^-1/2P_1P_2^-1/2))
        = ∑_i=1^nlog^2(λ_i(P_1P_2^-1)),

where P_1,P_2∈ℳ, ‖·‖ _F is the Frobenius norm, log(P) is the matrix logarithm, and λ_i(P) is the i-th eigenvalue of P. We additionally denote by ddtφ(t)=φ'(t)∈𝒯_φ(t)ℳ the velocity vector of the geodesic at t∈[0,1].
Figure <ref> presents an illustration of the geodesic curve and the Riemannian distance. The cone manifold of 2× 2 SPD matrices can be displayed in ℝ^3, since any symmetric matrix P=([ x y; y z ]) is positive if and only if x>0, z>0 and y^2<xz.





 §.§ Exponential and Logarithm maps

The Logarithm map, which projects an SPD matrix P_i∈ℳ to the tangent plane 𝒯_Pℳ at P∈ℳ, is given by 


    S_i  =Log_P(P_i)=P^1/2log(P^-1/2P_iP^-1/2)P^1/2∈𝒯_Pℳ.

The Exponential map, which projects a vector S_i∈𝒯_Pℳ back to the manifold ℳ is given by

    P_i =Exp_P(S_i)=P^1/2exp(P^-1/2S_iP^-1/2)P^1/2∈ℳ.

An important property relates the Logarithm and Exponential maps to the geodesic curve.
Formally, let P_1,P_2∈ℳ, and consider the (unique) geodesic φ(t) from P_1 to P_2. The initial velocity φ'(0)∈𝒯_P_1ℳ is given by the Logarithm map φ'(0)=Log_P_1(P_2).
Similarly, the Exponential map projects the initial velocity vector  φ'(0) back to P_2, namely, P_2=Exp_P_1(φ'(0)). 




 §.§ Riemannian mean


The Riemannian mean P of a set {P_i | P_i∈ℳ} is defined using the Fréchet mean:

    P≜min_P∈ℳ∑_id_R^2(P,P_i).

A special case is the Riemannian mean P of two SPD matrices P_1,P_2∈ℳ, which has a closed-form expression, and is located at the midpoint of the geodesic curve:

    P=φ(12)=P_1^1/2(P_1^-1/2P_2P_1^-1/2)^1/2P_1^1/2.

Generally, for more than two matrices, the solution of the optimization problem (<ref>) can be obtained by an iterative procedure. Barachant et al <cit.> presented an algorithm based on <cit.> for estimating the Riemannian mean. For completeness, we include their algorithm in Appendix <ref>.

Given a set {P_i|P_i∈ℳ} and its Riemannian mean P, there is a commonly used approximation of the Riemannian distances on ℳ in the neighborhood of P. 
Specifically, the approximation of the Riemannian distance d_R^2 is given by:

    d_R^2(P_i,P_j)≈‖Ŝ_i-Ŝ_j‖ _F^2  ,

where Ŝ_i=P^-12Log_P(P_i)P^-12.
For more details on the accuracy of this approximation, see <cit.>.




§ DOMAIN ADAPTATION WITH PARALLEL TRANSPORT




 §.§ Overview

Let 𝒳^(1)={x_i^(1)(t)} _i=1^N_1 and 𝒳^(2)={x_i^(2)(t)} _i=1^N_2 be two subsets of N_1 and N_2 high-dimensional time series, respectively, where x_i^(k)(t)∈ℝ^D. Suppose each subset lives in a particular domain, which could be related to the acquisition modality, session, deployment, and set of environmental conditions. In our notation, the superscript k denotes the index of the subset, the subscript i denotes the index of the time-series within each subset, and t represents the time axis of each time-series. 

Our exposition focuses only on two subsets, and the generalization for any number of subsets is discussed at the end of this section. In addition, we consider here time-series, but our derivation does not take the temporal order into account, and therefore, the extension to other types of data, where t is merely a sample index, e.g., images, is straight-forward.

Analyzing such data typically raises many challenges. For example, a long-standing problem is how to efficiently compare between high-dimensional point clouds, and particularly, time-series. When the data are measured signals, sample comparisons become even more challenging, since such high-dimensional measured data usually contain high levels of noise.

In particular, in our setting, we face an additional challenge, since the data is given in different domains; comparing time-series from the same subset is a difficult task by itself, even more so is comparing time-series from two subsets from different domains.

Our goal is to find a new joint representation of the two subsets in an unsupervised manner. Broadly, we aim to devise a low-dimensional representation in a Euclidean space that facilitates efficient and meaningful comparisons. 
For the purpose of evaluation, we associate the time-series x_i^(k)(t) with labels y_i^(k) and define “meaningful” comparisons with respect to these labels. 
More concretely, we evaluate the joint representation by the Euclidean distance between the new representation of any two time-series with similar corresponding labels, independently of the time-series respective domain.
We emphasize that the proposed approach is unsupervised and it does not depend on the labels, which are only considered for the purpose of evaluation.

Devising such a new representation will facilitate efficient and accurate domain adaptation schemes. 
Specifically, given a subset {x_i^(1)(t)} _i=1^N_1  with corresponding labels { y_i^(1)} _i=1^N_1, we could train a classifier based on the new derived representation of the subset. Then, when another unlabeled subset {x_i^(2)(t)} _i=1^N_2 becomes available, we could apply the trained classifier to the derived (joint) representation of the latter subset.



 §.§ Illustrative example


To put the problem setting and our proposed solution in context, throughout the paper, we will follow an illustrative example, taken from the brain computer interface (BCI) competition IV (dataset IIa) <cit.>.
Consider data from a BCI experiment of motor imagery comprising of recordings from D=22 Electroencephalography (EEG) electrodes. The dataset contains several subjects, where each subject was asked repeatedly to perform one out of four motor imagery tasks: movement of the right hand, the left hand, both feet, and the tongue.

Let 𝒳^(1)={x_i^(1)(t)} _i=1^N_1 be a subset of recordings acquired from a single subject, indexed (1), where the time-series x_i^(1)(t) consists of the signals, recorded simultaneously from the D EEG channels during the i-th repetition/trial. Each time series x_i^(1)(t) is attached with a label y_i^(1), denoting the imagery task performed at the the i-th trial.
Common practice is to train a classifier based on 𝒳^(1), so that the imagery task could be identified from new EEG recordings. This capability could then be the basis for devising brain computer interfaces, for example, to control prosthetics. 

Suppose a new subset 𝒳^(2)={x_i^(2)(t)} _i=1^N_2 of recordings acquired from another subject, indexed (2), becomes available. Applying the classifier, trained based on data from subject (1), to the new subset of recordings from subject (2) yields poor results, as we will demonstrate in Section <ref>. Indeed, most methods addressing this particular challenge, as well as related problems, exclusively analyze data from each individual subject separately. By constructing a joint representation for both 𝒳^(1) and 𝒳^(2), which is oblivious to the specific subject, we develop a classifier that is trained on data from one subject and applied to data from another subject without any calibration, i.e., without any labeled data from the new (test) subject.



 §.§ Covariance matrices as data features


As described before, we suggest looking at the data through the lens of covariance matrices. We denote the covariance matrices by: 

    P_i^(k)=𝔼[(x_i^(k)(t)-μ_i^(k))(x_i^(k)(t)-μ_i^(k))^T] ,

where μ_i^(k) = 𝔼[x_i^(k)(t)]. Typically, since the statistics of the data is unknown, we use estimates of the covariance, such as the sample covariance. We note that our approach is applicable to any kind of input data given as SPD matrices. For example, in machine learning, common practice is to use kernels which represent an inner product between features after some non-linear transformation <cit.>. 



By using covariance matrices as data features we enjoy a few key benefits. First, since covariance matrices are computed from data by averaging over time, they tend to be robust to noise. Second, covariance matrices can be seen as a low dimensional representation. Third, they have useful geometric properties and a well-developed Riemannian framework, as described in Section <ref>. Particularly, they have a Riemannian metric (<ref>), facilitating appropriate data samples comparisons, which is a basic ingredient of many analysis and learning techniques.
In this work, we build on and extend the latter.

Recently, the usefulness of covariance matrices has been demonstrated in the context of the BCI problem <cit.>.  
There, Barachant et al. considered data from a single subject and proposed to project the covariance matrices {P_i} of the recordings from each trial (after some whitening) into the tangent plane of the Riemannian mean P, namely compute S_i=Log_P(P_i). Then, a classifier was trained on the set {S_i}. Using this approach, state of the art results for motor imagery task classification were obtained.
However, when considering several subsets from multiple domains, such as different sessions or subjects, as reported in <cit.>, the covariance matrices convey a domain-specific content, which in turn poses limitations on task classification. 
For multiple sessions on different days, Barachant et al. proposed to subtract the Riemannian mean from each subset, namely, to project each subset 𝒫^(k)={P_i^(k)} to the tangent space at its own mean. Indeed, when the train set and the test set were obtained on different days, this mean normalization improved the task classification rate. However, in the case of multiple subjects, this approach is inadequate.  As mentioned before, given recordings from one subject as a train set and recordings from another subject as a test set, the classification of the different mental tasks based on covariance matrices fails completely. 


This illuminates the primary challenge addressed in this work – how to build a representation so that any two covariance matrices associated with the same mental task, but from possibly different sessions or subjects, will be given a similar representation.
Importantly, since the task labels are unknown, this objective cannot be directly imposed. In the sequel, we exploit the Riemannian geometry of covariance matrices, and devise such a representation in an unsupervised manner by preserving local geometric structures.



 §.§ Formulation


Consider two subsets 𝒫^(1) and 𝒫^(2) from two different domains consisting of N_1 and N_2 covariance matrices, respectively. Let  P^(1) and P^(2) be their respective Riemannian means. 
Let φ(t), given explicitly in (<ref>), denote the unique geodesic from  P^(2) to P^(1) such that φ(0)=P^(2) and φ(1)=P^(1).
Finally, let S_i^(k) be the symmetric matrix (or equivalently, the vector) in the tangent space 𝒯_P^(k)ℳ, obtained by projecting P_i^(k) to 𝒯_P^(k)ℳ:

    S_i^(k)=Log_P^(k)(P_i^(k)) ,

for k∈{1,2} and i∈{1,2,…,N_k} . 

Our goal now is to derive a new representation Γ(S_i^(2)) of S_i^(2) given by the map Γ:𝒯_P^(2)ℳ→𝒯_P^(1)ℳ, such that {S_i^(1)} and {Γ(S_i^(2)) } live in the same space. 
This allows us to relate samples from the two subsets, and compute quantities such as ⟨S_i^(1),Γ(S_j^(2))⟩ _P^(1). 
In addition, we require that the new representation will fulfill the following properties:



  * Zero mean:

    1N_2∑_i=1^N_2Γ(S_i^(2))=1N_1∑_i=1^N_1S_i^(1)=0



  * Inner product preservation:

    ⟨Γ(S_i^(2)),Γ(S_j^(2))⟩ _P^(1)=⟨S_i^(2),S_j^(2)⟩ _P^(2)

for all i,j ∈{1,…,N_2}.


  * Geodesic velocity preservation:

    Γ(φ'(0))=φ'(1)



Properties (1) and (2) imply that the new representation Γ preserves inter-sample relations, defined by the inner product. Note that a map Γ satisfying properties (1) and (2) is not unique;
for any Γ admitting to properties (1) and (2), the composition R∘Γ, where R is an arbitrary rotation within the subspace 𝒯_P^(1)ℳ, satisfies properties (1) and (2) as well.
To resolve this arbitrary degree of freedom, we use the geodesic between two points on the SPD manifold, which is unique <cit.>.
Concretely, in property (3), the two intrinsic symmetric matrices (vectors) φ'(0)∈𝒯_P^(2)ℳ and φ'(1)∈𝒯_P^(1)ℳ, induced by the velocity of the unique geodesic at the source and destination, are used to fix a rotation and to align the subset {Γ(S_i^(2))} with the subset {S_i^(1)}.  

We remark that the above properties imply that the subset {Γ(S_i^(2))} is embedded in the ⟨·,·⟩ _P^(1) inner product space. In the sequel, we will describe how to circumvent the dependence of the inner product space on P^(1) and make the new representation truly Euclidean by pre-whitening the data. Additionally, note that the mean subtraction presented in <cit.> admits only properties (1)-(2). 





 §.§ Domain adaptation

First, we explicitly provide the expression for parallel transport on the SPD cone manifold, and then we use it to define the map Γ.


Let A,B∈ℳ. The PT from B to A of any S∈𝒯_Bℳ is given by:

    Γ_B→A(S)≜ESE^T,

where E=(AB^-1)^12.

This lemma was presented in <cit.>. The proof of the lemma is given in Appendix <ref> and it is based on <cit.>.
An illustration of the PT on the SPD manifold is presented in Figure <ref>. 
Note that the inner products between the three vectors in the figure are preserved under the parallel transport Γ_B→A and the appearance could be misleading since the space is not Euclidean.




The representation Γ _P^(2)→P^(1)(S_i^(2)) given by the unique PT of S_i^(2) from P^(2) to P^(1) is well defined and satisfies properties (1)-(3).

The proof is given in Appendix <ref>.

Theorem <ref> sets the stage for domain adaptation. We propose a map Ψ:ℳ→ℳ that adapts the domain of the subset of SPD matrices 𝒫^(2) to the domain of the subset 𝒫^(1). For any P_i^(2)∈𝒫^(2), the map Ψ(P_i^(2)) is given by

    Ψ(P_i^(2))=Exp_P^(1)(Γ_P^(2)→P^(1)(Log_P^(2)(P_i^(2)))) .

To enhance the geometric insight, we explicitly describe the three steps comprising the construction of Ψ:
 


  * Project the SPD matrix P_i^(2) to the tangent plane 𝒯_P^(2)ℳ by S_i^(2)=Log_P^(2)(P_i^(2)).


  * Parallel transport S_i^(2) from  P^(2) to P^(1) by computing S_i^(2)→(1)=Γ_P^(2)→P^(1)(S_i^(2)).


  * Project the symmetric matrix  S_i^(2)→(1)∈𝒯_P^(1)ℳ back to the manifold using Exp_P^(1)(S_i^(2)→(1)).
 




The implementation of Ψ can be simplified and made more efficient by using the following theorem.



Let A,B,P∈ℳ and let S=Log_B(P)∈𝒯_Bℳ. Then,

    Exp_A(Γ_B→A(S))=EPE^T,

where E=(AB^-1)^12.


In words, the “parallel transport” of an SPD matrix P∈ℳ from B to A is given the same transformation applied to S=Log_B(P).
Namely, the “parallel transport” of the SPD matrix P from B to A is equal to projecting P to the tangent plane at B, parallel transporting the projection to the tangent plane at A, and then projecting back to the SPD manifold. 
As a consequence, we show in the sequel that the map Ψ in (<ref>) can be written simply in terms of Γ.

The proof of Theorem <ref> is given in Appendix <ref>. We note that we present the theorem in a general context, since we did not find such a result in the literature and believe it might be of independent interest.

Theorem <ref> enables us to efficiently compute Ψ(P_i^(2)), since it circumvents the computation of the Logarithm and Exponential maps of the SPD matrix in steps 1 and 3 above. Instead, the transformation defined by E is computed only once for the entire set, and (<ref>) can be recast as:

    Ψ(P_i^(2))=Γ_P^(2)→P^(1)(P_i^(2))=EP_i^(2)E^T

where E≜(P^(1)(P^(2))^-1)^12. Note that this equality is well defined since any tangent plane to the SPD manifold ℳ is the entire space of symmetric matrices <cit.>.

Thus far in the exposition, only the uniqueness of the geodesic curve on the manifold of SPD matrices was exploited, such that the PT along the geodesic admits the property in (<ref>), namely: Γ(φ'(0))=φ'(1). Importantly, PT specifically along the unique geodesic curve exhibits important invariance to the “relative” location on the manifold.



Two pairs (A_1,B_1) and (A_2, B_2), such that A_1,B_1,A_2,B_2 ∈ℳ, are equivalent if there exists an invertible matrix E such that

    A_2=Γ(A_1)=EA_1E^T 
    B_2=Γ(B_1)=EB_1E^T

We denote this relation by

    (A_1,B_1) ∼ (A_2, B_2)




The relation ∼ is an equivalence relation.

The proof is straight-forward as we show in the following.


  * Reflexivity is satisfied by setting E to be the identity matrix.

  * Symmetry: if  A_2=EA_1E^T then A_1=E^-1A_2E^-T and analogously for B_1, B_2.

  * Transitivity: if A_2=E_1A_1E_1^T and A_3=E_2A_2E_2^T then A_3=EA_1E^T where E=E_2E_1 and analogously for B_1, B_2.



In other words, two pairs are equivalent if the relation of the two matrices in the pair  is given by the same transformation Γ.
We interpret such equivalent pairs as matrices with equivalent intra-relations (e.g., if (A_1,B_1) ∼ (A_2, B_2), then d_R(A_1,B_1)=d_R(A_2,B_2)), but with a different global position on the manifold. For example, each two pairs in Figure <ref> are equivalent pairs.



Let (A_1,B_1) be a pair of SPD matrices A_1,B_1∈ℳ, and let [(A_1,B_1)] denote the equivalence class

    [(A_1,B_1)] = { (A_2, B_2) ∈ℳ×ℳ  |  (A_2,B_2) ∼(A_1,B_1) },

of all matrix pairs that are equivalent to (A_1,B_1).
Then, for any (A_2,B_2) ∈[(A_1,B_1)]:

    Γ∘Γ_B_1→A_1=Γ_B_2→A_2∘Γ ,

where Γ(P)=EPE^T and E is the transformation defined in Definition <ref>.


The proof is given in Appendix <ref>.

An immediate consequence of Proposition <ref> is that the domain adaptation via the representation Ψ is invariant to the relative position of P^(1) and P^(2) on the manifold, and is constructed equivalently for every pair in the equivalence class [(P^(1),P^(2))].

To demonstrate the importance of the property above, we revisit the illustrating BCI problem. Suppose (P_A1,P_B1) are the Riemannian means of the covariance matrices of Subject A and Subject B recorded in Session 1, and suppose (P_A2,P_B2) are the Riemannian means of the covariance matrices of Subject A and Subject B recorded in Session 2. If (P_A1,P_B1) ∼  (P_A2,P_B2), then there exists a transformation Γ such that Γ encodes the relation between Session 1 and Session 2 whereas the relation of the two subjects is encoded by Γ_P_B1→P_A1 or by Γ_P_B2→P_A2, depending on the session.

Proposition <ref> guarantees the consistence of the relation between Subject A and Subject B. Namely, the Riemannian mean of Subject B in Session 2 can be related to the Riemannian mean of Subject A in Session 1 using the relation between the sessions (given by Γ) and the relation between the two subjects (given either by Γ_P_B1→P_A1 or by Γ_P_B2→P_A2), independently of the relative location of the means on the manifold.






 §.§ Extension to K subsets


Overall, by Theorem <ref>, for a general number of subsets K≥ 2,  we can apply PT using Ψ  (<ref>) directly to the SPD matrices 𝒫^(k)={P_i^(k)} without projections to and from the tangent plane. 
Let  P̂  denote the Riemannian mean of Riemannian means (centroids) {P^(k)} _k=1^K of the subsets, namely,  
    P̂=min_P∑_k=1^Kd_R^2(P,P^(k)).

Each subset 𝒫^(k) is then parallel transported from its corresponding centroid P^(k) to P̂. Formally, let Γ_i^(k) denote P_i^(k) after applying PT, which is given by

    Γ_i^(k)=Γ_P^(k)→P̂(P_i^(k)),    ∀ i, k,

and let S̃_i^(k) be the projection of  Γ_i^(k) to the Euclidean tangent space  (<ref>):

    S̃_i^(k)=log(P̂^-12Γ_i^(k)P̂^-12).

This projection, which is further discussed in <cit.>, can be interpreted as (i) data whitening by  Γ̃_i^(k)=P̂^-12Γ_i^(k)P̂^-12, and (ii) projection to 𝒯_Iℳ where I is the identity matrix. The projected symmetric matrices (vectors) S̃_i^(k) indeed reside in a Euclidean space. 
The proposed algorithm is given in Algorithm <ref>.



We conclude this section with two remarks. First, since the matrices S̃_i^(k) are symmetric, only their upper (or lower) triangular part with a gain factor of √(2)  applied to all non-diagonal elements could be taken into account. 

Second, alternative choices of P̂ could also be used, for example, the identity matrix.
Indeed, recently <cit.> proposed to align datasets for transfer learning in a similar context using the identity matrix as P̂.
However in <cit.>, the alignment appeared as an empirical affine transformation, whereas in this work, we provide the geometrical justification and rigorous mathematical analysis.

Specifically, in the case of two subsets with means P_A∈ℳ and P_B∈ℳ, the affine transformation presented in <cit.> can be interpreted as two consecutive applications of PT: from P_B to I and then from I to P
_A. 
The arbitrary choice of I as an intermediate point introduces dependence of the algorithm on the global position on the manifold. Indeed, such a procedure, which can be expressed by Γ _I→P_A∘Γ _P_B→I does not admit the invariance property specified in Proposition <ref>.

Interestingly, the method proposed in <cit.> coincides with the present work, namely, Γ _P_B→P_A =Γ _I→P_A∘Γ _P_B→I when the identity matrix  I is on the geodesic φ between P_B and  P_A. In this case, the matrices P_A and P_B commute and they have the same eigenvectors (see Appendix <ref>).
From a data analysis perspective, when P_A and P_B are two covariance matrices of two subsets, this implies that the subsets have the same principal components.

We set P̂ as the Riemannian mean of the centroids so that the overall transport applied to the covariance matrices is minimal. This choice is motivated by the assumption that transporting accumulates distortions. 
This is a straight-forward generalization of the two subsets case, where the parallel transport is carried out along the shortest path (unique geodesic curve).





§ EXPERIMENTAL RESULTS


In this section we show the results of  Algorithm <ref> for both a synthetic example and for real data.



 §.§ Toy Problem

We generate time series in ℝ^2, so that their covariance matrices are in ℝ^2× 2. Since the covariance matrices are symmetric, this particular choice enables us to visualize them in ℝ^3. Concretely, any 2× 2 symmetric matrix A=([ x y; y z ]) can be visualized in ℝ^3 using (x,y,z)∈ℝ^3. A  is positive-definite if and only if: x,z>0 and y^2<x z. These conditions establish the cone manifold of 2× 2 SPD matrices.

Consider the set of hidden multi-dimensional times series {s_i[n]} _i=1^100, given by:

    s_i[n]=[     sin(2π f_0n/T); cos(2π f_0n/T+ϕ_i) ],     n=0,…,T-1

where f_0=10, T=500, and  ϕ_i is uniformly drawn from [-π2,0].
Namely, each time-series s_i[n] consists of two oscillatory signals and is governed by a 1-dimensional hidden variable ϕ_i, the initial phase of the oscillations.
Indeed, the population covariance of s_i[n] is

    1/2[         1 -sin(ϕ_i); -sin(ϕ_i)         1 ]

which depends only on ϕ_i, and therefore, when presenting the population covariances of the time-series s_i[n] in ℝ^3, two coordinates are fixed and only one varies with i.

We generate two observable subsets, 𝒳^(k)={x_i^(k)[n]} _i=1^100, k=1,2 such that: 
    x_i^(k)[n]=M^(k)s_i[n],

where M^(1) is randomly chosen, and M^(2)=1.5([ -1  0;  0  1 ])M^(1).

The two subsets  𝒳^(1) and 𝒳^(2) can be viewed as two different observations of {s_i[n]} through two unknown observation functions M^(1) and M^(2). For example, 𝒳^(1) and 𝒳^(2)  can represent two different batches, and M^(1) and M^(2) can represent the discrepancy between two different sessions of a particular experiment.
For each x_i^(k)[n], we compute its sample covariance matrix by

    P_i^(k)=1T∑_nx_i^(k)[n](x_i^(k)[n])^T=M^(k)P_s_i(M^(k))^T,
where  P_s_i denotes the inaccessible sample covariance of s_i[n], which is given by:

    P_s_i=1T∑_n=0^T-1s_i[n](s_i[n])^T.

Our goal is to obtain a new representation of the observed data both in  𝒳^(1) and 𝒳^(2), which circumvents the effect of  M^(1) and M^(2). 
Moreover, in the new representation, we aspire to associate two observations from possibly different subsets which have a similar initial phase ϕ_i. 



In Figure <ref> we plot the 2× 2 SPD matrices in ℝ^3, where the black points mark the boundaries of the cone manifold. The red line marks the center of the cone, given by αI  for α∈[0,2], and the blue point on the red line indicates the identity matrix I, namely, where α = 1. 

Figure <ref>(a) presents the two subsets 𝒫^(k)={P_i^(k)}, k=1,2 of accessible sample covariance matrices, colored by ϕ_i (left) and by k (right).
We observe that the two subsets 𝒫^(1) and 𝒫^(2) are completely separated, while each subset has a similar structure governed by the values of ϕ_i.

We apply Steps (1)-(3) of Algorithm <ref> to the subsets 𝒫^(k) to obtain {Γ_i^(k)}. 

Figure <ref>(b) presents {Γ_i^(k)}, colored by ϕ_i (left) and k (right). Now we observe that in the new representation, the two subsets are aligned, namely the discrepancy caused by M^(1) and M^(2) is removed while the intrinsic structure given by ϕ_i is preserved. As a result, we can associate covariance matrices from different batches but with similar underlying ϕ_i values. Note that this was accomplished by Algorithm <ref> in a completely unsupervised manner, without access to the hidden “labels” ϕ_i.









 §.§ BCI - Motor Imagery 


As described in Section <ref>, we use data from the BCI competition IV <cit.>.
The dataset contains EEG recordings acquired by 22 EEG electrodes from 9 subjects, where the data from each subject was recorded on 2 different days of experiments. 
The experiment protocol consists of repeated trials, where in each trial the subject was asked to imagine performing one out of four possible movements: (i) right hand, (ii) left hand, (iii) both feet, and (iv) tongue. Overall, in a single day, each movement was repeated 72 times by each subject, and therefore, the dataset contains 288 trials from each subject in each day of experiments.

We remark that all the algorithms participating in the competition reported on poor classification results for particular four subjects. 
Since our goal is not to improve the classification of the data from each subject,  we excluded these four subjects (indexed 2,4,5,6).


Initially, focusing on the data from a single subject, we show that Algorithm <ref> builds a representation of the data which enables us to train a classifier with data from one day of experiments and apply it to data from the other day of experiments. 
Then, we further show that Algorithm <ref> builds a representation that allows us to  train a classifier based on data from one subject and apply it to data from a different subject without any additional labeled trials. 
Finally, we extend the latter result and show the performance on multiple subjects.



  §.§.§ Single Subject – Different Days

In the first experiment, we process the recordings of Subject 8 (arbitrarily chosen) from the two days of experiments. We report that the results for the other 4 subjects were similar.
We denote the subsets of trial recordings from day k=1,2 by 𝒳^(k)={x_i^(k)} _i=1^288. From the recordings of each trial i, we compute the sample covariance matrix P_i^(k)∈ℝ^22×22, and denote 𝒫^(k)={P_i^(k)}. 

To highlight the challenge, we first compute P̂, the Riemannian mean of all covariance matrices (from both subsets). Then, we project the matrices onto 𝒯_P̂ℳ by computing B_i^(k)=Log_P̂(P_i^(k)).  For visualization purpose, we apply tSNE <cit.> to the vectors {B_i^(k)}. 
Figure <ref> (a) presents the two dimensional representation of the vectors obtained by the tSNE algorithm. Namely, each point in the figure is the representation of a vector B_i^(k). On the left, the points are colored according to the different days (indexed by k=1,2), and on the right, the points are colored according to the mental task. We observe that, similarly to the toy problem, the recordings from the different days are completely separated. 
This implies that one cannot train a classifier from the recordings from day 1 and apply it to the recordings from day 2. 

We apply Algorithm <ref> to the subsets 𝒫^(k) covariance matrices and obtain the subsets 𝒮̃^(k)={S̃_i^(k)}. Figure<ref> (b) presents the two dimensional representation of the vectors 𝒮̃^(k) obtained by the tSNE algorithm. On the left, the points are colored according to the different days, and on the right, the points are colored according to the mental task. We observe that the difference between the different days of experiments is completely removed.
More importantly, we further observe that the new representations of the two subsets are aligned, i.e., we obtain similar representations of two recordings associated with the same mental task, regardless of their respective sessions (days of experiments).














  §.§.§ Two Subjects


We repeat the evaluation, but now with the two subsets  𝒳^(k), k=1,2 which were recorded from two subjects, specifically, Subject 3 and Subject 8. We repeat the steps from the previous examination. We compute P̂, the Riemannian mean of all covariance matrices (from both subsets). Then, we project the covariance matrices onto 𝒯_P̂ℳ and apply tSNE to obtain two dimensional representations. Figure <ref> (a) presents the two dimensional representations obtained by the tSNE algorithm. On the left, the points are colored by the subject index, and on the right the points are colored according to the mental task. Similarly to the single-subject two-sessions case, we observe that the recordings from different subjects are completely separated in the obtained representation. 

We also apply the mean transport approach presented in <cit.>. The mean transport is obtained by projecting each subset of covariance matrices 𝒫^(k) to its own tangent plane 𝒯_P^(k)ℳ, where P^(k) is the Riemannian mean of the k-th subset. In other words, we compute S_i^(k)=Log_P^(k)(P_i^(k)). Figure <ref> (b) presents the two dimensional representation obtained by the tSNE algorithm.  We observe that indeed the two subsets are not separated as in Figure <ref>(a), however, the inner structure of each subset was not preserved. Thus, this scheme is insufficient and does not support training a classifier based on data from one subject and applying it to data from another subject. 

Finally, we apply Algorithm <ref> to the subsets 𝒫^(k), and obtain the subsets 𝒮̃^(k)={S̃_i^(k)}. Figure <ref>(c) presents the two dimensional representation of the subsets 𝒮̃^(k) obtained by the tSNE algorithm.  We observe that in this representation, the subsets are not separated. Moreover, the two subsets are aligned according to the mental tasks, and indeed points that correspond to the same mental task assumed a similar value in the new representation. This new representation allows us to train a classifier using recordings from one subject and apply it to recordings from another subject.







  §.§.§ Multiple Subjects


In the third experiment, we apply Algorithm <ref> to multiple subjects. We processed data from five subjects (indexed 1,3,7,8 and 9) and from all trials corresponding to only three of the mental tasks: left hand, foot, and tongue. Namely, we omit one class for visualization purposes. We denote the subsets of the recordings by 𝒳^(k)={x_i^(k)} for k=1,3,7,8,9. As before, we compute the covariance matrices 𝒫^(k)={P_i^(k)}.  We compute P̂, the Riemannian mean of all covariance matrices. Then, we project the covariance matrices onto 𝒯_P̂ℳ and apply tSNE to obtain a two dimensional representation. Figure <ref> (a) presents the two dimensional representation obtained by the tSNE algorithm. On the left, the points are colored by the subject index, and on the right the points are colored according to the mental task. As before, we observe that the recordings from different subjects are completely separated. 

Next, we apply Algorithm <ref> to the five subsets of covariance matrices 𝒫^(k) and obtain five subsets of new representations 𝒮̃^(k)={S̃_i^(k)}. Figure <ref> (b) presents the two dimensional representation of the subsets 𝒮̃^(k) obtained by the tSNE algorithm. We observe that also in the multiple subjects scenario, Algorithm <ref> was able to center and align the subsets. 



To provide quantitative results, we apply a leave-one-subject-out cross validation, namely, we trained a classifier based on 4 out of the 5 subjects and evaluated the classification accuracy for each one of the three methods mentioned in Section <ref>. We compare the classification accuracy of  Algorithm <ref>  to the two other approaches denoted by: (i) “Baseline (No Transport)”, and (ii) the “Mean Transport” approach proposed in <cit.>.






Figure <ref> presents the classification accuracy obtained using the three competing methods as a function of the evaluated subject. In all cases, applying Algorithm <ref> to the  data dramatically improved the classification accuracy.
In addition, Figure <ref> presents the confusion matrices when Subject 3 was tested. Figure <ref> (left) presents the confusion matrix obtained from the data without applying any transportation (“Baseline”). Figure <ref> (center) presents the confusion matrix obtained by the “Mean Transport” approach. Figure <ref> (right) presents the confusion matrix obtained by Algorithm <ref>. 
We observe that Algorithm <ref> obtains significantly better classification results compared with the “Baseline” and “Mean Transport” algorithms. In addition, the confusion matrices highlight the challenge in training a classifier from multiple subject data. For example, in both confusion matrices in Figure <ref> (left) and (center), the “foot” mental task falsely dominated the prediction (it was predicted 183 and 164 times, respectively, whereas is was performed only 72 times). This implies that the decision regions of the classifiers are completely misaligned with the data from a new unseen subject.












 §.§ Sleep Stage Identification

Here, we demonstrate the applicability of Algorithm <ref> to real medical signals. Specifically, we address the problem of sleep stage identification. Typically, for this purpose, data are collected in sleep clinics with multiple multimodal sensors, and then, analyzed by a human expert. There are six different sleep stages: awake, REM, and sleep stages 1-4, indicating shallow to deep sleep.

The data we used is available online in <cit.> and described in detail in <cit.>. A single patient's night recording contains several measurements including two EEG channels and one electrooculography (EOG) channel sampled at 100[Hz].
We used recordings from three subjects. We split each subject's night into non-overlapping 30 seconds windows. We omit the awake and sleep stage 4 windows due to too few occurrences. For visual purposes, we kept only windows corresponding to REM and stage 3.

We denote the i-th window of the k-th subject by x_i^(k)(t) with its corresponding covariance matrix P_i^(k)∈ℝ^3× 3. 
We first compute P̂, the Riemannian mean of all covariance matrices (from the three subsets). Then, we project the matrices onto 𝒯_P̂ℳ by computing S_i^(k)=Log_P̂(P_i^(k)). For visualization purposes, we apply PCA to the vectors {S_i^(k)} and present the first two principle components.
Since the covariance matrices in this experiment are of size 3 × 3, dimension reduction using PCA was sufficient. It was preferred here over tSNE since it preserves the global geometry of the data.





Figure <ref> (a) presents the two dimensional representation of the vectors obtained by PCA. Namely, each point in the figure is the representation of a vector S_i^(k). On the left, the points are colored according to the different subjects, and on the right, the points are colored according to the sleep stage. We observe that the points are clustered according to the different subjects.
We apply Algorithm <ref> to the three subsets {P_i^(k)} of covariance matrices and obtain the subsets {S̃_i^(k)}. Figure <ref> (b) presents the two dimensional representation of the vectors {S̃_i^k)} obtained by PCA. On the left, the points are colored according to the different subjects, and on the right, the points are colored according to the sleep stage. Now we observe that the data is clustered according to the sleep stage while the difference between the three subjects is completely removed.

As in the Subsection <ref>, to provide quantitative results, we train a classifier based on Subject 1 and Subject 2 and evaluate the classification accuracy on Subject 3. 
Figure <ref> presents the obtained confusion matrices. Figure <ref> (left) presents the confusion matrix obtained from the data without applying any adaptation (“Baseline”). Figure <ref> (center) presents the confusion matrix obtained by the “Mean Transport” approach. Figure <ref> (right) presents the confusion matrix obtained by Algorithm <ref>. 
We observe that using Algorithm <ref> demonstrates better classification results compared with the “Baseline” and the “Mean Transport” algorithms. 














§ CONCLUSIONS


Analyzing complex data in high-dimension is challenging, since such data do not live in a Euclidean space. Therefore, basic operations such as comparisons, additions, and subtractions, which are the basis of any analysis and learning technique, do not necessarily exist and are not appropriately  defined. 

In this work, we propose to view the complex data through the lens of SPD matrices, which reside on an analytic Riemannian manifold. 
Using the Riemannian geometry of SPD matrices, we presented an approach for multi-domain data representation. Based on this new representation, we proposed an algorithm for domain adaptation. 

We extend the existing results in the Riemannian geomery of SPD matrices and establish a framework for the justification and analysis of the proposed solution.

We demonstrated the usefulness of the presented domain adaptation method in applications to simulation and real recorded data.





§ PROOFS OF LEMMA <REF>


We prove lemma <ref>  presented in <cit.>.

The PT of S  along the geodesic between B and A is given by <cit.>:

    Γ_B→A(S)=MSM^T

where M=B^1/2exp(B^-1/21/2Log_B(A)B^-1/2)B^-1/2. Now we will show that M can be written more simply, proving a more efficient way to compute it.
We have

    M    =B^1/2exp(B^-1/21/2Log_B(A)B^-1/2)B^-1/2
        =B^1/2exp(1/2log(B^-1/2AB^-1/2))B^-1/2
        =B^1/2(B^-1/2AB^-1/2)^1/2B^-1/2

and also

    M^2=(B^1/2(B^-1/2AB^-1/2)^1/2B^-1/2)^2=AB^-1=E^2

and since AB^-1 is similar to B^-1/2AB^-1/2>0, it has only positive eigenvalues and the square root is unique, namely E=M.




§ PROOF OF THEOREM <REF>



For better readability we denote A=P^(1) and B=P^(2).
First we remark that Γ_B→A is well defined since 𝒯_Pℳ is the space of all symmetric matrices regardless the matrix P, and  if the input S is symmetric then by definition Γ_B→A(S) is symmetric as well. 
Any tangent plane to the SPD manifold ℳ is the entire space of symmetric matrices <cit.>.

Condition (1) is immediate since Γ is a linear operation, and therefore we have

    ∑_i=1^N_2Γ(S_i^(2))=∑_i=1^N_2ES_i^(2)E^T=E∑_i=1^N_2S_i^(2)_=0E^T=0

Conditions (2) and (3) are derived from Lemma <ref>, since these are properties of PT.
For completeness, we provide their explicit proofs.
Proof of condition (2):
Let A,B∈ℳ and S_1,S_2∈𝒯_Bℳ and denote E=(AB^-1)^1/2=B^1/2(B^-1/2AB^-1/2)^1/2B^-1/2.
We have

    A^-1E    =A^-1B^1/2(B^-1/2AB^-1/2)^1/2B^-1/2
        =A^-1B^1/2B^-1/2AB^-1/2(B^-1/2AB^-1/2)^-1/2B^-1/2
        =B^-1/2(B^-1/2AB^-1/2)^-1/2B^-1/2

Namely,  A^-1E is a symmetric matrix.
Thus, we get that

    E^TA^-1E    =E^TE^TA^-1
        =(AB^-1)^TA^-1
        =(B^-1A)A^-1
        =B^-1

and finally, we obtain

    ⟨ES_1E^T,ES_2E^T⟩ _A    =⟨ES_1E^TA^-1,A^-1ES_2E^T⟩
        =Tr{ES_1E^TA^-1ES_2E^TA^-1}
        =Tr{S_1E^TA^-1ES_2E^TA^-1E}
        =Tr{S_1B^-1S_2B^-1}
        =⟨S_1,S_2⟩ _B

Proof of  condition (3):
Let B∈ℳ be an SPD matrix with the following spectral decomposition:

    B=MΛM^T

Then, we have

    d/dtB^t    =d/dtMΛ^tM^T
        =MΛ^tlog(Λ)M^T
        =MΛ^tM^TMlog(Λ)M^T
        =B^tlog(B)

Consider the geodesic φ(t) from B to A:

    φ(t)=B^1/2(B^-1/2AB^-1/2)^tB^1/2

Thus, its velocity at t=0 is given by

    φ'(0)=B^1/2log(B^-1/2AB^1/2)B^1/2=Log_B(A)

and similarly, the velocity at t=1 is given by

    φ'(1)     =B^1/2(B^-1/2AB^-1/2)^1log(B^-1/2AB^-1/2)B^1/2
        =AB^-1/2log(B^-1/2AB^-1/2)B^1/2
        =-AB^-1/2log(B^1/2A^-1B^1/2)B^1/2
        =-AB^-1/2log(B^1/2A^-1/2A^-1/2BA^-1/2A^1/2B^-1/2)B^1/2
       =_(*)-AB^-1/2B^1/2A^-1/2log(A^-1/2BA^-1/2)A^1/2B^-1/2B^1/2
        =-A^1/2log(A^-1/2BA^-1/2)A^1/2
        =-Log_A(B)

where in (*) we pull out V=B^1/2A^-1/2 and V^-1=A^1/2B^-1/2 from the log , since it is a scalar function: log(VPV^-1)=Vlog(P)V^-1. 

Let U be the following unitary matrix

    U=A^-1/2B^1/2(B^-1/2AB^-1/2)^1/2

Using U, we can rewrite E as:

    E=(AB^-1)^1/2=B^1/2(B^-1/2AB^-1/2)^1/2B^-1/2=A^1/2UB^-1/2

Finally, by combining (<ref>), (<ref>) and <ref>, we have

    Eφ'(0)E^T    =A^1/2UB^-1/2B^1/2log(B^-1/2AB^-1/2)B^1/2B^-1/2U^TA^1/2
        =A^1/2log(UB^-1/2AB^-1/2U^T)A^1/2
        =A^1/2log(A^1/2B^-1A^1/2)A^1/2
        =-A^1/2log(A^-1/2BA^-1/2)A^1/2
        =-Log_A(B)
        =φ'(1)





§ PROOF OF THEOREM <REF>




    Ψ(P)     =Exp_A(Γ_B→A(S))
        =Exp_A(ELog_B(P)E^T)
        =Exp_A(EB^1/2log(B^-1/2PB^-1/2)B^1/2E^T)
        =A^1/2exp(A^-1/2EB^1/2log(B^-1/2PB^-1/2)B^1/2E^TA^-1/2)A^1/2
        =A^1/2exp(Ulog(B^-1/2PB^-1/2)U^T)A^1/2
        =A^1/2Uexp(log(B^-1/2PB^-1/2))U^TA^1/2
        =A^1/2UB^-1/2PB^-1/2U^TA^1/2
        =A^1/2A^-1/2EB^1/2B^-1/2PB^-1/2B^1/2E^TA^-1/2A^1/2
        =EPE^T


where U=A^-1/2EB^1/2 is a unitary matrix, and therefore can be pulled out of the scalar exp function..




§ PROOF OF PROPOSITION <REF>





Let E_1=(A_1B_1^-1)^1/2=B_1^1/2(B_1^-1/2A_1B_1^-1/2)B_1^-1/2 and E_2=(A_2B_2^-1)^1/2=B_2^1/2(B_2^-1/2A_2B_2^-1/2)B_2^-1/2.

Since

    Γ(Γ_B_1→A_1(S))=EE_1SE_1^TE^T ,

 and 

    Γ_B_2→A_2(Γ(S))=E_2ESE^TE_2^T ,

it is enough to show that

    E_2E=EE_1.


First, note that K=B_2^1/2E^-TB_1^-1/2 is unitary:

    KK^T    =B_2^1/2E^-TB_1^-1/2(B_2^1/2E^-TB_1^-1/2)^T
        =B_2^1/2E^-TB_1^-1/2B_1^-1/2E^-1B_2^1/2
        =B_2^1/2E^-TB_1^-1E^-1B_2^1/2
        =B_2^1/2B_2^-1B_2^1/2
        =I .

Now, we have (<ref>).
Finally, we have

    E_2E    =B_2^1/2(B_2^-1/2A_2B_2^-1/2)^1/2B_2^-1/2E
        =B_2^1/2(B_2^-1/2A_2B_2^-1/2)^1/2B_2^1/2B_2^-1E
        =EB_1^1/2(B_1^-1/2A_1B_1^-1/2)^1/2B_1^1/2E^TB_2^-1E
        =EB_1^1/2(B_1^-1/2A_1B_1^-1/2)^1/2B_1^-1/2
        =EE_1 .





§  


In this appendix we proof that if the identity matrix I is on the geodesic φ between A and  B, then, the matrices A and B commute and they have the same eigenvectors.

if I is on the geodesic φ(t), then there exist some t_0∈(0,1) (if t_0=0 or t_0=1 the result is trivial) such that:

    φ(t_0)=A^1/2(A^-1/2BA^-1/2)^t_0A^1/2=I

Now, consider A=VΛV^T the eigenvalue-decomposition of A. By multiplying both from the right and the left, we have

    A^1/2(A^-1/2BA^-1/2)^t_0A^1/2    =I
    (A^-1/2BA^-1/2)^t_0    =A^-1

By raising to the power of 1/t_0, and then multiplying both from the right and the left again, we have

    A^-1/2BA^-1/2    =A^-1/t_0
    B    =A^1-1/t_0
    B    =BΛ^1-1/t_0V^T



Thus, B has the same eigenvectors as A and they commute

    AB    =VΛV^TVΛ^1-1/t_0V^T
        =VΛΛ^1-1/t_0V^T
        =VΛ^1-1/t_0V^TVΛV^T
        =BA





§ RIEMANNIAN MEAN ALGORITHM











IEEEtran














