This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2496166
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, FEBRUARY 2015

1

Principal Curves on Riemannian Manifolds
Søren Hauberg
Abstract—Euclidean statistics are often generalized to Riemannian manifolds by replacing straight-line interpolations withgeodesic ones. While these Riemannian models are familiar-looking, they are restricted by the inflexibility of geodesics, and they
rely on constructions which are optimal only in Euclidean domains. We consider extensions of Principal Component Analysis(PCA) to Riemannian manifolds. Classic Riemannian approaches seek a geodesic curve passing through the mean that optimize
a criteria of interest. The requirements that the solution both is geodesic and must pass through the mean tend to imply that themethods only work well when the manifold is mostly flat within the support of the generating distribution. We argue that instead ofgeneralizing linear Euclidean models, it is more fruitful to generalize non-linear Euclidean models. Specifically, we extend theclassic Principal Curves from Hastie & Stuetzle to data residing on a complete Riemannian manifold. We show that for elliptical
distributions in the tangent of spaces of constant curvature, the standard principal geodesic is a principal curve. The proposedmodel is simple to compute and avoids many of the pitfalls of traditional geodesic approaches. We empirically demonstrate theeffectiveness of the Riemannian principal curves on several manifolds and datasets.
Index Terms—Principal component analysis, Principal curves, Differential geometry, Riemannian manifolds

F

1

I NTRODUCTIONP

RINCIPAL COMPONENT ANALYSIS ( PCA ) is a core tool
for statistical analysis of vectorial data in RD . Given adata set x1:N = {x1 , . . . , xN }, PCA estimates the subspace
through the mean that represent the data as well as possible,

vPCA = arg minNXv, kvk=1 n=1→ (xn ) − xn k2 ,k proj−µv

NXv, kvk=1 n=1→ (xn ) − µk2 .
k proj−µv(2)While PCA is an essential tool for data analysis in RD ,not all data is Euclidean. Attention is growing towards data
residing on Riemannian manifolds, e.g. shape data [1]–[6], DTI images [7]–[9], image features [10]–[12], motion
models [13], [14], human poses [15]–[18], and more.
Statistical models on Riemannian manifolds often generalize Euclidean counterparts by replacing straight lines withgeodesic curves in the problem formulations. The resulting
models are generally referred to as Principal Geodesic
Analysis (PGA). In this line of thinking, Eq. 1 becomesvPGA1 = arg min

NXv∈Tµ M n=1


dist2 projγµ (v) (xn ), xn ,Poincaré half-plane

Example geodesic

Poincaré half-plane(1)→ denotes projection onto the line spanned bywhere proj−µv
the data mean µ and principal component v. Expanding
the squared norm shows that minimizing the reconstruction
error (1) is equivalent to maximizing the projected variancevPCA = arg max

Data
Intrinsic mean
Principal geodesic
Principal curve

Fig. 1. Left: An initial example: observations (bluedots) are given in the Poincaré half-plane. The intrinsic
mean (red diamond) and the principal geodesic (yellow
dashed line) both do a poor job of summarizing themain trend of the data due to the restrictive structureof geodesic curves. The principal curve (purple line),
in contrast, captures the main mode of variation well.
Right: Example geodesics in the Poincaré half-plane
illustrate the geometry of this space; geodesics areeither vertical lines or half-circles.
distance; and projc denotes projection onto the curve c:
projc (x) = c(t̂(x))

where
t̂(x) = arg min dist(c(t), x) .
t

Similarly, Eq. 2 becomesvPGA2 = arg max(3)

where γµ (v) : t 7→ Expµ (t · v) denotes the geodesic withstarting point µ and initial velocity v; dist is the geodesic• S. Hauberg is with the Section for Cognitive Systems at the Technical
University of Denmark.
E-mail: sohau@dtu.dk

(4)

N
X

v∈Tµ M n=1


dist2 projγµ (v) (xn ), µ .

(5)

It is worth noting that while Eq. 1 and 2 are equivalent
in the Euclidean domain, the corresponding Riemannian
energies (3, 5) may differ in solutions [19]. We shall returnto these models in Sec. 2; for now it suffice to think of these
Riemannian generalizations to PCA as replacing straightlines with geodesics to give methods in non-linear spaces.
Figure 1 provides an initial example with data in the
Poincaré half-plane. PGA is unable to fit the main trend of

Copyright (c) 2015 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2496166
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, FEBRUARY 2015

the data for two reasons: 1) The solution curve is required
to pass through the data mean, which is outside the supportof the data; and 2) the solution is required to be a geodesic,
but no suitable geodesic exists in the Poincaré half-plane.
We argue for an alternative: Rather than extending linear
Euclidean models, we find it simpler and more flexible
to extend non-linear Euclidean techniques to Riemanniandomains. In particular, this allows us to drop the two above
assumptions. Figure 1 shows that this improves upon PGA.
This paper is organized as follows: The next section
briefly introduce Riemannian geometry and related statisticaloperations. Here notation is also established. Section 3
describe principal curves and the proposed Riemannian
extension. Experimental results are provided in Sec. 4, andthe paper is concluded with a discussion in Sec. 5. The
online supplements contain both proofs and animations.

2

BACKGROUND AND R ELATED W ORK

In this section, we establish elementary prerequisites. Moredetails can be found in standard references, e.g. [20]–[22].

2where wn denotes optional weights for each observation. In
Euclidean domains, this problem is uniquely solved by theempirical average. On Riemannian manifolds the solution
may, however, no longer be unique; e.g., consider uniformly
distributed data on a unit sphere (the entire sphere solves
Eq. 7). Whether µ is unique or not depends on both the
spread of data and the curvature of the manifold [23], [24].
Here, we recall the uniqueness result from Kendall [24].
Definition 1 (Regular geodesic ball). The ball B(x̃, r) ={x ∈ M | dist(x, x̃) < r} is “geodesic” if there exista minimizing geodesic from the√ center x̃ to any pointin B(x̃, r). This is regular if 2r κ < π, where κ is the
maximum of the Riemannian curvature in the ball.
Theorem 1 (Existence and uniqueness of means [24]). Let
the data be drawn from a distribution with support inside
the regular geodesic ball B(x̃, r), then there exist one andonly one intrinsic mean inside B(x̃, r).
In practice, any concern about uniqueness is ignored andEq. 7 is solved using gradient descent; see Algorithm 1 [22].
Empirical experience indicates that multiple modes of Eq. 7
often exist for real data, so it is potentially problematic torely on the uniqueness of a mean.

2.1 Elements of Riemannian Geometry
Riemannian geometry studies smooth manifolds endowedwith a smoothly changing metric. A metric on a manifold
M is an inner product in the tangent space Tx M at each Algorithm 1 Intrinsic mean.
point x ∈ M. If the inner product varies smoothly it isInitialize µ, e.g. as a random data point.
called a Riemannian metric. The tangent space Tx M isrepeat PN
a Euclidean space that locally approximates the manifold,∇ = n=1 wn Logµ (xn ).
and the tangential norm can be integrated to provide curveµ = Expµ (α · ∇).
lengths. Any curve that locally minimizes its length is known
until convergence
as a geodesic. We assume that all points can be joined bya geodesic, i.e. we only consider complete Riemannian
manifolds. Distances between points are then defined as the
2.3 Principal Geodesics
length of the shortest connecting geodesic.
Many operations are defined in the Euclidean tangent As discussed in Sec. 1, standard generalizations of PCA tospace Tx M, and there are mappings between this approxima- Riemannian manifolds seeks a geodesic curve through thetion of the manifold and the manifold itself. The exponential mean that maximizes the projected variance (5) or minimizes
map at a point x takes a tangent vector v ∈ Tx M to
the reconstruction error (3). This assumes a unique mean.
Fletcher et al. [25] suggest an approximate solution toy = Expx (v) ∈ M(6)(3) that computes the principal geodesic as the leading
such that the curve t 7→ Expx (t·v) is a geodesic originating eigenvector of the covariance matrix
at x with initial direction v. In general, the exponential map
N
1 X
is only defined in a neighborhood of x. The inverse mapping,Σ=
Logµ (xn ) Logµ (xn )| .
(8)
N − 1 n=1
which takes y to v ∈ Tx M, is known as the logarithm map
and is denoted Logx (y).
In essence, Fletcher et al. perform Euclidean PCA in thetangent space at the mean. Sommer et al. [26] provide a
2.2 Intrinsic Meansnumerical optimization scheme for solving Eq. 3 directly
Expectations are defined in the Euclidean domain through with no approximations. In practice, this is rarely used as it
integration. This, however, generalize poorly to Riemannian requires numerical integration of Jacobi fields and second
domains as the expectation of a manifold-valued random order derivatives of geodesic families. Fletcher et al. [27]variable may easily violate the constraints of the manifold also consider Eq. 5, and again propose an approximation[8]. It is therefore common to instead consider intrinsic scheme that performs Euclidean PCA in the tangent spacemeans which are defined as the point with minimal sum-of- at the mean. The present author is unaware of any exact
squared distances to the data,
numerical solution to Eq. 5.
N
It is worth noting that both Eq. 3 and 5 rely uponX
µ = arg min
wn dist2 (xn , x),(7) projections of points onto a geodesic curve (4). Thisxn=1
projection need not be unique; e.g. consider projecting the

Copyright (c) 2015 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2496166
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, FEBRUARY 20153

north pole on a sphere onto the equator (all points on theequator are solutions). Practical implementations ignore thispotential issue and perform the projection through a onedimensional optimization of Eq. 4 along the curve.
Huckemann et al. [19] consider Riemannian manifolds
modulo isometric Lie group actions, and study an extension
Fig. 2. Self-consistency (9) implies that the expectationof Eq. 3 where the principal geodesic need not pass throughof all points that project to the same point on a curve
the intrinsic mean. The current numerical implementation ofequals the point of projection.
the resulting geodesic PCA, however, relies on brute-force
techniques that are not expected to scale gracefully.
Algorithm 2 Riemannian principal curve.
Panaretos et al. [28] drop the geodesic requirement andInitialize the principal curve as c1:T = {c(1), . . . , c(T )}.
consider smooth curves where any point on the curve has arepeattangent velocity vector that fits the leading eigenvector oft̂n = arg min
 t (dist(c(t),xn ), ∀n.
a tangent space PCA. This increases flexibility, though thedist(c(t), c(t̂n )w=k, ∀n, t.
t,ncurve is required to pass through the intrinsic mean.
σ
PN
Jung et al. [29] consider data on the unit sphere andc(t) = arg minx n=1 wt,n dist2 (xn , x), ∀t.
recursively discard the least interesting dimension of theuntil convergencedata, thereby getting principal nested spheres. This elegantly
removes the dependency on the intrinsic mean as thecomponents are estimated before the mean. The approach 3.2 Principal Curves for Riemannian Datais, however, dependent on the embedding of the sphere andit is unclear if the approach can generalize to other spaces. Self-consistency (9) is extended to the Riemannian domain
The above methods are all derived from a least-squares by replacing the expectation with an intrinsic mean (7)Nperspective. Recent interest is growing toward more generalXc(t) = arg minδt=t̂n dist2 (xn , x),(10)probabilistic models, such as the Probabilistic PGA fromxn=1
Zhang & Fletcher [30]. Here a maximum likelihood estimateof a tangential normal distribution is estimated using a where δt=t̂ is the Kronecker delta function, i.e. δt=t̂ =nnMonte Carlo EM algorithm. Attractively, this provides a 1 for t = t̂n and 0 otherwise. Note that this definition
noise model of the data, but the reliance on Monte Carlo is not enough to ensure a continuous principal curve as
makes the approach computationally challenging.
intrinsic means need not change continuously when the datachanges continuously. In practice, the issue is resolved byrepresenting the curve as a discrete set of points.
3 P RINCIPAL C URVES
From the intrinsic self-consistency (10), Hastie & StuetIn this section we first review the classical Euclidean zle’s simple algorithm generalizes directly to data residing
principal curves for non-linear PCA. This is then extended to on Riemannian manifolds. Next we describe the needed
Riemannian manifolds in Sec. 3.2, and some links between modifications, which are summarized in Algorithm 2.
the proposed model and PGA are established in Sec. 3.3.
3.1

Principal Curves for Euclidean Data

Informally, a principal curve in RD is a smooth curve
passing through the “middle” of a distribution [31]. Let
c : R → RD be a one-dimensional curve, and let t̂(x) ∈ R
be the index of the projection of x ∈ RD onto the curvec (see Eq. 4). A principal curve is then defined as a curvethat satisfies the self-consistency property (see Fig. 2)c(t) = E(x | t̂(x) = t).

(9)

3.2.1 Curve Representation
Like Hastie & Stuetzle [31], we represent a curve as aset of points, c1:T = {c(1), . . . , c(T )}, joined by geodesic
segments. This does, in principle, not give a smooth curve1 ,but a fine discretization alleviates this in practice.
A point xn can be projected onto such a discrete curveby computing the geodesic distance to all discrete points in
c1:T , and picking the closest point as the projection, i.e.
projc (xn ) = c(t̂n )

where

(11)
t̂n = arg min(dist(c(t), xn ).
t=1,...,TIn words: If we consider all points that project onto c(t)then their expectation must equal c(t).
Hastie & Stuetzle [31] suggest a simple iterative algorithm 3.2.2 Kernel Smoothing
for computing a principal curve from finite data. The The basic algorithm for computing principal curves firstalgorithm represents the principal curve as a discrete set of projects the data x1:N onto the discrete curve c1:T , c.f .
points, and proceeds in two steps: first the data is projected Eq. 11. Each point c(t) on the curve is then updated as theonto the current curve estimate, and then each discrete point average of the points that project onto the point, c.f . Eq. 10.
on the curve is updated to the average of the data which
1. This is also the case for geodesics on manifolds where geodesics
projects onto the point. To ensure that the resulting curve is must be estimated numerically; e.g. Euler’s method provides a polygonalsmooth, a smoothing kernel is used to estimate the average. curve which must be smoothed post hoc to ensure a smooth geodesic.
Copyright (c) 2015 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2496166
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, FEBRUARY 2015

4This intrinsic mean can be computed numerically using
Algorithm 1. Unsurprisingly, this provides a poor estimationof the principal curve as not enough data projects onto
individual points c(t) to provide a stable average. We, thus,
replace δ in Eq. 10 with a smooth kernel k,probabilistic definition of principal curves that solves thisissue by modeling the data distribution as a mixture-ofGaussians. However, standard generalizations of Gaussian
distributions to Riemannian manifolds [22] does not have
known normalization constants so mixture models arechallenging. We, thus, stay with the original definition (9).
NX
2The original principal curves [31] are critical paths of thec(t) ← arg min
wt,n dist (xn , x),xenergymeasuring sum-of-square-deviations from the data.
n=1
(12)


Duchamp& Stuetzle [34], however, showed that they aredist(c(t), c(t̂n )wt,n = k.
only saddle points and not minimizers. Gerber & Whitakerσ[35] provide an excellent discussion of this phenomena andWe use a quartic kernel [32] k(∆) = (1 − ∆2 )2 · δ|∆|≤1 , show that principal curves as defined by Hastie & Stuetzle
which is smooth and has finite support, but other choices are only well-defined due to implicit regularization, such as
are equally possible.
limiting the number of points used to represent the curve as
well as the kernel smoother used for computing averages.
3.2.3 Closed Principal Curves
These concerns carry over to the Riemannian domain,
Algorithm 2 is an iterative procedure that preserve the where such issues are common; e.g. the principal geodesic
connectivity of the points used to represent the principal is generally not an optimizer of either energies (3, 5) due tocurve. A closed principal curve can, thus, be estimated tangent space approximations. In practice, the main concern
by initializing the algorithm with a closed curve. This is with PGA is the dependency on a global intrinsic mean,in contrast to PGA, which only provides a closed curve, which generally does not exist, and even if it does exist
when geodesics happen to be closed. As we shall see, this may easily fall outside the support of the data distribution.
is a helpful feature for modeling periodic data, which is Principal curves avoid these issues as the finite support ofunavailable in alternative Riemannian models.
the kernel smoother implies that the average in Eq. 12 only
operates on a subset of the data close to c(t), i.e. a local
3.2.4 Initialization
average. If the entire dataset x1:N is in a regular geodesic
Algorithm 2 relies on an initial estimate of the principal ball B(x̃, r), then the data included in the local average
curve. We employ two different strategies, though alterna- (12) is in another regular geodesic ball B(x̂, r̂) ⊆ B(x̃, r)tives are possible. When we seek an open-ended principal of smaller radius (r̂ ≤ r). The local averages needed for
curve, we initialize with the principal geodesic, computed computing a principal curve are then more likely to beusing the approximations from Fletcher et al. [27]. When unique than the global intrinsic average, c.f . Theorem 1.
we seek a closed principal curve, we express the data in its
two leading (Euclidean) principal components in the tangent 4 E XPERIMENTS
space at the intrinsic mean. We fit a circle to this data, andIn this section, we evaluate the empirical performance ofinitialize our algorithm with its image under Expµ (·).
the principal curve on several Riemannian manifolds anddatasets. We compare with the standard principal geodesic
3.3 Relationship Between Principal Geodesics
computed via the leading eigenvector of the covarianceand Principal Curves
matrix in the tangent space of the intrinsic mean [27]. In
For Euclidean data, principal components are only known all experiments, the smoothing parameter is cross-validated.
to be principal curves when the data follows an elliptical
distribution [31]. We provide a similar result in the classic
4.1 Poincaré Half-Plane
model spaces of constant curvature.
As a first illustration, we consider synthetic data distributed
Theorem 2. Let x ∈ Tµ M, where M is either the Euon the two-dimensional Poincaré half-plane H2 [21]. This
clidean plane, the D-sphere or the hyperbolic plane, follow
Riemannian manifold has constant negative curvature, anda zero-mean elliptical distribution with finite covariancea Riemannian metric given by y −1 I at (x, y) ∈ H2 , y > 0.
Σ, and let v1 , . . . , vD be the eigenvectors of Σ in order
Geodesics in H2 are either vertical lines or circular arcsof decreasing eigenvalues. Then, the principal geodesic
perpendicular to the x-axis. The Poincaré metric appears inc(t) = Expµ (t · v1 ) is a principal curve.
information geometry as the Fisher information metric in
the space of univariate normal distributions [36].
Proof: See Appendix A in the supplements.
Figure 1 shows 100 samples from

3.4 Concerns and Benefits2x ∼ U[−3,3] ; y ∼ 3/4|x| + N 10, (1/4) ,(13)
While self-consistency is intuitively appealing, the principal
curve does not necessarily go through the “middle” of where UI denotes the uniform distribution on the intervala dataset. Hastie & Stuetzle [31] note that c need not I. As geodesics are generally circular arcs which are
be a principal curve for data xn = c(tn ) + n , when perpendicular to the x-axis, it is not possible to make any
the noise n is zero-mean. Tibshirani [33] proposes a geodesic follow the trend of the generating distribution (13).
Copyright (c) 2015 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2496166
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, FEBRUARY 2015

As a consequence, the principal geodesic provides a poor fit
to the data as is evident from Fig. 1. It is worth noting thatthe intrinsic mean is not within the support of the distribution(13), and no smooth curve through this point will be able
to capture the trend. In contrast, the principal curve follows
the underlying trend of the generating distribution.
We repeat the experiment for different distributions:2x ∼ U[−3,3] ; y ∼ N 10, (1/4)
(14)


2
x ∼ U[−3,3] ; y ∼ N 5, (1/4) − sin(x)(15)
x ∼ U[−3,3] ;y ∼ 3/4 · x2 + N (10, 1) .

(16)

Figure 3 shows the data along with the principal geodesic
and the principal curve: it is evident that the first fail at
capturing the main trend, while the latter succeeds.
These synthetic experiments illustrate a practical issue
with the principal geodesic construction. Just as straightlines often provide a poor fit to data in Euclidean domains,
geodesics cannot, in general, capture the main trend of agenerating distribution. Even when a geodesic exists that
captures the trend, it may not pass through the intrinsic
mean, which is a practical requirement in PGA. For any but
the simplest generating distributions, the principal geodesiccan, thus, not be expected to capture the main trend.

4.3

Orthogonal Matrices

In the previous section, we only considered a single bonefrom human motion capture data. We consider the same
sequence as before, but now represent the entire skeleton asa collection of 32 rotation matrices describing the rotationof a bone relative to its parent in the kinematic tree [15].
Each data point is, thus, a member of SO(3)32 . We remove
the global position and orientation from the data prior to
analysis. As the underlying walking motion is periodic we
estimate a closed principal curve.
Figure 5 shows example data points along with projections
onto both the principal geodesic and the principal curve.
The supplements contain the corresponding animation. It
is evident that the principal geodesic fails to represent thewalking motion. The principal curve captures the postureand the leg-movement, though the arm-movement is a bit
too simplified. Figure 4 shows the data in the two leading
principal components in the tangent space at the intrinsic
mean. The figure further shows the principal geodesic andthe principal curve projected into this space as well. It isevident that the ability to fit a closed (non-geodesic) curveto the data is ideal for periodic data.
4.4

4.2

5

Smooth Metric Learning

The Unit Sphere

In the previous experiments, there was little a priori reason
Next we consider motion capture data of a person walking to expect a geodesic curve to be able to fit the data, which
in a circular pattern [37], [38]. First, we consider the global gives PGA a disadvantage. To give the principal geodesicorientation of the left femur (i.e. thigh bone) represented ideal conditions, we study data on a Riemannian manifold
as a unit vector, i.e. xn ∈ S2 . The data is periodic and where geodesics by definition follow local trends of theD
represented in a space with periodic geodesics, so the data [40]. Let xn ∈ R denote Euclidean observations. WeD
now equip R with a Riemannian metric by learning metricprincipal geodesic is expected to capture the main trend.
Figure 4 shows that the data approximately follow a small tensors at select locations and smoothly interpolating thesecircle, giving an intrinsic mean at the corresponding pole. to form a Riemannian metric; see [40] for details. HereAny smooth curve passing through this pole will fail to we learn local metrics as local inverse covariance matrices.
model the trend of the data, and the principal geodesic This implies that geodesics locally follow the data.
We consider observations of hand-written one-digits fromdrastically misfits the data. This illustrates a general issue
with PGA. Consider data distributed uniformly over a “band” the MNIST [41] dataset (Fig. 6). For visualization, wearound the equator of a sphere. The intrinsic mean is given project these into the two leading Euclidean principalby either pole, and the principal geodesic will provide a fit components, and learn a smoothly changing metric. The
which is orthogonal to the underlying trend of the data. In left panel of Fig. 7 shows the data along with geodesics
terms of least-squares reconstruction error, the (geodesic) to the intrinsic mean. It is clear that geodesics follow theequator provide a better fit to the data than the principal trend of the data. These are the ideal conditions for PGA.
geodesic. The strong reliance on the intrinsic mean in PGA, The center panel of Fig. 7 shows outgoing geodesics from
thus, prevents the principal geodesic from fitting the data. the intrinsic mean, and the right panel shows the estimatedAs the data captures a periodic phenomena, we fit a principal geodesic. This initially follows the trend of theclosed principal geodesic. The Fig. 4(left) shows that the data, but eventually falls outside the support of the data.
principal curve benefits from not being forced through the The principal curve, again, provides a better summary of
global intrinsic mean. We choose the smoothing parameter the main trend of the data.
σ of the principal curve with cross-validation measuring the
As we have previously seen, one problem for the principal
sum-of-squared distances to the estimated curve on held-out geodesic is that it is sensitive to the location of the intrinsic
data (with a small additive bias towards large values of mean. Here, the mean is slightly outside the support of theσ). Figure 4(center) shows that the error measure has a data, which appears to prevent the principal geodesic fromwell-defined optimum. While cross-validation is known to fitting the data. The sensitivity of PGA to the location of thebe unreliable for principal curves [34], we have found it to mean is concerning, when considering that intrinsic means
work well enough in our experiments. Alternatively model are often not unique. The principal curve avoids this issue
selection could rely on e.g. minimum entropy [39].
to some extent as it only rely on local means.
Copyright (c) 2015 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2496166
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, FEBRUARY 2015

Data
Intrinsic mean
Principal geodesic
Principal curve

14

4.5

9.6

Data
Intrinsic mean
Principal geodesic
Principal curve

Data
Intrinsic mean
Principal geodesic
Principal curve
8.6-5

Poincaré half-plane

22

Poincaré half-plane

6.5

Poincaré half-plane

10.6

6

0

5

2.5

-6

-1

4

6
-3

0

3

Fig. 3. Synthetic data on the Poincaré half-plane along with the principal geodesic and the principal curve. The
principal curve consistently better at capturing the trend of the data.
Data
Intrinsic mean
Principal geodesic
Principal curve

0

1.5

0.08

Data
OriginPrincipal geodesic
Principal curve

Error minimum

Test Error

1

0.04
0

1

-1
1

0

0
-1

-1

0
0

0.2
Smoother bandwidth (σ)

0.4
-1.5

0

1.5

Projection on
principal geodesic

Projection on
principal curve

Data

Fig. 4. Left: Motion capture data is represented on the power manifold SO(3)32 ; here the data is shown in the firsttwo (Euclidean) principal components in the tangent space at the intrinsic mean. The horizontal (yellow) line isthen the principal geodesic. The (purple) closed curve is the principal curve mapped to the tangent space andprojected into the first two principal components. The principal curve provide the best description of the data.
Center: The cross-validation error of the smoothing parameter in the experiment corresponding to the left panel.
Right: The relative orientation of the left femur of a person walking is represented as a point on the unit sphere.
Here, the periodicity of walking implies that the data roughly lies on a small circle; this imply an intrinsic mean at
the pole of the small circle. The principal geodesic, thus, fail to capture the trend of the data. The (closed) principal
curve does better.

Fig. 5. Motion capture data along with the corresponding projection onto both the principal geodesic and theprincipal curve. The principal curve captures the structure of data much better than the principal geodesic.
Copyright (c) 2015 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2496166
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, FEBRUARY 2015

Local PCA metric

0.6

Local PCA metric

0.6

0

0

0

Data
Intrinsic mean
Principal geodesic
Principal curve

Data
Intrinsic mean
Geodesics

Data
Intrinsic mean
Geodesics

-0.8
-1.5

Local PCA metric

0.6

0

2

-0.8
-1.5

7

0

2

-0.8
-1.5

0

2

Fig. 7. MNIST images of the digit 1 is projected onto its first two (Euclidean) principal components, and a smoothly
changing local diagonal metric is estimated as a local inverse covariance matrix. Left panel show geodesics
between data and the intrinsic mean; geodesics generally stay with the support of the data. Center panel show
geodesics originating at the intrinsic mean, shooting in different directions. Right panel show the principal geodesic
and the principal curve. Even when geodesics follow the trend of the data, the principal curve captures the overall
trend better than the principal geodesic.

Fig. 6. Example images from the MNIST dataset. The
orientation constitute the main source of variation.

5

S UMMARY AND D ISCUSSION

Principal geodesic analysis (PGA) is one of the most widelyused tools for statistical analysis of data on Riemannian
manifolds. This approach seeks a geodesic through the
intrinsic mean, which captures the main trend of the data.
This construction is a straight-forward generalization of
Euclidean PCA, which makes PGA both easy to understandand implement. The requirement that the estimated curve is
geodesic and passes through the mean is at the same timethe main strength and the main weakness of PGA. Thesestrong constraints on the solution provide much needed regularization when only small quantities of data are available.
However, the same constraints quickly prevent PGA from
capturing the trend; in particular the requirement that theestimated curve passes through the mean is detrimental.
Even when a suitable geodesic exists, it may not be foundas the mean may fall outside the support of the data.
We have proposed to extend the classic principal curves
from Hastie & Stuetzle [31] to Riemannian manifolds. The
resulting model is easy to implement as it only relieson the computation of locally weighted means. As the
local means are more likely to be unique than the global
counterpart (Theorem 1), the model is less sensitive tocurvature than PGA. The principal curve also forgoes theneed for a geodesic solution as the resulting curve is flexibly
represented as a piece-wise geodesic curve. This is both astrength and a limitation of the proposed approach. Fromsmall quantities of data it may be unrealistic to avoid
overfitting, but when more data is available, the increased
flexibility generally yields substantial improvements overthe inflexible geodesic solution. This is a standard dilemma.

An added benefit of the principal curve is that it cantrivially fit a closed curve to the data, even in spaces that do
not have closed geodesics. We have seen how this gracefully
allows for modeling periodic phenomena.
The main limitation of the principal curve framework is
the curve representation. By representing the curve as adiscrete set of points it is difficult to scale the approach toestimate higher dimensional surfaces. In principle, it is trivial
to extend the method to estimate a D-dimensional surface[31], but since the number of needed discretization points
grows exponentially with D, this quickly get impractical.
The running time of the method exceeds that of PGA as
we use this for initialization. The repeated calculation of
local intrinsic means may seem daunting from a computational perspective, but in practice this does not appear tobe an issue. As the means are only local, they tend to notinclude many data points, but more importantly, they tend
to be less sensitive to curvature. This implies that only few
iterations of the gradient descent scheme (Algorithm 1) areneeded to reach convergence. When numerical solvers areused to compute geodesics it is often more time-consuming
to compute geodesics between far-away points than nearby
ones. This also influence the computation of a local mean
as it only relies on nearby points.
Finally, we argue that it is beneficial to generalize nonlinear Euclidean statistics to Riemannian manifolds, rather
than linear statistics. These non-linear models often only
rely on local computations, which make them ideal for nonlinear representation spaces. We further expect the principal
curves to be extendable to more general metric spaces [42],[43] as they do not depend on tangent space constructions.

ACKNOWLEDGMENTS
Søren Hauberg is funded in part by the Danish Council for
Independent Research (DFF), Natural Sciences. The author
is grateful to Aasa Feragen and anonymous reviewers forhelpful feedback on early drafts of the manuscript.

Copyright (c) 2015 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.

This is the author's version of an article that has been published in this journal. Changes were made to this version by the publisher prior to publication.
The final version of record is available at http://dx.doi.org/10.1109/TPAMI.2015.2496166
SUBMITTED TO IEEE TRANSACTIONS ON PATTERN ANALYSIS AND MACHINE INTELLIGENCE, VOL. XXX, NO. XXX, FEBRUARY 2015

R EFERENCES
D. Kendall, “Shape manifolds, Procrustean metrics, and complex
projective spaces,” Bulletin of the London Mathematical Society,vol. 16, no. 2, pp. 81–121, 1984.
[2] S. Taheri, P. Turaga, and R. Chellappa, “Towards view-invariant
expression analysis using analytic shape manifolds,” in IEEE International Conference on Automatic Face & Gesture Recognition, 2011,
pp. 306–313.
[3] O. Freifeld and M. Black, “Lie bodies: A manifold representationof 3D human shape,” in European Conference on Computer Vision(ECCV), ser. Part I, LNCS 7572, A. Fitzgibbon et al. (Eds.), Ed.
Springer-Verlag, 2012, pp. 1–14.
[4] A. Srivastava, S. H. Joshi, W. Mio, and X. Liu, “Statistical shapeanalysis: Clustering, learning, and testing,” IEEE Transactions on
Pattern Analysis and Machine Intelligence (TPAMI), vol. 27, no. 4,
pp. 590–602, 2005.
[5] L. Younes, “Spaces and manifolds of shapes in computer vision: An
overview,” Image and Vision Computing, vol. 30, no. 6, pp. 389–397,
2012.
[6] S. Kurtek, E. Klassen, J. C. Gore, Z. Ding, and A. Srivastava, “Elastic
geodesic paths in shape space of parameterized surfaces,” IEEE
Transactions on Pattern Analysis and Machine Intelligence (TPAMI),
vol. 34, no. 9, pp. 1717–1730, 2012.
[7] C. Lenglet, R. Deriche, and O. Faugeras, “Inferring white matter
geometry from diffusion tensor MRI: Application to connectivity
mapping,” in European Conference on Computer Vision (ECCV), ser.
Lecture Notes in Computer Science, vol. 3024, 2004, pp. 127–140.
[8] X. Pennec, P. Fillard, and N. Ayache, “A Riemannian framework fortensor computing,” International Journal of Computer Vision (IJCV),
vol. 66, no. 1, pp. 41–66, 2006.
[9] Y. Wang, H. Salehian, G. Cheng, and B. C. Vemuri, “Tracking on
the product manifold of shape and orientation for tractography fromdiffusion MRI,” in IEEE Conference on Computer Vision and Pattern
Recognition (CVPR), 2014, pp. 3051–3056.
[10] O. Tuzel, F. Porikli, and P. Meer, “Region covariance: A fast descriptor
for detection and classification,” in European Conference on Computer
Vision (ECCV). Springer, 2006, pp. 589–600.
[11] F. Porikli, O. Tuzel, and P. Meer, “Covariance tracking using model
update based on Lie algebra,” in IEEE Conference on Computer
Vision and Pattern Recognition (CVPR), vol. 1, 2006, pp. 728–735.
[12] O. Freifeld, S. Hauberg, and M. J. Black, “Model transport: Towards
scalable transfer learning on manifolds,” in IEEE Conference on
Computer Vision and Pattern Recognition (CVPR), 2014.
[13] P. K. Turaga, A. Veeraraghavan, A. Srivastava, and R. Chellappa,“Statistical computations on Grassmann and Stiefel manifolds forimage and video-based recognition,” IEEE Transactions on PatternAnalysis and Machine Intelligence (TPAMI), vol. 33, no. 11, pp.
2273–2286, 2011.
[14] H. Çetingul and R. Vidal, “Intrinsic mean shift for clustering on
stiefel and grassmann manifolds,” in IEEE Conference on Computer
Vision and Pattern Recognition (CVPR), 2009, pp. 1896–1902.
[15] S. Said, N. Courty, N. Le Bihan, and S. J. Sangwine, “Exact principal
geodesic analysis for data on SO(3),” in Proceedings of the 15th
European Signal Processing Conference, 2007, pp. 1700–1705.
[16] S. Hauberg, S. Sommer, and K. S. Pedersen, “Natural metrics andleast-committed priors for articulated tracking,” Image and Vision
Computing, vol. 30, no. 6-7, pp. 453–461, 2012.
[17] ——, “Gaussian-like Spatial Priors for Articulated Tracking,” in
ECCV, ser. LNCS, K. Daniilidis, P. Maragos, and N. Paragios, Eds.,
vol. 6311. Springer, 2010, pp. 425–437.
[18] S. Hauberg, F. Lauze, and K. S. Pedersen, “Unscented Kalman
Filtering on Riemannian Manifolds,” Journal of Mathematical
Imaging and Vision, 2011.
[19] S. Huckemann, T. Hotz, and A. Munk, “Intrinsic shape analysis:
geodesic PCA for Riemannian manifolds modulo isometric Lie group
actions,” Statistica Sinica, vol. 20, no. 1, pp. 1–58, 2010.
[20] J. Lee, Introduction to Smooth Manifolds. Springer, 2002.
[21] M. Berger, A Panoramic View of Riemannian Geometry. Springer,
2007.
[22] X. Pennec, “Intrinsic statistics on Riemannian manifolds: Basic tools
for geometric measurements,” Journal of Mathematical Imaging andVision, vol. 25, no. 1, pp. 127–154, 2006.
[23] H. Karcher, “Riemannian center of mass and mollifier smoothing,”
Communications on pure and applied mathematics, vol. 30, no. 5,
pp. 509–541, 1977.
[1]8[24] W. S. Kendall, “Probability, convexity, and harmonic maps with small
image i: uniqueness and fine existence,” Proceedings of the London
Mathematical Society, vol. 3, no. 2, pp. 371–406, 1990.
[25] P. T. Fletcher, C. Lu, and S. Joshi, “Statistics of shape via principal
geodesic analysis on lie groups,” in IEEE Conference on Computer
Vision and Pattern Recognition (CVPR), vol. 1, 2003, pp. 1–7.
[26] S. Sommer, F. Lauze, and M. Nielsen, “Optimization over geodesics
for exact principal geodesic analysis,” Advances in Computational
Mathematics, vol. 40, no. 2, pp. 283–313, 2014.
[27] P. T. Fletcher, C. Lu, S. M. Pizer, and S. Joshi, “Principal Geodesic
Analysis for the study of Nonlinear Statistics of Shape,” IEEE
Transactions on Medical Imaging (TMI), vol. 23, no. 8, pp. 995–1005,
2004.
[28] V. M. Panaretos, T. Pham, and Z. Yao, “Principal flows,” Journal of
the American Statistical Association (JASA), vol. 109, no. 505, pp.
424–436, 2014.
[29] S. Jung, I. L. Dryden, and J. Marron, “Analysis of principal nested
spheres,” Biometrika, vol. 99, no. 3, pp. 551–568, 2012.
[30] M. Zhang and P. Fletcher, “Probabilistic principal geodesic analysis,”in Advances in Neural Information Processing Systems (NIPS), 2013,
pp. 1178–1186.
[31] T. Hastie and W. Stuetzle, “Principal curves,” Journal of the American
Statistical Association, vol. 84, no. 406, pp. 502–516, 1989.
[32] P. Diggle, “A kernel method for smoothing point process data,”
Journal of the Royal Statistical Society. Series C (Applied Statistics),
vol. 34, no. 2, pp. 138–147, 1985.
[33] R. Tibshirani, “Principal curves revisited,” Statistics and Computing,vol. 2, no. 4, pp. 183–190, 1992.
[34] T. Duchamp and W. Stuetzle, “Extremal properties of principal curvesin the plane,” The Annals of Statistics, vol. 24, no. 4, pp. 1511–1520,
08 1996.
[35] S. Gerber and R. Whitaker, “Regularization-free principal curve
estimation,” Journal of Machine Learning Research (JMLR), vol. 14,
pp. 1285–1302, 2013.
[36] S. Amari and H. Nagaoka, Methods of information geometry, ser.
Translations of mathematical monographs; v. 191. AMS, 2000.
[37] C. Ionescu, D. Papava, V. Olaru, and C. Sminchisescu, “Human3.6m:
Large scale datasets and predictive methods for 3d human sensing in
natural environments,” IEEE Transactions on Pattern Analysis andMachine Intelligence (TPAMI), vol. 36, no. 7, pp. 1325–1339, 2014.
[38] C. S. Catalin Ionescu, Fuxin Li, “Latent structured models for humanpose estimation,” in International Conference on Computer Vision(ICCV), 2011.
[39] H. Wang and T. C. M. Lee, “Automatic parameter selection fora k-segments algorithm for computing principal curves,” Pattern
Recognition Letters, vol. 27, no. 10, pp. 1142–1150, Jul. 2006.
[40] S. Hauberg, O. Freifeld, and M. J. Black, “A geometric take on
metric learning,” in Advances in Neural Information Processing
Systems (NIPS) 25, P. Bartlett, F. Pereira, C. Burges, L. Bottou, and
K. Weinberger, Eds. MIT Press, 2012, pp. 2033–2041.
[41] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner, “Gradient-basedlearning applied to document recognition,” Proceedings of the IEEE,
vol. 86, no. 11, pp. 2278–2324, 1998.
[42] A. Feragen, M. Owen, J. Petersen, M. Wille, L. Thomsen, A. Dirksen,
and M. de Bruijne, “Tree-space statistics and approximations forlarge-scale analysis of anatomical trees,” in Information Processingin Medical Imaging (IPMI), ser. Lecture Notes in Computer Science.
Springer, 2013, vol. 7917, pp. 74–85.
[43] S. Gerber, T. Tasdizen, P. T. Fletcher, S. Joshi, R. Whitaker, A. D. N.
Initiative et al., “Manifold modeling for brain population analysis,”
Medical image analysis, vol. 14, no. 5, pp. 643–653, 2010.

Søren Hauberg received his Ph.D. in computer science from the Univ. of Copenhagenin 2011. He has been a visiting scholar at
UC Berkeley (2010), and a post doc at thePerceiving Systems department at the Max
Planck Institute for Intelligent Systems (2012–2014). He is currently a post doc at the Section for Cognitive Systems, at the technical
Univ. of Denmark. He works in the span ofgeometry and statistics.

Copyright (c) 2015 IEEE. Personal use is permitted. For any other purposes, permission must be obtained from the IEEE by emailing pubs-permissions@ieee.org.