# Approximating Sets of Symmetric and Positive-Definite Matrices by Geodesics

Authors: L. Machado; F. Silva Leite
Journal: Conference Papers in Mathematics
 (2013)

---

## Abstract

We formulate a generalized version of the classical linear regression problem on Riemannian manifolds and derive the counterpart to the normal equations for the manifold of symmetric and positive definite matrices, equipped with the only metric that is invariant under the natural action of the general linear group.


---

## Body

## 1. Introduction

The geometry of the set of symmetric and positive definite (SPD) matrices is in the focus of intensive research activity involving tensor analysis. The importance of SPD matrices lies in the fact that they encode image information. As a consequence, they appear in several contexts of computer vision, such as, for instance, in [1, 2], in medical image analysis to interpolate scattered diffusion tensor magnetic resonance imaging, [3, 4],  but also in the area of continuum physics related to averaging methods for the case of the elasticity tensor of the generalized Hooke’s law [5]. Since the set of SPD matrices has a natural structure of Riemannian manifold, the rich theory of differential geometry can be used to solve real problems that may be formulated on this manifold.One particular problem of interest, that as far as we know has not been studied before, is that of approximating a set of data points in the SPD manifold by a geodesic. In the present paper, we first formulate the problem of finding the geodesic that best fits a given set of time-labelled points on a general Riemannian manifold. This corresponds to the natural generalization of the classical linear regression problem inℝn. Solving this problem on a Riemannian manifold requires knowing explicit formulas for geodesics and for the geodesic distance between two points. Such is the case of connected and compact Lie groups, where geodesics are one-parameter subgroups or their translations, or of Euclidean spheres, where geodesics are the great circles. The geodesic regression problem has been studied for these two cases in [6], and the numerical implementation of the spherical case will appear soon in [7].Our main objective is to derive the counterpart of the normal equations when the given data lies in the SPD manifold, equipped with a particular Riemannian metric that is affine-invariant. The paper is organized as follows. In Section2, we start with the formulation of the geodesic fitting problem for the general case of a geodesically complete Riemannian manifold. In Section 3, we specialize to the case of the SPD manifold endowed with its natural affine-invariant Riemannian metric and gather all the necessary background to achieve our goal. The main result appears in Section 3.1, where, in particular, we derive the normal equations for the geodesic fitting problem. We conclude the paper with some remarks in Section 4.


## 2. Geometric Notions

We start with basic notions in differential and Riemannian geometry that are important throughout the paper. For more details, we refer to classical books on the subject, as for instance, [8–10].LetM be an n-dimensional Riemannian manifold with Riemannian metric 〈·,·〉 and ∇ its Levi-Civita connection. Given p∈M, TpM denotes the tangent space of M at p, and TM stands for the tangent bundle of M.A vector fieldV along a curve c:I⊂ℝ→M is a mapping that assigns to each t∈I, the vector V(t)∈Tc(t)M. If c is smooth, the velocity vector field, dc/dt, is an example of a smooth vector field in M. If V is induced by the vector field X:M→TM, that is, if V(t)=Xc(t), then the covariant derivative of V along c is defined by
(1)DVdt=∇dc/dtX.
A vector field V along a curve c is said to be parallel if
(2)DVdt=0.
By definition, a geodesic is a smooth curve c whose velocity vector field is parallel along c. That is,
(3)Ddt(dcdt)=0.
This is a second order differential equation and according to the theory of existence and uniqueness for ordinary differential equations, given p∈M and v∈TpM, there exists a unique geodesic c:[0,1]→M, satisfying c(0)=p and (dc/dt)(0)=v. This geodesic can be parameterized explicitly by
(4)c(t)=expp(tv),
where expp(v) is the point on the geodesic c that is at a distance equal to ∥v∥ away from p.In spite of the fact that the exponential map,exp:TM→M, is only a terminology, there are some special Riemannian manifolds where it can be explicitly defined. For instance, in Euclidean spaces the geodesic starting at p with initial velocity v are the straight lines:
(5)expp(tv)=p+tv.
For the unit sphere Sn⊂ℝn+1, equipped with the Riemannian metric induced from the Euclidean metric of the embedding space, geodesics are great circles. In particular, the geodesic starting at p with velocity v is given by
(6)expp(tv)=pcos(t∥v∥)+v∥v∥sin(t∥v∥).
The manifold M is said to be geodesically complete if every pair of points in M can be joined by a geodesic. In this case, (see, for instance [11]) the geodesic that joins p at t=0 to q at t=1 is parameterized explicitly by
(7)c(t)=expp(texpp-1(q)),
and the distance between p and q, which is simply the length of this minimizing geodesic, is given by
(8)d(p,q)=〈expp-1(q),expp-1(q)〉1/2.


### 2.1. Statement of the Problem

We are now in conditions to formulate a natural generalization of the linear regression problem on a Riemannian manifold.Given a finite set of points inM, p0,p1,…,pN, and a set of instants of time t0<t1<⋯<tN, for simplicity assumed to form a partition of the unit time interval [0,1], find a parameterized geodesic
(9)γ:[0,1]⟶Mt↦γ(t)=expp(tv),
that yields the minimum value for the functional
(10)E(γ)=∑i=0N‍d2(pi,γ(ti)),
where d is the distance function defined by (8).Solving this optimization problem is only possible for manifolds where explicit forms for geodesics are available. Theoretical results for the cases whenM is the unit n-sphere or a connected and compact Lie group appeared in [6] and, more recently, analogous results were obtained for the Grassmann manifold in [12] and for the SPD manifold equipped with the Log-Euclidean metric in [13]. A numerical approach to solve the problem on the unit n-sphere has been undertaken in [7]. In the next section we study the geodesic fitting problem for the SPD manifold equipped with its natural affine-invariant Riemannian metric.


## 2.1. Statement of the Problem

We are now in conditions to formulate a natural generalization of the linear regression problem on a Riemannian manifold.Given a finite set of points inM, p0,p1,…,pN, and a set of instants of time t0<t1<⋯<tN, for simplicity assumed to form a partition of the unit time interval [0,1], find a parameterized geodesic
(9)γ:[0,1]⟶Mt↦γ(t)=expp(tv),
that yields the minimum value for the functional
(10)E(γ)=∑i=0N‍d2(pi,γ(ti)),
where d is the distance function defined by (8).Solving this optimization problem is only possible for manifolds where explicit forms for geodesics are available. Theoretical results for the cases whenM is the unit n-sphere or a connected and compact Lie group appeared in [6] and, more recently, analogous results were obtained for the Grassmann manifold in [12] and for the SPD manifold equipped with the Log-Euclidean metric in [13]. A numerical approach to solve the problem on the unit n-sphere has been undertaken in [7]. In the next section we study the geodesic fitting problem for the SPD manifold equipped with its natural affine-invariant Riemannian metric.


## 3. The SPD Manifold

For the purpose of this paper, the SPD manifold consisting of all symmetric and positive definite matrices will be equipped with a particular Riemannian metric which is affine invariant. This metric has been used in [14, 15] where other important results concerning the SPD manifold also appeared. At this point we introduce some more notations. Let 𝔤𝔩(n,ℝ) denote the vector space of all n×n real matrices equipped with the metric induced by the Frobenius norm; that is,
(11)〈X,Y〉:=Tr(X⊤Y),X,Y∈𝔤𝔩(n,ℝ).
It turns out that 𝔤𝔩(n,ℝ) is the Lie algebra of the general linear group GL(n,ℝ), consisting of all n×n invertible matrices with real entries.Consider the natural action of the Lie groupGL(n,ℝ) on its Lie algebra:
(12)σ:GL(n,ℝ)×𝔤𝔩(n,ℝ)⟶𝔤𝔩(n,ℝ)(A,X)↦σ(A,X)=AXA⊤.
The orbit of the identity matrix I with respect to σ,
(13)𝒪I={AA⊤:A∈GL(n,ℝ)},
is the set of all SPD matrices, hereafter denoted by 𝔰+(n). With this characterization, it follows immediately that one may look at 𝔰+(n) as an embedding manifold of the manifold consisting of all symmetric matrices of order n, 𝔰(n).In order to characterize the tangent space of𝔰+(n) at a point P,  TP𝔰+(n), one considers any smooth curve t↦α(t)∈𝔰+(n) satisfying α(0)=P and derives conditions for α˙(0). It can be easily checked that
(14)TP𝔰+(n)={P1/2SP1/2:S∈𝔰(n)},
which is isomorphic to 𝔰(n).The Frobenius inner product in𝔰(n) induces the following inner product in TP𝔰+(n):
(15)〈X,Y〉P=Tr(P-1XP-1Y),∀X,Y∈TP𝔰+(n).
Endowed with the metric (15), 𝔰+(n) has the structure of an embedded Riemannian submanifold of 𝔰(n).SinceTP𝔰+(n)≅𝔰(n), it is immediate to conclude that the orthogonal space TP⊥𝔰+(n) consists of the null matrix only.The unique geodesicγ, with respect to the metric (15), satisfying γ(0)=P∈𝔰+(n) and (dγ/dt)(0)=W∈TP𝔰+(n) is given explicitly by
(16)γ(t)=P1/2etP-1/2WP-1/2P1/2,
where eX is the usual matrix exponential.According to the above expression for geodesics, the geodesic distance between two SPD matricesP and Q is simply given by
(17)d(P,Q)=〈log(P-1/2QP-1/2),log(P-1/2QP-1/2)〉1/2.
Since the exponential map e:𝔰(n)→𝔰+(n) is a primary matrix function, [16], it satisfies
(18)P-1eXP=eP-1XP,
and, therefore, the parameterized geodesic (16) may be rewritten as
(19)γ(t)=PetP-1W.
Moreover, the exponential map is a global diffeomorphism, and thus its inverse, denoted by log, is defined in the whole space 𝔰+(n).Analogously, the matrix logarithm also satisfies(20)log(PAP-1)=P(logA)P-1,
for nonsingular matrices A and P.The above property together with the expression for the derivative of the exponential given in Lemma1 plays an important role in the results derived in this paper.Before stating the next result, define the adjoint operator atX∈𝔤𝔩(n,ℝ) by
(21)adX:𝔤𝔩(n,ℝ)⟶𝔤𝔩(n,ℝ)Y↦adX(Y)=[X,Y]=XY-YX,
and note that
(22)eXYe-X=eadX(Y).Lemma 1 (see [17]).
Ift↦A(t) is a matrix valued function taking values in 𝔰(n), then
(23)ddteA(t)=eu-1u|u=adA(t)(A˙(t))eA(t),
where A˙=dA/dt and the operator (eu-1)/u denotes the sum of the power series ∑m=0+∞‍um/(m+1)!.The next two lemmas contain properties of the operator(eu-1)/u that will be used to derive our main result.Lemma 2.
IfA is a symmetric matrix, then ((eu-1)/u)|u=adA is a self-adjoint operator.Proof.
According to the Taylor series expansion of the operator(eu-1)/u, it is enough to prove that the operator adAm is self-adjoint, for all m∈ℕ; that is,
(24)〈adAm(X),Y〉=〈X,adAm(Y)〉,
for all X,Y∈𝔤l(n,ℝ). This fact can be proved using induction on m.
Whenm=1, using the properties of the trace, one gets
(25)〈adA(X),Y〉=Tr([A,X]⊤Y)=Tr((X⊤A-AX⊤)Y)=Tr(X⊤AY-X⊤YA)=Tr(X⊤[A,Y])=〈X,adA(Y)〉.
Now, assume that (24) is satisfied for a positive integer m and prove that it is also satisfied for m+1. To do that, one uses the fact adAm+1=adA(adAm)=adAm(adA), the basic step, and also the inductive step:
(26)〈adAm+1(X),Y〉=〈adA(adAm)(X),Y〉=〈adAm(X),adA(Y)〉=〈X,adAm(adA(Y))〉=〈X,adAm+1(Y)〉.Lemma 3.
LetP∈𝔰+(n) and X,Z∈𝔰(n). Then the following holds:
(27)eu-1u|u=adP-1X(Z)=P-1/2eu-1u|u=adP-1/2XP-1/2(P1/2ZP-1/2)P1/2.Proof.
Again, taking into account the Taylor series expansion of the operator(eu-1)/u, it is enough to prove that
(28)adP-1Xm(Z)=P-1/2adP-1/2XP-1/2m(P1/2ZP-1/2)P1/2,
for all positive integers m.
It is immediate to check that assertion (28) holds for m=1. In fact,
(29)P-1/2adP-1/2XP-1/2(P1/2ZP-1/2)P1/2=P-1/2P-1/2XP-1/2P1/2ZP-1/2P1/2-P-1/2P1/2ZP-1/2P-1/2XP-1/2P1/2=P-1XZ-ZP-1X=adP-1X(Z).
To prove the induction step, assume that (28) is satisfied for m and prove that it is also satisfied for m+1:
(30)adP-1Xm+1(Z)=adP-1X(adP-1Xm(Z))=adP-1X(P-1/2adP-1/2XP-1/2m(P1/2ZP-1/2)P1/2)=P-1/2adP-1/2XP-1/2(adP-1/2XP-1/2m(P1/2ZP-1/2))P1/2=P-1/2adP-1/2XP-1/2m+1(P1/2ZP-1/2)P1/2.To finish this section, we state in the next lemma one interesting property concerning the derivative of the matrix logarithm function.Lemma 4 (see [15]).
Lett↦P(t) be a matrix valued function taking values in 𝔰+(n). Then,
(31)ddtTr(log2P(t))=2Tr((logP(t))P-1(t)P˙(t)).


### 3.1. Normal Equations for the Geodesic Fitting Problem

In order to derive the first order necessary optimality conditions for the optimization problem stated in Section2.1 for the particular case of the SPD manifold 𝔰+(n), equipped with the Riemannian metric (15), we start by finding the critical points of the function
(32)F:𝔰+(n)×𝔰(n)⟶ℝ,(P,X)↦F(P,X)=12∑i=0N‍Tr(log2(Pi-1/2PetiP-1XPi-1/2)).
For this particular purpose, let s↦P(s) be a smooth curve in 𝔰+(n) satisfying P(0)=P and P˙(0)=Y, and let X:s↦X(s) be a smooth curve in 𝔰(n) satisfying X(0)=X and X˙(0)=Z. Then, (P,X) is a critical point of F if and only if
(33)(dF)(P,X)(Y,Z)=0,∀Y,Z∈𝔰(n).
In order to simplify the writing, let us define
(34)Wi(s)=Pi-1/2P(s)etiP(s)-1X(s)Pi-1/2.
Then, using Lemma 4,
(35)(dF)(P,X)(Y,Z)=dds|s=0F(P(s),X(s))=∑i=0N‍Tr((logWi(0))Wi-1(0)W˙i(0)).
In order to compute W˙i(0), we use the expression for the derivative of the matrix exponential given by (23) to write
(36)W˙i(0)=Pi-1/2YetiP-1XPi-1/2+Pi-1/2Peu-1u|u=adtiP-1X(-tiP-1YP-1X+tiP-1Z)etiP-1XPi-1/2.
Plugging the above expression into (35), one gets
(37)(dF)(P,X)(Y,Z)=∑i=0N‍〈W˙i(0)Wi-1(0),logWi(0)〉=∑i=0N〈Pi-1/2(Y-Peu-1u|u=adtiP-1X(tiP-1YP-1X))‍×P-1Pi1/2,logWi(0)(Y-Peu-1u|u=adtiP-1X(tiP-1YP-1X))〉+∑i=0N‍〈Pi-1/2Peu-1u|u=adtiP-1X(tiP-1Z)P-1Pi1/2,logWi(0)eu-1u|u=adtiP-1X〉=∑i=0N‍〈(Y-Peu-1u|u=adtiP-1X(tiP-1YP-1X))×P-1,Pi-1/2logWi(0)Pi1/2eu-1u|u=adtiP-1/2XP-1/2〉+∑i=0N‍〈Peu-1u|u=adtiP-1X(tiP-1Z)P-1,Pi-1/2logWi(0)Pi1/2eu-1u|u=adtiP-1X〉.
We now state the main result of the paper.Theorem 5.
A necessary condition for the geodesicγ(t)=P1/2etP-1/2XP-1/2P1/2 to be a local minimizer of the function F defined by (32) is that the pair (P,X)∈𝔰+(n)×𝔰(n) satisfies the following system of equations:
(38)∑i=0N‍log(γ(ti)Pi-1)=0,∑i=0N‍tieu-1u|u=adtiP-1/2XP-1/2(P1/2log(Pi-1γ(ti))P-1/2)=0.Proof.
By definition,(P,X) is a critical point of F if and only if (dF)(P,X)(Y,Z)=0, for all Y,Z∈𝔰(n). Therefore, according to (37), one must have
(39)∑i=0N‍〈Peu-1u|u=adtiP-1X(tiP-1Z)P-1,Pi1/2logWi(0)Pi-1/2eu-1u|u=adtiP-1X〉=0,
for all Z∈𝔰(n), and, simultaneously,
(40)∑i=0N‍〈(Y-Peu-1u|u=adtiP-1X(tiP-1YP-1X))×P-1,Pi1/2logWi(0)Pi-1/2eu-1u|u=adtiP-1/2XP-1/2〉=0,
for all Y∈𝔰(n).
Using Lemma3  and the fact that the operator (eu-1)/u is self-adjoint when restricted to 𝔰(n), (39) can be rewritten as
(41)∑i=0N‍〈eu-1u|u=adtiP-1/2XP-1/2(tiP-1/2ZP-1/2),P1/2Pi-1/2logWi(0)Pi1/2P-1/2eu-1u|u=adtiP-1/2XP-1/2〉=0⟺〈∑i=0N‍tieu-1u|u=adtiP-1/2XP-1/2P-1/2ZP-1/2,∑i=0N‍tieu-1u|u=adtiP-1/2XP-1/2(P1/2Pi-1/2×logWi(0)Pi1/2P-1/2)eu-1u|u=adtiP-1/2XP-1/2〉=0⟺〈∑i=0N‍tieu-1u|u=adtiP-1/2XP-1/2P1/2ZP1/2,∑i=0N‍tieu-1u|u=adtiP-1/2XP-1/2(P1/2Pi-1/2×logWi(0)Pi1/2P-1/2)eu-1u|u=adtiP-1/2XP-1/2〉P=0.
But, since TP⊥𝔰+(n)={0}, one must have
(42)∑i=0N‍tieu-1u|u=adtiP-1/2XP-1/2(P1/2Pi-1/2logWi(0)Pi1/2P-1/2)=0.
Now, using the definition of Wi(s) given by (34) and the logarithm property (20), it is immediate to check that (42) is equivalent to the second equation in (38).
Using similar arguments, (40) is equivalent to
(43)〈P-1/2YP-1/2,∑i=0N‍P1/2Pi1/2logWi(0)Pi-1/2P-1/2〉=〈∑i=0N‍tieu-1u|u=adtiP-1/2XP-1/2P-1/2YP-1XP-1/2,∑i=0N‍tieu-1u|u=adtiP-1/2XP-1/2(P1/2Pi-1/2×logWi(0)Pi1/2P-1/2)eu-1u|u=adtiP-1/2XP-1/2〉=0.
Taking into account (42), the last identity reduces to
(44)〈P-1/2YP-1/2,∑i=0N‍P1/2Pi1/2logWi(0)Pi-1/2P-1/2〉=0.
But,
(45)〈P-1/2YP-1/2,∑i=0N‍P1/2Pi1/2logWi(0)Pi-1/2P-1/2〉=0⟺〈P1/2YP1/2,∑i=0N‍P1/2Pi1/2logWi(0)P1/2YP1/2,×Pi-1/2P-1/2∑i=0N‍〉P=0,
meaning that
(46)∑i=0N‍P1/2Pi1/2logWi(0)Pi-1/2P-1/2=0.
Using similar arguments to those used above, it can be easily checked that (46) is equivalent to the first equation of (38), and this completes the proof.Remark 6.
Similarly to the Euclidean counterpart, we name (38) the normal equations for the geodesic fitting problem in 𝔰+(n). However, in contrast to the classical case where the corresponding normal equations can be solved explicitly, the normal equations for the geodesic fitting problem on the SPD manifold are highly non-linear, and only approximate solutions are expected to be obtained using appropriate numerical methods.


## 3.1. Normal Equations for the Geodesic Fitting Problem

In order to derive the first order necessary optimality conditions for the optimization problem stated in Section2.1 for the particular case of the SPD manifold 𝔰+(n), equipped with the Riemannian metric (15), we start by finding the critical points of the function
(32)F:𝔰+(n)×𝔰(n)⟶ℝ,(P,X)↦F(P,X)=12∑i=0N‍Tr(log2(Pi-1/2PetiP-1XPi-1/2)).
For this particular purpose, let s↦P(s) be a smooth curve in 𝔰+(n) satisfying P(0)=P and P˙(0)=Y, and let X:s↦X(s) be a smooth curve in 𝔰(n) satisfying X(0)=X and X˙(0)=Z. Then, (P,X) is a critical point of F if and only if
(33)(dF)(P,X)(Y,Z)=0,∀Y,Z∈𝔰(n).
In order to simplify the writing, let us define
(34)Wi(s)=Pi-1/2P(s)etiP(s)-1X(s)Pi-1/2.
Then, using Lemma 4,
(35)(dF)(P,X)(Y,Z)=dds|s=0F(P(s),X(s))=∑i=0N‍Tr((logWi(0))Wi-1(0)W˙i(0)).
In order to compute W˙i(0), we use the expression for the derivative of the matrix exponential given by (23) to write
(36)W˙i(0)=Pi-1/2YetiP-1XPi-1/2+Pi-1/2Peu-1u|u=adtiP-1X(-tiP-1YP-1X+tiP-1Z)etiP-1XPi-1/2.
Plugging the above expression into (35), one gets
(37)(dF)(P,X)(Y,Z)=∑i=0N‍〈W˙i(0)Wi-1(0),logWi(0)〉=∑i=0N〈Pi-1/2(Y-Peu-1u|u=adtiP-1X(tiP-1YP-1X))‍×P-1Pi1/2,logWi(0)(Y-Peu-1u|u=adtiP-1X(tiP-1YP-1X))〉+∑i=0N‍〈Pi-1/2Peu-1u|u=adtiP-1X(tiP-1Z)P-1Pi1/2,logWi(0)eu-1u|u=adtiP-1X〉=∑i=0N‍〈(Y-Peu-1u|u=adtiP-1X(tiP-1YP-1X))×P-1,Pi-1/2logWi(0)Pi1/2eu-1u|u=adtiP-1/2XP-1/2〉+∑i=0N‍〈Peu-1u|u=adtiP-1X(tiP-1Z)P-1,Pi-1/2logWi(0)Pi1/2eu-1u|u=adtiP-1X〉.
We now state the main result of the paper.Theorem 5.
A necessary condition for the geodesicγ(t)=P1/2etP-1/2XP-1/2P1/2 to be a local minimizer of the function F defined by (32) is that the pair (P,X)∈𝔰+(n)×𝔰(n) satisfies the following system of equations:
(38)∑i=0N‍log(γ(ti)Pi-1)=0,∑i=0N‍tieu-1u|u=adtiP-1/2XP-1/2(P1/2log(Pi-1γ(ti))P-1/2)=0.Proof.
By definition,(P,X) is a critical point of F if and only if (dF)(P,X)(Y,Z)=0, for all Y,Z∈𝔰(n). Therefore, according to (37), one must have
(39)∑i=0N‍〈Peu-1u|u=adtiP-1X(tiP-1Z)P-1,Pi1/2logWi(0)Pi-1/2eu-1u|u=adtiP-1X〉=0,
for all Z∈𝔰(n), and, simultaneously,
(40)∑i=0N‍〈(Y-Peu-1u|u=adtiP-1X(tiP-1YP-1X))×P-1,Pi1/2logWi(0)Pi-1/2eu-1u|u=adtiP-1/2XP-1/2〉=0,
for all Y∈𝔰(n).
Using Lemma3  and the fact that the operator (eu-1)/u is self-adjoint when restricted to 𝔰(n), (39) can be rewritten as
(41)∑i=0N‍〈eu-1u|u=adtiP-1/2XP-1/2(tiP-1/2ZP-1/2),P1/2Pi-1/2logWi(0)Pi1/2P-1/2eu-1u|u=adtiP-1/2XP-1/2〉=0⟺〈∑i=0N‍tieu-1u|u=adtiP-1/2XP-1/2P-1/2ZP-1/2,∑i=0N‍tieu-1u|u=adtiP-1/2XP-1/2(P1/2Pi-1/2×logWi(0)Pi1/2P-1/2)eu-1u|u=adtiP-1/2XP-1/2〉=0⟺〈∑i=0N‍tieu-1u|u=adtiP-1/2XP-1/2P1/2ZP1/2,∑i=0N‍tieu-1u|u=adtiP-1/2XP-1/2(P1/2Pi-1/2×logWi(0)Pi1/2P-1/2)eu-1u|u=adtiP-1/2XP-1/2〉P=0.
But, since TP⊥𝔰+(n)={0}, one must have
(42)∑i=0N‍tieu-1u|u=adtiP-1/2XP-1/2(P1/2Pi-1/2logWi(0)Pi1/2P-1/2)=0.
Now, using the definition of Wi(s) given by (34) and the logarithm property (20), it is immediate to check that (42) is equivalent to the second equation in (38).
Using similar arguments, (40) is equivalent to
(43)〈P-1/2YP-1/2,∑i=0N‍P1/2Pi1/2logWi(0)Pi-1/2P-1/2〉=〈∑i=0N‍tieu-1u|u=adtiP-1/2XP-1/2P-1/2YP-1XP-1/2,∑i=0N‍tieu-1u|u=adtiP-1/2XP-1/2(P1/2Pi-1/2×logWi(0)Pi1/2P-1/2)eu-1u|u=adtiP-1/2XP-1/2〉=0.
Taking into account (42), the last identity reduces to
(44)〈P-1/2YP-1/2,∑i=0N‍P1/2Pi1/2logWi(0)Pi-1/2P-1/2〉=0.
But,
(45)〈P-1/2YP-1/2,∑i=0N‍P1/2Pi1/2logWi(0)Pi-1/2P-1/2〉=0⟺〈P1/2YP1/2,∑i=0N‍P1/2Pi1/2logWi(0)P1/2YP1/2,×Pi-1/2P-1/2∑i=0N‍〉P=0,
meaning that
(46)∑i=0N‍P1/2Pi1/2logWi(0)Pi-1/2P-1/2=0.
Using similar arguments to those used above, it can be easily checked that (46) is equivalent to the first equation of (38), and this completes the proof.Remark 6.
Similarly to the Euclidean counterpart, we name (38) the normal equations for the geodesic fitting problem in 𝔰+(n). However, in contrast to the classical case where the corresponding normal equations can be solved explicitly, the normal equations for the geodesic fitting problem on the SPD manifold are highly non-linear, and only approximate solutions are expected to be obtained using appropriate numerical methods.


## 4. Concluding Remarks

In this paper we studied the geodesic fitting problem for the SPD manifold endowed with the canonical Riemannian metric (15). This optimization problem is the natural generalization of the classical linear regression problem. The adopted metric has the particularity of being the unique metric that is GL(n,ℝ)-invariant. Also, the sectional curvature of this metric is nonnegative, an extremely convenient geometric feature since, in this case, the square of the distance function has the particularity of being convex.The optimization problem studied here turned out to be more complex than the corresponding problem for the SPD manifold endowed with the Log-Euclidean metric. This metric was defined in [18] as being the Euclidean metric in the space of the logarithms of SPD matrices. With respect to the Log-Euclidean metric, the geodesic fitting problem could be solved explicitly in [13].In the main result of this paper, Theorem5, we have derived the counterpart of the normal equations of the classical linear regression problem. However, contrary to the classical situation, we obtained highly nonlinear equations that require numerical methods on manifolds in order to obtain approximate solutions. These numerical issues will be analysed in future work.Generalizations of other classical curve fitting methods, such as polynomial or spline approximation, are also important but very challenging. These problems started to gain ground since the pioneer work of Jupp and Kent [19], for spherical smoothing splines. Since then, there is an increasing interest in this subject, but the results are not yet satisfactory. One of the drawbacks is that no explicit formulas for generalized polynomials on manifolds are known. The variational approach appearing in [6, 20–22] and the analytical approach considered in [23] to generalize polynomials on manifolds end up facing the problem of solving highly nonlinear differential equations. Nevertheless, a variational approach to solve more general least square problems on manifolds that appeared in [21] may inspire some readers, at least theoretically.



---

*Source: 425608-2013-08-26.xml*