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.

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 . 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 Section 2, 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, [810].

Let be an dimensional Riemannian manifold with Riemannian metric and its Levi-Civita connection. Given , denotes the tangent space of at , and stands for the tangent bundle of .

A vector field along a curve is a mapping that assigns to each , the vector . If is smooth, the velocity vector field, , is an example of a smooth vector field in . If is induced by the vector field , that is, if , then the covariant derivative of along is defined by A vector field along a curve is said to be parallel if By definition, a geodesic is a smooth curve whose velocity vector field is parallel along . That is, This is a second order differential equation and according to the theory of existence and uniqueness for ordinary differential equations, given and , there exists a unique geodesic , satisfying and . This geodesic can be parameterized explicitly by where is the point on the geodesic that is at a distance equal to away from .

In spite of the fact that the exponential map, , 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 with initial velocity are the straight lines: For the unit sphere , equipped with the Riemannian metric induced from the Euclidean metric of the embedding space, geodesics are great circles. In particular, the geodesic starting at with velocity is given by The manifold is said to be geodesically complete if every pair of points in can be joined by a geodesic. In this case, (see, for instance [11]) the geodesic that joins at to at is parameterized explicitly by and the distance between and , which is simply the length of this minimizing geodesic, is given by

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 in , , and a set of instants of time , for simplicity assumed to form a partition of the unit time interval , find a parameterized geodesic that yields the minimum value for the functional where 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 when is the unit 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 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 denote the vector space of all real matrices equipped with the metric induced by the Frobenius norm; that is, It turns out that is the Lie algebra of the general linear group , consisting of all invertible matrices with real entries.

Consider the natural action of the Lie group on its Lie algebra: The orbit of the identity matrix with respect to , is the set of all SPD matrices, hereafter denoted by . With this characterization, it follows immediately that one may look at as an embedding manifold of the manifold consisting of all symmetric matrices of order , .

In order to characterize the tangent space of at a point , , one considers any smooth curve satisfying and derives conditions for . It can be easily checked that which is isomorphic to .

The Frobenius inner product in induces the following inner product in : Endowed with the metric (15), has the structure of an embedded Riemannian submanifold of .

Since , it is immediate to conclude that the orthogonal space consists of the null matrix only.

The unique geodesic , with respect to the metric (15), satisfying and is given explicitly by where is the usual matrix exponential.

According to the above expression for geodesics, the geodesic distance between two SPD matrices and is simply given by Since the exponential map is a primary matrix function, [16], it satisfies and, therefore, the parameterized geodesic (16) may be rewritten as Moreover, the exponential map is a global diffeomorphism, and thus its inverse, denoted by , is defined in the whole space .

Analogously, the matrix logarithm also satisfies for nonsingular matrices and .

The above property together with the expression for the derivative of the exponential given in Lemma 1 plays an important role in the results derived in this paper.

Before stating the next result, define the adjoint operator at by and note that

Lemma 1 (see [17]). If is a matrix valued function taking values in , then where and the operator denotes the sum of the power series .

The next two lemmas contain properties of the operator that will be used to derive our main result.

Lemma 2. If is a symmetric matrix, then is a self-adjoint operator.

Proof. According to the Taylor series expansion of the operator , it is enough to prove that the operator is self-adjoint, for all ; that is, for all . This fact can be proved using induction on .
When , using the properties of the trace, one gets Now, assume that (24) is satisfied for a positive integer and prove that it is also satisfied for . To do that, one uses the fact , the basic step, and also the inductive step:

Lemma 3. Let and . Then the following holds:

Proof. Again, taking into account the Taylor series expansion of the operator , it is enough to prove that for all positive integers .
It is immediate to check that assertion (28) holds for . In fact, To prove the induction step, assume that (28) is satisfied for and prove that it is also satisfied for :

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]). Let be a matrix valued function taking values in . Then,

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 Section 2.1 for the particular case of the SPD manifold , equipped with the Riemannian metric (15), we start by finding the critical points of the function For this particular purpose, let be a smooth curve in satisfying and , and let be a smooth curve in satisfying and . Then, is a critical point of if and only if In order to simplify the writing, let us define Then, using Lemma 4, In order to compute , we use the expression for the derivative of the matrix exponential given by (23) to write Plugging the above expression into (35), one gets We now state the main result of the paper.

Theorem 5. A necessary condition for the geodesic to be a local minimizer of the function defined by (32) is that the pair satisfies the following system of equations:

Proof. By definition, is a critical point of if and only if , for all . Therefore, according to (37), one must have for all , and, simultaneously, for all .
Using Lemma 3  and the fact that the operator is self-adjoint when restricted to , (39) can be rewritten as But, since , one must have Now, using the definition of 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 Taking into account (42), the last identity reduces to But, meaning that 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 . 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 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, Theorem 5, 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, 2022] 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.

Acknowledgment

L. Machado and F. Silva Leite work is partially supported by Project PTDC/EEA-CRO/122812/2010.