Noname manuscript No.
(will be inserted by the editor)

3D Human Pose Tracking Priors using Geodesic Mixture
ModelsEdgar Simo-Serra · Carme Torras · Francesc Moreno-Noguer

Received: date / Accepted: date

Abstract We present a novel approach for learninga finite mixture model on a Riemannian manifold inwhich Euclidean metrics are not applicable and oneneeds to resort to geodesic distances consistent withthe manifold geometry. For this purpose, we draw inspiration on a variant of the expectation-maximization
algorithm, that uses a minimum message length criterion to automatically estimate the optimal numberof components from multivariate data lying on an Euclidean space. In order to use this approach on Riemannian manifolds, we propose a formulation in which
each component is defined on a different tangent space,thus avoiding the problems associated with the loss of
accuracy produced when linearizing the manifold witha single tangent space. Our approach can be applied
to any type of manifold for which it is possible to estimate its tangent space. Additionally, we consider using
shrinkage covariance estimation to improve the robustness of the method, especially when dealing with very
sparsely distributed samples. We evaluate the approachon a number of situations, going from data clustering on
manifolds to combining pose and kinematics of articulated bodies for 3D human pose tracking. In all cases,
we demonstrate remarkable improvement compared toseveral chosen baselines.
Keywords Probabilistic Priors · Mixture Modelling ·Riemannian Manifolds · 3D Human Pose · Human
Kinematics
E. Simo-Serra (Corresponding author)
Waseda University
Okubo 3-4-1, Shinjuku, 169-8555 Tokyo, Japan
E-mail: esimo@aoni.waseda.jp
C. Torras · F. Moreno-Noguer
Institut de Robòtica i Informàtica Industrial (CSIC-UPC)
08028, Barcelona, Spain
E-mail: {torras,fmoreno}@iri.upc.edu

1 Introduction
The use of Riemannian manifolds and their statistics
has recently gained popularity in a wide range of applications involving non-linear data modeling. For instance, they have been used to model shape changes
in the brain [13], diffusion tensor imaging [44], deformations of anatomical parts [18] and human motion[6,60]. In this work we tackle the problem of approximating the Probability Density Function (PDF) of apotentially large dataset that lies on a known Riemannian manifold. We address this by creating a completely
data-driven algorithm consistent with the manifold, i.e.,an algorithm that yields a PDF defined on the manifold. This PDF can then be used as a prior in higherorder models, by combining it with image evidence inhybrid discriminative-generative models [53], or by exploiting it to constrain the search space in a tracking
framework [1]. We will show particular applications of
the proposed prior in the case of 3D human pose estimation, demonstrating a remarkable improvement compared to other widely used models.
A standard procedure to operate on a manifold is to
use the logarithmic map to project the data points onto
the tangent space of the mean point on the manifold [18,24,60]. After this linearization, Euclidean statistics arecomputed and projected back to the manifold using the
exponential map. This process is iteratively repeateduntil convergence of the computed statistics. Unfortunately, while this approximation is effective to model
data with a reduced extent, it is prone to fail when
dealing with data that covers wide regions of the manifold.
In the proposed finite mixture model, we overcome
this limitation by simultaneously considering multiple
tangent spaces, distributed along the whole manifold

2Fig. 1: Illustration of the proposed mixture
model approach. Each mixture component has its
own tangent space, ensuring the consistency of themodel while minimizing accuracy loss.

samples are sparsely distributed on the manifold (because the manifold has a very large dimensionality or
because the number of samples is small, or a combination of both effects). This makes our approach both
appropriate to handle situations with either large orsmall amounts of data, while our previous versions were
mostly effective when dealing with large datasets. We
finally unify and extend the evaluation in [55,56] to consider more manifolds and the improvements proposed inthis paper. Results will show that our manifold-based
finite mixture model can be used to exploit the known
structure of the data, outperforming approaches that
do not. We provide the source code1 of our approach.

2 Related Workas seen in Fig. 1. We draw inspiration on the unsupervised algorithm from [17], which given data lying in an
Euclidean space, automatically computes the numberof model components that minimizes a message lengthcost. By representing each component as a distribution on the tangent space at its corresponding meanon the manifold, we are then able to generalize thealgorithm to Riemannian manifolds and at the sametime mitigate the accuracy loss produced when using
a single tangent space. Furthermore, since our model
is semi-parametric, we can handle an arbitrarily large
number of samples. This is in contrast to existing nonparametric approaches [41] whose complexity grows withthe training set size.
As an example of practical application of our mixture model, we will consider the 3D human pose tracking problem, which has been traditionally addressedwith kinematic priors based on Gaussian diffusion [16,20, 23,52,58]. This consists in simply searching in a small
area defined by a Gaussian distribution centered on theprevious pose, i.e., xt = xt−1 +, where xt would be thepose at time t and  would be a Gaussian perturbation
with 0 mean and diagonal covariance. However, thissimple model does not constrain the pose to lie on its
underlying manifold, and does indeed explore a muchhigher dimensional space than it should be strictly necessary. We will show that using our model as a kinematic prior we can effectively focus our solution on theactual manifold, greatly outperforming standard Gaussian diffusion models.
A preliminary version of this work appeared in [55]with an application to introducing kinematic priors later
presented in [56]. We extend this work by considering
improvements for the covariance estimation. In particular, we consider shrinkage estimators that are shownto outperform empirical covariance estimation when the

Manifolds have always been very important in computer vision [47]. The two more widely used manifolds
have been the one of symmetric semi-definite matrices [21,22,28,29,43,57], and the Grassman manifolds [27,50,67]. However, most of these approaches focus on exploiting very specific manifolds and do not generalize toother manifolds. In contrast, our approach is applicable
to all Riemannian manifolds with explicit exponential
and logarithmic maps.
Recently, there has been an influx of theoretical results in statistics on Riemannian manifolds [42] thathave allowed for their widespread usage in modeling.
For example, there exists several PCA generalizations
to non-linear data such as the Principal Geodesic Analysis [18,60] and the Geodesic PCA [24]. Yet, these methods only use one single tangent space located at thegeodesic mean, which can lead them to have significant
accuracy error when input data is spread out widely on
the manifold. Other algorithms based on kernels [13]and Markov Chain Monte Carlo sampling [6] have been
specifically used in regression tasks on manifolds, butthey have not been applied to stochastic modeling problems. There have been recent attempts at removing thetangent space linearization [59,75], which, however, can
not yet scale to the large amounts of data we consider
in this work.
Other approaches address classification models on
Riemannian manifolds [47,64,65,68]. For binary cases,
the classifier is usually built in a “flattened” version of
the manifold, obtained via the tangent space [68]. Multiple category classification problems have been tackled
by replacing the tangent space mapping with rolling
maps [9], and by using extensions of the Kernel methods
to Riemannian manifolds [28,29]. In any event, theseapproaches have been exclusively used for classification1

http://hi.cs.waseda.ac.jp/~esimo/code/gfmm/

3

Model

Complexity

Scales?

Consistent?

PDF?

Gaussian diff.
GMM [51]
PGA [18]
GPLVM [31]GPDM [73]hGPLVM [32]CRBM [63]
GCMFA [36]

Low
Low
Low
Low
MediumMedium
High
High

Yes
Yes
Yes
No
No
No
Yes
No

No
No
Yes
No
No
No
No
No

Yes
Yes
No
Yes
Yes
Yes
Yes
Yes

GFMM (Ours)

Low

Yes

Yes

Yes

Table 1: Comparison of several commonly usedhuman pose models. Models are considered to scale
if they can handle well large amounts of data (∼100K
samples) and to be consistent if they use geodesic distances instead of other metrics. The last column reflects
whether or not a model is actually modeling the Probability Density Function (PDF) of the data.

problems, which are out of the scope of the current paper, focused on PDF modeling for use as priors.
With regards to density estimation on Riemannian
manifolds, various non-parametric approaches [40,41]
have been proven to be appropriate. However, as theircomplexity is dependent on the number of training samples, they scale poorly for large datasets. In contrast,
semi-parametric models such as the mixture model we
propose here can handle large amounts of data efficiently. Dirichlet Processes have been used for fitting
mixture models [10], and recently modified to handle
the case of the sphere manifold [62], although, compared to our approach, they have not been extended toarbitrary Riemannian manifolds. Another widely studied manifold is that of tensor fields [35], for which a nonparametric Kernel Density Estimation approach was recently proposed [8]. In [39], individuals on the tangent
bundle are modeled and populations are compared withgeneralized statistical hypothesis tests, but no parametric model is learned. The interesting approach in [2]is similar to ours in spirit, as it proposes a techniqueto perform Expectation Maximization (EM) on manifolds. However, this work considers data-driven manifolds, resulting in a high computational overhead forlarge training sets. In addition, it neither estimates thenumber of clusters, nor makes use of the tangent space,which allows our model to be defined “on the manifold”.
As for human pose, it has been traditionally modeled as a tree of connected joints [26,37,38]. There havebeen many different ways of modeling this. One of themost straightforward approaches is to make use of a
Gaussian Mixture Model (GMM) [51]. Another popular trend is to use Gaussian Processes (GP)-based ap-

proaches, such as GP-Latent Variable Models [31] andthe constrained GP [71]. These have been extended toconsider dynamics in the Gaussian Process Dynamic
Model (GPDM) [69,73,74], and also to consider topological constraints [70]. Hierarchical variants [32] (hGPLVM) have also been used in tracking-by-detection [1].
However, Gaussian Processes do not scale well to large
datasets due to their O(n3 ) complexity for prediction.
Sparse approximations do exist [45], but in general do
not perform as well. In contrast, once our model has
been estimated it has O(1) complexity for sampling.
There have been other approaches for modeling human pose such as learning Conditional Restricted Boltzmann Machines (CRBM) [63]. However, these methods
hold on a complex learning procedure that uses several
approximations, and make the training of good models
harder. Li et al. [36] proposed the Globally CoordinatedMixture of Factor Analyzers (GCMFA) model which issimilar to the GPLVM ones in the sense it is performing a strong non-linear dimensionality reduction. Yet,as GPLVM, it does not scale well to large datasets suchas the ones we consider in this work.
We would like to point out that none of the aforementioned approaches is consistent with the manifoldof human motion. Some of them use directly the 3D position of the joints while others use angles. In the case of
considering 3D points the limb length may vary during
the tracking, which is neither realistic nor desirable. In
the case of angle representations, they have an inherent
periodicity and thus are not a vector space even though
they are usually treated as such. Two nearby angles may
have very different values, e.g., 0 and 2π. In this case thedistance using the angular value would be 2π while thetrue geodesic distance is 0. Our model can handle boththese limitations. Another model that can handle thisis the Principal Geodesic Analysis (PGA) [18]. However, this model uses a single tangent space and does
not model the Probability Density Function (PDF).
Table 1 summarizes the properties of the models wehave commented in terms of their complexity, ability
to scale, manifold consistence and if they provide or
not a PDF. In particular, our approach scales well, isconsistent with the manifold, has low complexity (i.e.,
it just considers a single hyperparameter), and can beeasily learned using an Expectation-Maximization algorithm. It is worth noting here that our model is also
the fastest of them for sampling (it is O(1)). For instance, our Matlab implementation yields over 100,000
samples per second. An example of our model can beseen in Fig. 2.

4Ground Truth

ZoomedSampling

3
2
1
0
−1
−2
−3

0.2

Right Hand Tangent Space

0.08

25

0.15

0.06
20

0.1

0.04

0.05

15

0.02

10

0
0

5

−0.05
−0.02

−0.1

0

−0.04

−5

−0.2

−0.06

−10

−0.25
−0.4

−0.08
−0.15

−0.15

−0.3

−0.2

−0.1

0

0.1

0.2

0.3

−3 −2 −1

0

1

2

3

Right Elbow Tangent Space

−0.1

−0.05

0

0.05

0.1

Fig. 2: Example of our motion prior. Top: For aparticular pose, 100 motion samples of the predicteddistribution are shown. For visualization purposes themagnitude of the samples is multiplied by 3. Bottom:
Visualization of some of the joint samples with theirassociated log-likelihood. The ground truth is shownwith a black diamond.

Fig. 3: Representation of geodesics on the
S 2 manifold. The tangent space ensures that
→ However,k logp (x)k is the true geodesic distance of −px.
−→
k logp (a) − logp (b)k is not the geodesic distance of ab.
a smooth inner product on the tangent space Tp M ateach point p on the manifold. Consider a parametrized∂
curve γ : [0, 1] → M with velocity γ̇(t) = ∂tγ(t). Ageodesic is a curve that minimizes the distance between
the two points p = γ(0) and x = γ(1). More formally, ageodesic is a curve with null acceleration along M , i.e.,Dthe covariant derivative ∂tγ̇(t) is 0 for all t ∈ [0, 1].
We will denote the length of the geodesic or geodesic
distance asZ 1qd(p, x) =gγ (t) (γ̇(t), γ̇(t)) dt .
(1)
0

3 Geodesic Finite Mixture Model
We next describe our approach, starting with some basic notions on Riemannian geometry and statistics on
manifolds. We then integrate these ingredients in a mixture modeling algorithm to build manifold-based generative models.

3.1 Manifolds, Geodesics and Tangent Spaces
Manifolds arise naturally in many real-world problems.
One of the most well-known is the manifold representing spatial rotations. For example, when studying human motion, it is a common practice to use the spatial
rotations of the different body parts to obtain a subjectagnostic representation of the whole body pose. Thisis usually done with angle representations that havean inherent periodicity and thus are not a vector space.
By considering the Riemannian manifold of spatial rotations it is possible to use tangent spaces as a local
vector space representation, and use powerful statistical tools based on Euclidean metrics. For an in depth
description of Riemannian manifolds we refer the reader
to [4,7].
A Riemannian manifold (M, g) is a differentiable
manifold M equipped with a metric g, that provides

We can now define the exponential map expp at p =γ(0) and its inverse, the logarithmic map logp as
expp :

Tp M −→ M
v 7−→ expp (v) = γ(1) = x

logp :

M −→ Tp M
x 7−→ logp (x) = v

.

,
(2)

The exponential map is locally diffeomorphic ontoa neighborhood of p. Let V (p) be the largest such neighborhood, then logp (x) is defined for any point x ∈ V (p).
Geodesics γ(x,v) (t) = expx (tv) from t = 0 to infinity
can either be minimizing all the way or only up to atime t0 < ∞ and not any further. In this latter case,the point z = γ(x,v) (t0 ) is called a cut point. The setof all cut points forms the cut locus, and the corresponding vectors the tangential cut locus. The maximal
domain of V (p) will be the domain containing 0 anddelimited by the tangential cut locus. The geodesic distance can also be written using the logarithmic map asd(p, x) = k logp (x)k (see Fig. 3).
In general there is no closed-form of the expp andlogp maps for an arbitrary manifold. There are, though,approximations for computing them in Riemannian manifolds [14,61]. Additionally, efficient closed-form solutions exist for certain manifolds [46]. We will discuss5

some of these manifolds and their associated expp andlogp maps in the next section.

approximate normal distributions on the tangent space
Tµ M at the mean µ with covariance matrix:

3.2 Statistics on Tangent Spaces

1 X
Σ=logµ (xi ) logµ (xi )> ,
N i=1

N

While it is possible to define distributions on manifolds [42], we will focus on approximating Gaussian
PDFs of data on a manifold using the tangent space.
For instance, the mean of N points xi on a manifold
can be calculated as [30]:µ = arg min
p

N
X

2

d (xi , p) .

(3)

i=1

This is iteratively optimized using the expp and logp
maps,
!
N
δ X
µ(t + 1) = expµ(t)logµ(t) (xi ) ,
(4)
N i=1
until kµ(t + 1) − µ(t)k <  for some threshold , with δ
being the step size parameter.
Knowing the mean value µ and the concentration
matrix Γ we can write the distribution that maximizes
entropy on the tangent space as a normal distribution
centered on the point µ ∈ M, corresponding to theorigin (ν = 0) in the tangent space:!
logµ (x)> Γ logµ (x)Nµ (ν, Γ ) = a exp −,(5)2where the normalization constant a and covariance matrix Σ are related to the concentration matrix by:
!
Z
>
log
(x)
Γ
log
(x)
µ
µ
a−1 =
exp −
dM(x) ,
(6)
2
M
and
Z
Σ=a

logµ (x)> logµ (x)Mlogµ (x)> Γ logµ (x)
exp −
2

!
dM(x) .

(7)

Note that this integral is over the manifold M and that
not all points of the tangent space Tµ M correspond to
one single point on the manifold, i.e., the tangential cut
locus. In particular, for the S 2 sphere, the tangent spaceis defined inside a circle and not the whole R2 plane.
The circle and the area outside of it forms the tangential
cut locus, and for any point on the tangential cut locus
there exists more than one minimizing geodesic to theorigin. As a simplification, in our formulation we will
not consider the tangential cut locus, and will directly

(8)

and
a−1 =

q

(2π)D det(Σ) ,(9)

By not taking into account the tangential cut locus weare underestimating the true normalization parameter.
We next perform a simple analysis of the error incurred by obviating the tangential cut locus for the S 2
sphere, which is at a distance π from the origin. We cancompute the exact normalization term by calculating
the integral of the 2D Gaussian on the tangent space
using Eq. (6). We consider the scenario of a normal
distribution centered at the origin and with diagonalcovariance matrix Σ = diag(σ, σ). Using polar coordinates we can write: 2
Z 2π Z π−r∗ −1(a ) =exp
r dr dθ
2σ
0
0



−π 2
= 2πσ 1 − exp
.
(10)
2σ
If we do not take into account the tangential cut locus,
the normalization term of Eq. (9) becomes a = 2πσ. As
expected, we are underestimating the true constant by afactor of a∗ /a = 1−exp(−π 2 /(2σ)). To make an estima−π 2tion error over a 1%, σ ≥ 2 log(0.01)≈ 1.072. Therefore,
unless the distribution has an extremely large covariance, for the S 2 manifold, the estimation error will beless than 1%. In the case of human articulations modeled with S 2 joints, most of them do not even have amovement range of 1.072 radians, and their covarianceis much smaller, making this error negligible. Furthermore, since we model the data as a mixture of Gaussians, each of the components of the mixture will havea much smaller covariance than the total covariance ofthe data. Experimentally, we found no difference between considering or not the tangential cut locus.
For an alternate approach to estimate a normal distribution on a Riemannian manifold where the Taylor
expansion of the Riemannian metric is used, please refer to [42].

3.3 Improving Covariance Estimations
We next discuss alternative approximations to estimate
covariance matrices. These new estimates will later beused in the place of the empirical covariance matrix Σ.

6The standard approach to compute a covariance matrix S = [sjk ], defines its entries as:
N

sjk =

1 X
(xij − µj )(xik − µk )> ,
N i=1

(11)

with µj and µk being the j-th and k-th element of themean of the N samples x with dimensionality D. Notethat this is what was used in Eq. (8).
This empirical covariance matrix is known to bea poor estimation of the true covariance matrix when
the number of samples is small compared to their dimensionality, yielding samples which are very sparsely
distributed. Several approaches have been proposed for
improving this approximation [34,49]. In this paper, we
will focus on Ledoit-Wolf (LW) [33] and Oracle Approximating Shrinkage (OAS) [11] techniques, as besides being accurate for small datasets, they are also efficient
to compute for large training sets. They both belong to
the so-called family of “shrinkage estimators”.
In linear shrinkage problems, the estimation of thecovariance matrix is formulated as a constrained MSE
minimization w.r.t. a true covariance Σ:Another interesting case is when input samples are
subject to noise, e.g. due to measurement uncertainty.
In this situation we can introduce a prior on the structure of the input data and, for instance, parameterize
every sample xi by a specific Gaussian distribution withcovariance Σxi :xi = yi + N (0, Σx−1),iwhere yi is the mean value of xi .
For the sake of completion we will consider eachsample xi to be weighted by wi (we will use this in thefollowing subsection, for the EM computation). WithPN
out loss of generality we assumei=1 wi = N . The
mean of the N samples can then be written as:
N

1 Xµ = E[x] = E [E[x|y]] = E[wy] =wi yi .
N i=1

(16)

By using the law of total covariance we can then write
the biased weighted sample covariance asΣ = cov(x) = E[cov(x|y)] + cov(E[x|y])= E[cov wN (x, Σx−1 ) ] + cov(wy)= E[w2 Σx ] + cov(wy)min E [kΣ ∗ − ΣkF /D]
ρ

N

tr(S)
I + (1 − ρ)S ,(12)Dp
where kAkF = tr(AA> ) is the Frobenius norm, ρ isthe shrinking coefficient, and I is the identity matrix.
Ledoit-Wolf [33] proposed using the following shrinking coefficient:!
PN>
i=1 kxi xi − SkF /DρLW = min,1 ,(13)
N 2 (tr(S 2 ) − tr2 (S)/D)
s.t. Σ ∗ = ρwhich is proven to converge to the optimal solution
when N, D → ∞ and D/N → c, 0 < c < ∞.
On the other hand, Chen et al. [11] propose to use
the Oracle Approximating Shrinkage (OAS):
ρOAS = min(15)!
(1 − 2/D)tr(S 2 ) + tr2 (S), 1 ,(N + 1 − 2/D) tr(S 2 ) − tr2 (S)/D(14)

which is the limiting form of the optimal oracle estimator or ideal value of ρ.
Both the LW and OAS shrinkage estimators havethe desirable property that the estimated covariance Σ
is in general invertible, unlike the empiric covarianceestimation. In the experimental section we will showthat they also provide a better covariance estimation
for a wide variety of problems.

=


1 X
wi Σxi wi + (yi − µ)(yi − µ)> .
N i=1

(17)

In the results section we will show how it is possible totreat the sample noise Σxi as a hyperparameter whenestimating mixtures from few samples, and that it behaves as a regularization parameter that helps to improve performance.

3.4 Unsupervised Finite Mixture Modeling
Recall that our ultimate goal is to fit a mixture modelon Riemannian manifolds. For this, we will draw inspiration on [17], a variant of the EM algorithm [15] thatuses the Minimum Message Length criterion (MML) [72]to estimate the number of clusters and their parametersin an unsupervised manner.
Given an input dataset, this algorithm starts by randomly initializing a large number of mixtures. During
the Maximization (M) step, a MML criterion is used toannihilate components that are not well supported by
the data. In addition, upon EM convergence, the least
probable mixture component is also forcibly annihilated
and the algorithm continues until a minimum numberof components is reached. The approach in [17] is designed to work with data in an Euclidean space. To use
it in Riemannian manifolds, we modify the M-step asfollows.

7

We define each mixture component with a mean µk
and a concentration matrix Γk = Σk−1 as a normal distribution on its own tangent space Tµk M:p(x|θk ) ≈ Nµk 0, Σk−1 ,(18)with θk = (µk , Σk−1 ). Remember that the mean µk
is defined on the manifold M, while the concentration matrix Γk is defined on the tangent space Tµk M
at the mean νk = 0. Also note that the dimensionality of the space embedding the manifold is larger
than the actual dimension of the manifold, which inturn is equal to the dimension of the tangent space.
That is, for an arbitrary embedding of the manifold,dim(Embedding(M)) > dim(Tp M) = dim(M) = D.
This dimensionality determines the total number of parameters Dθ specifying each component, and, as we will
explain below, plays an important role during the component annihilation process. For full covariance matrices it can be easily found that Dθ = D + D(D + 1)/2.
We next describe how the EM algorithm is extended
from Euclidean to Riemmanian manifolds. For a fullderivation of the algorithm please see Appendix A. Specifically, let us assume that K components survived after
iteration t − 1. Then, in the E-step we compute the responsibility that each component k takes for every sample xi :αk (t − 1)p(xi |θk (t − 1))(i)wk = PK,
k=1 αk (t − 1)p(xi |θk (t − 1))

Fig. 4: Representation of a conditioned distribution. In the left we show a joint Gaussian distribution
over two variables xA and xB . We then on the rightillustrate the resulting distribution p(xA |xB ) for a particular point xB = b. We can see that this is also a
Gaussian distribution, albeit one-dimensional.

on the Minimum Message Length is then computed.
This EM process is repeated until the convergence ofthe score or until reaching a minimum number of components Kmin . If this number is not reached, the component with the least responsibility is eliminated (evenif it is larger than Dθ /2) and the EM process is repeated. Finally, the configuration with minimum scoreis retained (see [17] for details), yielding a resulting distribution with the form

(19)

for k = 1, . . . , K and i = 1, . . . , N , and where αk (t − 1)
are the relative weights of each component k.
In the M-step we update the weight αk , the mean
µk and covariance Σk for each of the components asfollows:p(x|θ) =KX

αk p(x|θk ) .

(21)

k=1

N

1 X (i)wk
αk (t) =
w =
,
N i k
N

2
N
X
N (i) (i)µk (t) = arg mind
wk x , p,wkp
i=1Σk (t) =3.5 Conditional Distribution

N

>
1 X
(i)
logµk (t) (x(i) ) logµk (t) (x(i) ) wk ,wk i=1(20)
If we wish to augment the data with noise covarianceassociated to each sample, it is as simple as adding theweighted average of the noise covariance to Σk (t) as per
Eq. (17).
After each M-step, we follow the same annihilation
criterion as in [17], and eliminate those components
whose accumulated responsibility wk is below a Dθ /2
threshold. A score for the remaining components based

It is possible to use the generative model just describedin prediction tasks by estimating the distribution of asubset of variables given another subset of variables.
Essentially, given a joint distribution, the distributionof a subset of variables conditioned on another subset
of variables is computed as illustrated in Fig. 4. To dothis we need to split the dimensions of the manifold into
two subsets, x = (xA , xB ). Each of these components
can be expressed in terms of the mixture as:θk = (µk , Σk−1 ) =Σk,A Σk,AB(µk,A , µk,B ),Σk,BA Σk,B

−1 !
.
(22)

8

The conditional distribution of one subset given theother (i.e., the regression function), can be written as:p(xA , xB |θ)p(xB |θB )PK
αk p(xB |θk,B )p(xA |xB , θk,A )= k=1 PK.
k=1 αk p(xB |θk,B )p(xA |xB , θ) =(23)
Observe that this is a new mixture model:p(xA |xB , θ) =KX

πk p(xA |xB , θk ) ,

(24)

k=1

with weights:
αk p(xB |θk,B )πk = PK,
j=1 αj p(xB |θj,B )(25)

In the case of the Euclidean space, p(xA |xB , θk,A )can be written as:−1
p(xA |xB , θk,A ) = N (µk,A|B , Σk,A|B),−1µk,A|B = µk,A + Σk,AB Σk,B
(xB − µk,B ) ,−1
Σk,A|B = Σk,A − Σk,AB Σk,B
Σk,BA .

(26)

In our case we assume that the tangent spaces are notbeing moved and that they are centered on the mean.
This allows us to write the conditional probabilities as:−1
p(xA |xB , θk,A ) = Nµk,A|B (νk,A|B , Σk,A|B),−1
νk,A|B = Σk,AB Σk,B
logµk,B (xB ) ,−1
Σk,A|B = Σk,A − Σk,AB Σk,B
Σk,BA ,

(27)

where logµk,B is the subspace of the tangent space logµk
at the subset of the mean µk,B for the mixture k. Notethat the tangent spaces are not being moved, i.e., they
still remain centered on µ, although we are only looking
at a subspace of the tangent space.
In the results section we will show how we can use
this to predict the kinematics of a person given her/his
pose.

3.6 Implementation Considerations
While using the tangent spaces allows representing PDFs
of data on manifolds, this comes at a price of higher
computational cost as the data must be repeatedly projected back and forth from the tangent space to the
manifold. There are, though, several implementationconsiderations that can be taken into account to improve the efficiency.
For instance, as mentioned in [17], we might consider
using less expressive covariance matrices (e.g. diagonalones). However, when using tangent spaces, there is notnecessarily a global coordinate frame representation, as
the orthonormal basis of the tangent space depends on
the logp map, and thus, depends on the point p at which
it is calculated. When running the EM algorithm, thetangent spaces at step t + 1 may have a completely different basis than those at step t, and thus, the data
likelihood can change drastically, making the optimization much more non-linear. Since (in contrast to fullcovariance matrices) diagonal matrices can not adapt
to arbitrary rotations of the coordinate system, theirperformance is in general quite poor. This limitation
can only be bypassed when there exist a global coordinate frame that can be defined for all tangent spaces
independent of the point they are centered on, e.g., theEuclidean spaces Rn .
Nevertheless, when working with a manifold whichis the Cartesian product of other manifolds such as
SAB = SA × SB , it is possible to use a block-diagonal
matrix of the form:  −1

ΓSA 0
ΣSA 0
ΣSAB === ΓS−1
,
(28)
AB
0 ΣSB
0 ΓS−1
B
which by construction is a valid covariance matrix andavoids the issue of dealing with arbitrary orthonormal
basis. The null row and column elements highly simplify
the computational cost. By using less expressive covariance matrices, the model has fewer degrees of freedom
and generally converges in fewer iterations, besides requiring less training data.
Furthermore, while in the Euclidean case the meanof a set of points can be computed in closed form, whenworking with manifolds we need to do this iteratively.
In our implementation this is required in the M-step,
where the parameters θk are estimated as a weighted
combination of terms.
only a subset of
P By considering(i)samples S such that i∈S wk > s for a certain threshold s , it is possible to improve the computational efficiency without sacrificing accuracy.
In order to improve the initialization, we consider
using the k-means algorithm in the embedding space
as a coarse initialization for the algorithm. This ensures a fair spread of the initial clusters over the datain contrast to a random sampling that may fail when
the data is unevenly distributed. In the results section
we will show that manifolds with large dimension leadto fewer clusters being annihilated in the initial iteration and more stable convergence properties.
4 Manifolds
We will now describe several manifolds that we will use
in the experimental section. For each manifold we will

9

briefly discuss their structure and the expp and logp
map implementations.

previous one. The whole pose will thus be represented
as the Cartesian product of all the relative rotations between consecutive joints. We will write the human posemanifold H as

4.1 Quadratic surfaces
H = S2 × S2 × · · · × S2 .
In general, there is no closed-form of the expp and logp
maps for an arbitrary Riemannian manifold. There are,though, approximations for computing them [14,61]. In
particular we will consider implicitly defined surfaces.
Computing the expp map corresponds to solving an
initial value ordinary differential equation problem. Refer to [14] for a neat numerical algorithm to obtain solutions to the map.
The computation of the logp map is harder and isbased on a shooting algorithm. This relies on the expp
map to iteratively refine an initial guess, and is a natural generalization of error correction from the Euclidean
space to manifolds. We use the implementation proposed in [61].
4.2 S 2 Manifold
There is an explicit mapping between the unit sphere S 2and its tangent space Tp S 2 [46]. Let x = (x1 , x2 , x3 )> ,y = (y1 , y2 , y3 )> be two unit spoke directions in S 2 andv = (v1 , v2 )> a point in Tp S 2 . The expp and logp maps
are in this case:
sin kvk
sin kvk−1
expp (v) = Rp
v1
, v2
, cos kvk ,
kvk
kvk


θ
θ
logp (x) = y1
, y2
,
(29)
sin θ
sin θ
where Rp is the rotation of p to the north pole, kvk=
1
(v12 + v22 ) 2 , y = Rp x and θ = arccos(y3 ).

(30)

In this case the expp and the logp maps for the humanpose manifold will consist of the Cartesian product of
the expp and logp maps for all consecutive joints, which
is one less than the total number of joints. Note thatwhile there exist other manifolds for 3D human pose,such as one defined by using forward kinematics [23],they do not have closed form solutions for the expp andlogp operators.
4.4 Joint Pose and Kinematic Manifold
We can obtain a joint human pose and kinematics manifold with the tangent bundle of the human pose manifold, which we equip with a Riemannian metric called
the Sasaki metric [48]. The tangent bundle is definedas:
T M = {(x, v) | x ∈ M, v ∈ Tx M} .

(31)

Let (u, w) be a vector tangent to T M at a point(x, v). Both u and w must be lifted from the tangent
space Tx M to T(x,v) T M. The lift of the u component
is called the horizontal lift of u and denoted uh . Likewise, the lift of the w component is called the verticallift of w and denoted wv . Geodesics along uh move x
while parallelly translating u, whereas geodesics along
wv move v linearly while keeping x fixed. Given twoelements a = (x1 , v1 ), b = (x2 , v2 ) ∈ T M, the Sasaki
metric ĝ(a, b) is given asĝ(xh1 , xh2 ) = g(x1 , x2 ) ,4.3 Human Pose Manifold
The human pose is commonly represented using a discrete set of points in 3D space that correspond to different articulations [25,53,54]. For a specific individual,
the distance between two consecutive joints, e.g., elbow
and hand, is fixed. It is therefore common to separate
the specific characteristics of the individual given by the
distances between two consecutive joints, from his/herpose, i.e., the relative rotation between two consecutive
joints [26]. We will therefore consider the human poseas the set of relative rotations between all pairs of consecutive joints which forms a tree structure [60,66]. We
will further represent this relative motion as points ona sphere. That is, given a specific joint, the next consecutive joint will lay on the S 2 sphere centered on theĝ(xh1 , v2h ) = ĝ(xh2 , v1h ) = 0 ,ĝ(v1h , v2h ) = g(v1 , v2 ) ,(32)

where g is the metric on M. Notice that there are no
cross-terms between x and v.
For two consecutive poses x1 and x2 acquired at aconstant framerate, we can compute the velocity v12 between x2 and x1 directly on the tangent space throughthe logarithmic map at x1 and thus define the joint poseand kinematic manifold as:
T H = (x1 , v12 ) = x1 , logx1 (x2 ) | x1 , x2 ∈ M ,(33)
where kv12 k is the geodesic distance between both points
as shown in Fig. 5.

103

6 Clusters
Clusters Correct2

Ours 6.09 (0.32)
1-TM 7.05 (1.38)
vMF 16.59 (1.71)
DP-GMM [10] 7.83 (1.27)
DP-TGMM [62] 4.49 (1.98)

1
0
−1
−2

0.92
0.46
0.00
0.08
0.55

8 Clusters
Clusters Correct
8.00 (0.00)
15.25 (2.17)
19.86 (2.35)
3.43 (3.26)
8.34 (1.02)

1.00
0.00
0.00
0.26
0.51

−3
−3 −2

−1

0

1

2

3

Fig. 5: Visualization of the velocity. Velocities correspond to points on the tangent space at x1 . Givena point x2 (corresponding to the pose acquired rightafter x1 ), the velocity is the curve going from x1 tox2 which is equivalent to a straight line in the tangent
space. The Euclidean norm of v12 on the tangent spacecorresponds to the geodesic distance from x1 to x2 .
We next define the exponential and logarithmic maps
for the elements (x, v) ∈ T H. For a local neighbourhood of x, the elements v ∈ Tx H form a vector space
and thus the operators can be simply defined as in the
Euclidean case withlogv1 (v2 ) = v2 − v1

andexpv1 (v2 ) = v2 + v1 .

(34)

This ensures that the mean of the data on the tangent
space at the geodesic mean will be 0, i.e, the mean of{logp (xi )}xi will be 0 if p is the geodesic mean.
We can also further extend this manifold by including more poses acquired sequentially at a constantframe rate. For example, given three consecutive poses
x1 , x2 and x3 it is possible to map these points to amanifold in which one pose is the reference and theother two poses are considered offsets from that reference or tangent vectors. That is (x1 , x2 , x3 ) is mappedto (logx2 (x1 ), x2 , logx2 (x3 )) where x2 would be the local reference.

5 Results

Table 2: Recovering distributions on the Spheremanifold. We show results on recovering an original distribution on the Sphere manifold with 6 and
8 clusters. Results obtained from 100 different evaluations and 1000 samples per cluster. We compare ourmethod with using a single tangent space (1-TM), vonMises-Fisher distributions (vMF), and two variants of
Dirichlet Process Mixture Models (DP-GMM and DPTGMM).
4. Tracking prior by extending the manifold with its
tangent bundle.

5.1 Recovering Distributions on Quadratic Surfaces
We first present experiments on two 2D manifolds defined by
M = (x, y, z) | cx2 + y 2 + z 2 = 1 .
(35)
For the first example we generate 1800 points on thesphere S 2 (i.e. c = 1) using 6 clusters with parametersµi = [cos(iπ/3), 0, sin(iπ)/3] and Σi = diag(0.2, 0.3)
for i = 1, . . . , 6. The algorithm is initialized with 30
clusters. This manifold has the closed-form solution for
the expp and logp operators given by Eq. (29), allowing
the method to execute in under a minute. Fig. 6 shows
how the final solution retrieves the 6 clusters used togenerate the data.
For the second example we generate 1500 points
from a mixture of 5 Gaussians on the manifold of Eq. (35)with c = −2, as shown in Fig. 7a. The Gaussian parameters used in this case are:In this section we will extensively evaluate our model
on both synthetic and real data. We consider a number of manifolds and evaluation procedures to show theflexibility and diversity of our model. The section is
structured as follows:µ1 = [0.83, −1.09, 1.09]µ2 = [0, 0, 1]µ3 = [0.30, −0.77, 0.77]µ4 = [0, 0, 1]µ5 = [0, −0.89, 0.44]

1. Evaluation on the recovery of synthetic distributions
on quadratic surfaces.
2. In depth results on modeling synthetic distributions
on the S 2 sphere.
3. Modeling human pose and kinematics on the largescale Human3.6M dataset [26].

It is worth to highlight two important details about thismixture. First, the covariances depend on the local orthonormal basis of the tangent plane, thus even if they
are diagonal, in practice they are not when projectedback onto the manifold, as shown in Fig. 7a. Second,clusters 2 and 4 share the same mean. In this example,

Σ1 = diag(0.20, 0.30)
Σ2 = diag(0.25, 0.10)
Σ3 = diag(0.20, 0.10)
Σ4 = diag(0.10, 0.25)
Σ5 = diag(0.25, 0.10)

11

30

4200

Clusters

Cost

4000
3800
3600

20

10

3400
3200

0

100

(a) Input Data

200

300

400

100

200

300

400

Iteration

Iteration(b) Cost Function(c) Number of Clusters(d) Solution Cluster #1

Fig. 6: Sphere Example. (a) Input data on the manifold. It consists of 1800 points sampled from a mixture of
6 Gaussians. Note that they are colored only for visualization purposes; our algorithm does not know a priorithese clusters. (b) Evolution of the cost function, based on the Minimum Message Length of all componentsof the mixture. Vertical lines represent iterations in which a cluster is annihilated. The optimal mixture (with 6
components) is highlighted with a green dot. (c) Evolution of the number of clusters. (d) The points projected onto
the tangent space of one specific cluster from the solution mixture. Each point is colored by the value of Eq. (18).
The cluster on the opposite side of the point the tangent space is centered on is seen to be spread around the cut
locus, which is a circle of radius π. Best viewed in color.
1400

15

Clusters

Cost

1300
1200
1100

10

5

1000
900

0
50

(a) Input Data(e) Solution Cluster #2

100 150 200

50

100 150 200

Iteration

Iteration(b) Cost Function(c) Number of Clusters(f) Solution Cluster #3(g) Solution Cluster #4(d) Solution Cluster #1(h) Solution Cluster #5

Fig. 7: Quadratic Surface Example. (a) Section of the manifold with input data generated from 1500 points
sampled from a mixture of 5 Gaussian distributions. (b), (c) refer to Fig. 6. All retrieved clusters are shown in (d)to (h). Best viewed in color.
the algorithm is initialized with 15 clusters and uses
the generic forms of the expp and logp operators that
rely on the derivative of the implicit manifold equation
as detailed in [14,61]. Additionally, the threshold s described in Sect. 3.6 is set to 0.999 to speed up the computations. In this scenario the underlying distribution
is recovered as shown in Fig. 7.
In these synthetic experiments, we also analyze theeffects of the non-linearities of the expp and logp operators by evaluating the method on a sphere (Eq. (35)with c = 1) using 6 clusters with mean and covariance parameters µi = [cos(iπ/3), 0, sin(iπ)/3] and Σi =diag(σ, σ) for i = 1, . . . , 6, and for increasing valuesof σ. Several examples of input distributions with different covariances are shown in Fig. 8a. The effect of
the number of input samples is seen by testing with
N = {600, 1800, 6000}. The algorithm parameters used
are the same as in the aforementioned sphere example.
The results are shown in Fig. 8b-d. With less spread
Gaussians and little overlap, having more data can be

12
25
N=600
N=1800N=600020

σ = 0.15

σ = 0.30

0.15(a) Example Input Data

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

5

σ = 0.60

1

15

10

σ = 0.45

1

0.30

0.45

0.60

0
0.15

0.30

0.45

0.60

0.15

0.30

0.45

0.60

σ

σ

σ

(b) Mean Clusters(c) K = 6 Solutions(d) K ≤ 9 Solutions

Fig. 8: Effect of cluster size and number of samples. Evaluation of the proposed algorithm for increasing
covariance sizes diag(σ, σ) and three amounts of samples points N on the sphere S 2 . For each combination of theseparameters we run 100 different experiments and we report the average results. (a) Examples of the different input
data distributions we consider. (b) Mean number of solution clusters found. (c), (d) Ratio of estimated solutions
with a number of clusters subject to different constraints. Best viewed in color.
seen to be beneficial. However, with more overlap andsamples, generally the data gets partitioned into moreclusters. These results seem to indicate that the algorithm tends to implicitly favor smaller Gaussians overlarger ones, suggesting that there shouldn’t be problems with approximating distributions. It is also worth
mentioning that these results are for clustering. Whenestimating a probability density function, the numberof clusters is not particularly important as long as theunderlying density is properly approximated.
Finally, in order to evaluate the benefit of using multiple tangent spaces over a single one, we perform acomparison on the sphere manifold, in two situations:the same 6 clusters as in Fig. 8 with Σi = diag(0.2, 0.3),
and when fully covering the sphere with two additional
clusters centered at (0, 1, 0) and (0, −1, 0). We also compare against an approach that uses von Mises-Fisher(vMF) distributions [3], which is specifically designed tocluster data on a sphere, and two approaches based on
Dirichlet Processes (DP-GMM and DP-TGMM). DPGMM uses split/merge proposals with parallel sampling in order to estimate a mixture model. This wasthen extended to sphere manifolds to use multiple tangent spaces (DP-TGMM). For the vMF distributions
we use our own improved implementation based on [17](see Appendix B), while we use the authors implementation of DP-GMM and DP-TGMM from [62]. We use
the default parameters for all approaches except DPGMM and DP-TGMM in which we augment the number of allowed iterations to improve their results.
The results are summarized in Table 2. In the 6cluster case our algorithm retrieves the correct numberof clusters in a 92% of the experiments, while one single
tangent plane only provides a 46% of success. Note that
we evaluate the performance of the methods based only

on the number of clusters, and not comparing the entire density probability. In the following subsection wewill show that the distributions obtained with one single
tangent plane are also much less representative as thoseobtained with the proposed approach. In the 8-cluster
case our algorithm’s performance improves and it always finds the correct clusters, while a single tangent
space always fails, with an average of 15 clusters found.
This is likely caused by the fact that the 8-clusters areevenly distributed around the sphere causing a single
tangent space to suffer from extreme deformation. This
amount of deformation on the contrary helps our optimization process to place mixtures, and thus tangent
spaces around the sphere exploring more of the solution space and finding the near optimal solution. On
the other hand, the 6-clusters do not have such a large
amount of deformation and our approach is sometimes
unable to properly locate the cluster means. Using vMF
distributions results in an oversegmentation of the datain both experiments. This is due to the fact that thevMF distributions use a single parameter κ for the concentration of the data, while our model allows for much
more expressive covariances in the form of matrices.
The Dirichlet process-based approaches show promising results, especially the sphere-specific approach DPTGMM, which is able to find the correct distribution
roughly half of the times. However, in all cases this approach is outperformed by our approach. These results
clearly show the advantage of using multiple tangent
planes to better approximate manifold distributions. In
the next subsection we will evaluate more into detailour method on this manifold.

13

5.2 Estimating Distributions on the S 2 sphere
In order to evaluate the different hyperparameters of
our model, we perform much more in depth analysis ofclustering on the S 2 sphere manifold. In particular, we
focus on the influence of the number of samples used tolearn the distribution, the specific covariance estimation
algorithm, k-means initialization, and the sample noise.
We consider a S 2 sphere with a synthetic distribution
formed by 6 clusters as shown in Fig. 6a. As a metric
we will consider the log-likelihood of 60,000 test samples, randomly generated from the distribution (10,000
samples per cluster). ThatPis, for each test sample x wefirst compute l(x) = log( k αk p(x|Θk )), and then wereport the average log-likelihood for all samples. The
larger this value, the better the samples fit the estimated distribution.
Additionally, unless otherwise specified, for the restof parameters we use the same values as in the previous subsection. For all experiments, we randomly sample from the distributions 100 times to obtain different
training sets and report the mean and standard deviation of the results for all 100 learned models, each evaluated on the 60,000 test samples.
In our first experiment we simultaneously considered the effect of the number of training samples, typeof covariance estimator and whether or not k-means isused for initialization. We report the results in Table 3.
We can see that the shrinkage-based estimators outperform the empirical-based ones, especially when thenumber of training samples is small compared to thedimensionality of the manifold. Also note that our approach consistently outperforms methods relying on asingle tangent space and the approach using von MisesFisher distributions. Regarding the use of k-means for
initialization, we did not observe that much difference,
likely due by the low-dimensionality of the manifold.
Finally, no matter the approach used, performance degrades with decreased amounts of training data.
We also looked into the effect of sample noise whenclustering. Even though the true sample noise is notknown, there are applications in which it may be possible to obtain it, such as when clustering data obtained
from sensors with known properties. In particular welook at the extreme case of small amounts of training
data given the manifold dimensionality. We summarize
the results of this analysis in Table 4. Note that all approaches in this case benefit from this added sample
noise until it becomes too large. The best result is obtained with Ledoit-Wolf covariance shrinkage approachand Σxi = 10−2 sample noise (the units of this noise
could be interpreted as an angle in radians).

5.3 Human Pose Prior

To illustrate a practical utility of our approach, we used
it to model human pose on the Human3.6M Dataset [25,26] which has different subjects performing a variety ofactivities. We consider a simplified model of the humanbody with 15 joints, represented in a 24-dimensionalpose manifold. The corresponding block-diagonal covariance (as in [55]) has 46 non-zero elements, and thefull covariance matrix has 576 non-zero elements.
We split the dataset by using all the 15 categories ofactions, each comprising two subcategories, for actors
5, 6, 7, 8, 9, and 11 for the training set, and use actor 1
as the test set. We perform an in-depth analysis of ourmodel using this split, and finally show results on generalization using other actors. The diversity of the actions
makes the dataset very challenging to learn. This gives
us 465,325 frames for training and 62,064 frames fortesting. Since the frames are highly correlated because
motions are smooth, we perform a random subsampling
before training our model. To assess its influence, wewill consider four different subsampling levels.
We consider three baselines: standard Gaussian mixture model directly on 3D joint positions, using a single tangent space for clustering (1-TM), and clustering with von Mises-Fisher distributions (vMF). We also
compare five different variants of our model. The blockdiagonal covariance matrix approach of [55], without
covariance shrinkage estimators and with them. Furthermore, we observed that when not using block matrices, the empirical covariance estimation run into numerical issues, while shrinkage estimators performed stable.
Therefore, we also report results for non-block-diagonalcovariance matrices. All these cases are summarized inTable 5.
Results show that, as expected, the non-manifold
based approach (GMM) fails. The one tangent planemodel (1-TM) performs poorly due to the non-linearitiesof the approximation, which gets worse with more training data. The von Mises-Fisher performs consistently
better. We see that our models perform best, althoughwith too much heavily correlated training data (larger
subsampling ratios) they tend to suffer from overfitting
and performance degrades. In particular, the best performance is obtained using the Oracle Approximating
Shrinkage algorithm for covariance estimation with 5%of the training data. This model greatly outperforms
the previous work and confirms the importance of morerobust covariance estimation algorithms.

14

Cov. Est.
k-means

Ours
LW [33]

Empirical
No
Yes

No

Yes

OAS [11]
No
Yes

1-TM
Empirical
No

vMF
vMF
No

6000 training samples

mean
std

-0.7474
0.0044-0.7474
0.0041-0.7466
0.0035-0.7466
0.0035-0.7473
0.0041-0.7472
0.0043-10.9331
0.6564-1.77220.0040600 training samples

mean
std

-0.8092
0.0279-0.8097
0.0288-0.7773
0.0127-0.77780.0128-0.80120.0239-0.7953
0.0214-10.3445
0.8049-1.8386
0.014560 training samples

mean
std

-1.48200.3094-1.4826
0.2795-1.1064
0.1121-1.10310.1095-1.18290.1702-1.1849
0.1434-10.0075
2.0754-2.2944
0.1662

Table 3: Comparison of the effect of number of training samples, k-means initialization, and different
covariance estimators. Evaluation of various settings on a synthetic distribution on the S 2 sphere using the loglikelihood of 60,000 random test samples. The best result for each fixed number of training samples is highlighted
in bold.
Cov. Est.
k-means

Sample noise

Empirical
No
Yes

LW [33]No
Yes

OAS [11]
No
Yes

0

mean
std

-1.48200.3094-1.4826
0.2795-1.1064
0.1121-1.10310.1095-1.18290.1702-1.1849
0.1434

10−3

mean
std

-1.3673
0.2038-1.40930.2529-1.08900.1071-1.08110.1050-1.1473
0.1270-1.1572
0.1423

10−2

mean
std

-1.0877
0.0987-1.0916
0.0979-1.0265
0.0730-1.0157
0.0609-1.0138
0.0781-1.0057
0.0772

10−1

mean
std

-1.2261
0.0395-1.2287
0.0399-1.2374
0.0437-1.24400.0440-1.2595
0.0432-1.26210.0437

Table 4: Comparison of the effect of the sample noise, k-means initialization, and different covarianceestimators. Evaluation on a synthetic distribution on the S 2 sphere with 60 training samples (10 per cluster)using the log-likelihood of 60,000 random test samples. The best result is shown in bold.
Left Hand

Left Hip

Right Hand

0.08

0.08
0.06 60

60

Left Foot

0.025

0.05

0.02

0.0460

60

0.06
0.04
50

50

0.04

0 40

40

0.03

0.015

0.02

−0.02

30

30

0
−0.01020−0.00520−0.06−0.08 10

−0.02

0
−0.02

−0.01

0

0.01

0.02

0.03

−0.12
0.04 −0.06

−20
−0.03

−0.015

−0.04
−40

0
−0.04

−0.02

Left Hand

0

0.02

0.04

0.06

0.08

−0.02
−0.03

0.1

−0.02

Left Foot

0.5

−0.01

0

0.01

0.02

−20

−40

−0.05
−0.05

0

Right Elbow

0.08
60

0.0440

40

0.02
20

20

0.05

Right Hand

0.1

60
0.06

0.3

0.08

0.4

60

60

50

50

40

40

30

30

20

20

0.2

0.06

0.3

0.1
0.04

0.2

0
0

0.1

0
−0.02
−0.04−20

−20

0

0.02
0

−0.1

−0.02

0

−0.2
−0.06
−40

−0.04
−40

−0.1
−0.08
−0.2
−0.2

0

−0.02

−0.01
10

−0.1
−0.04
−0.03

20

0

−0.04
0

40

0.01
20

0.01
0.005

0.02

40
0.02

−0.15

−0.1

−0.05

0

0.05

0.1

0.15

−60
−0.1
−0.1

−0.3

−0.06
−60
−0.05

0

0.05

−0.08
−0.1

0.1

−0.05

0

0.05

0.1

−0.4
−0.5

10

10

0

0
−0.4

−0.3

−0.2

−0.1

0

0.1

0.2

0.3

Fig. 9: Qualitative examples. Several examples from the test set using the 15% subsampled model on the V1
manifold. We visualize the ground truth and 100 samples from our model in 3D. For visualization purposes thevelocity is scaled by a factor 10× and the samples are scaled by 3×. We also show the distribution of the samples
on the tangent space for some of the joints, scored by their log-likelihood with the ground truth as a black diamond.
5.4 Tangent bundle-based Tracking Prior
We also evaluate the proposed algorithm in a tracking
task, in which, given several consecutive frames we seek
to predict the next one. We will consider three mani-folds:

V1 = xt , logxt (xt+1 ) = T H ,
V2 = logxt (xt−1 ), xt , logxt (xt+1 ) ,
V3 = logxt (xt−2 ), logxt (xt−1 ), xt , logxt (xt+1 ) ,

(36)

15

Ours

Model

0.01

Subsampling
0.05
0.15

0.30

GMM
1-TM
vMF

-708.4-1.6533.413-708.4-2.6473.414-708.4-3.9593.412-708.4-3.975
3.416

Block
LW block
OAS block
LW
OAS

4.995
5.285
4.937
7.272
7.250

4.570
5.581
4.621
8.039
8.411

3.130
3.614
3.076
5.225
7.276

1.388
1.560
1.289
2.765
4.624

Table 5: Modeling 3D Human Pose. Loglikelihood of test data for various models. Evaluating models learned with different subsamplings ofthe training data. We consider Gaussian Mixture Models (GMM), using a single tangent plane (1-TM), vonMises-Fisher distributions (vMF), and five different
variants of our model (using block covariance matrices or not, and what type of covariance estimation algorithm is used). The best result is highlighted in bold.

where xt is the pose at frame t.
For each manifold we will estimate a mixture thatlearns p(V). Once this is learned it is then possible topredict the next pose by using the conditional distribution p(logxt (xt+1 )|V ∗ ) where V ∗ is the manifold resulting from removing logxt (xt+1 ) from V. Recall that thismarginal is indeed another mixture model.
We first perform a quantitative evaluation by looking at the log-likelihood of p(V) for the test and train
sets. As we did in the previous subsection, we subsample the training data to gain in efficiency. This is possible with no performance loss, due to the large degreeof redundancy in the training data. A subsampling percentage of 15% corresponds to 69,799 training samples,roughly the same number as the test set. The parameters are set to the same values as in the pose-onlycase. Table 6 summarizes the results, in which for eachmethod and experimental setup, we display the number of components of the estimated mixture and the loglikelihood of the train and test sets. Note that the number of estimated components increases with the amountof data. As the dimension of the manifold is increased,
it is harder to learn the models due to the additional
degrees of freedom of each covariance matrix. It is important to note that we are unable to use full covariancematrices on the V2 and V3 manifolds due to the large
number of degrees of freedom (5,184 and 9,216 respectively), while it is possible to use them with 0.15 and0.30 subsampling on the V1 manifold (2,304 degrees of
freedom). Furthermore, it is only possible to estimate
the covariance of such large matrices with shrinkage
covariance estimators. Using a full covariance matrix
provides a large increase in performance over the block-

diagonal when applicable. Additionally, adding more
temporal information increases performance, although,more data is needed to learn these models.
We next evaluate the model to predict future positions, that is, the log-likelihood of p(logxt (xt+1 )|V ∗ ), instead of the global likelihood p(V). We compare against
a Gaussian Diffusion (GD) approach both trained on a
global level (a single Gaussian is averaged for all joints)
and on a local level (a single Gaussian is averaged foreach joint independently)2 . Recall that both these approaches consists of simply defining the motion as a
Gaussian distribution centered on the previous frame,
and they operate directly on the Euclidean space, andnot on the manifold. We train several kinematic models with different degrees of subsampling of the training
data, and report the results in Table 7. We can seethat the Local Gaussian Diffusion (LGD) model outperforms the standard Gaussian Diffusion (GD) model.
Yet, our model outperforms both of them by a considerable margin. These results largely agree with Table 6.
In order to assess the generalization capability of thealgorithm, we evaluate our approach with different subject splits and summarize the results in Table 8. We
use the leave-one-out strategy: using all the subjects
except one for training, which is used for testing. We
evaluate as many times as there are subjects, changing the subject which is being left out for testing eachtime. For fairness with the Gaussian Diffusion (GD) approaches that only account for pose (and not velocity)
information, we just consider the V1 manifold. For all
approaches a 0.15 subsampling ratio is used. Regarding
our approach, we use it both with block diagonal matrices, as in [56], and with the improved version basedon the Oracle Approximating Shrinkage (OAS). Observe that the latter yields a large performance gain.
On average, both alternatives outperform the GaussianDiffusion baselines, except for subject 7. In this dataset
in particular, the actors were given a lot of freedom toperform the actions. It is likely that subject’s 7 motion
largely deviates from other subjects. It is also interesting to note that subjects 1, 8, and 11 have better performance on the test set rather than the train set. Thisis likely due to the fact that there is correlation acrosssubjects.
Finally, we show some qualitative examples in Fig. 9.
For this, we directly sample from the conditional distribution for several frames. It is worth noting that we
can obtain 100,000 samples in 0.85 seconds on a Intel
Core i7 2.93GHz CPU using a Matlab implementation.
2 Since vMF-distributions are not directly applicable to predicting velocities on the hypersphere we do not include them inthis experiment.

16
Subsampling

0.05

0.15

0.30

Manifold

Method

#Mix

Train

Test

#Mix

Train

Test

#Mix

Train

Test

V1

Block
LW block
OAS block
LW
OAS

38
25
40
-

18.70816.34718.701-

16.76715.197
16.773
-

173
161
177
9
8

24.23222.41724.38224.03625.97617.28116.99117.50220.08722.134

365
335
370
22
22

28.27626.48628.31828.26731.84217.03017.251
16.93521.71124.105V2Block
LW block
OAS block

19
17
17

25.87424.52225.18430.24928.98229.58775
72
78

32.44731.34832.66934.69433.81234.807

180
170
167

36.54535.31336.07835.27734.83435.133V3Block
LW block
OAS block

13
11
11

17.82516.01316.77227.29425.10625.1064337
47

26.57524.75927.07534.81933.00434.95998
84
98

31.38329.56531.17337.48936.47837.085

Table 6: Log-likelihood of joint pose and kinematic models p(V). We evaluate on several different manifolds
for various degrees of subsampling of the training data. For each case, we plot the number of estimated mixture
components (#Mix), and the log-likelihood values for the train and test sets. Testing is performed on subject 1
while training is performed on the rest of the subjects.

Subsampling
0.15
0.30

Manifold

Method

0.05

-

GD
LGD

5.450
6.428

5.425
6.412

5.431
6.410

V1

Block
LW block
OAS block
LW
OAS

11.84210.50311.713-11.73911.097
11.82912.11513.474

11.49310.99511.36312.93714.100V2Block
LW block
OAS block14.71514.00814.39816.50115.99316.64416.81916.40216.725Block
LW block
OAS block14.02413.44513.591

16.57016.05316.51717.34117.05017.259V3

Table 7: Evaluation of velocity estimationp(logxt (xt+1 )|V ∗ ). Evaluation on different manifolds
and subsampling levels of the training data. We report
the log-likelihood of the velocity prediction of the testing set. Testing is done on subject 1 while training isdone on the rest of the subjects. Gaussian Diffusion(GD) and Local Gaussian Diffusion (LGD) are used asbaselines. Both approaches operate on the Euclidean
space and not on the manifold.

validated on various synthetic examples that highlighttheir ability to both correctly approximate manifold
distributions and discover the underlying data structure. Furthermore, the approach has been tested on alarge and complex dataset, where it is shown to outperform the traditionally used Euclidean Gaussian Mixture
Model, von Mises-Fisher distributions and an approach
using a single tangent space.
As a particular example, we have deeply studied
the use of the model as a 3D pose tracking prior, andhave shown it greatly outperforms the standard Gaussian diffusion prior. Additionally, by using shrinkage
covariance estimation algorithms we are able to gain
both robustness to poor data, and use more expressive
covariance matrices.
Future works include exploiting the proposed algorithm on different manifolds and datasets. We have presented results using Gaussian distributions and havefocused on the S 2 manifold. However, the algorithm
presented here should work with any distribution andon any manifold for which the exponential and logarithmic map operators are provided, as shown on aquadratic surface. For example, it could be possible
to initially estimate unknown and non-parameterizable
manifolds [5], and use approximate operators [19].

6 Conclusions
7 Acknowledgments
We have presented a novel data-driven approach formodeling the probability density function of data located on a Riemannian manifold. By using a mixtureof distributions, each with its own tangent space, we areable to ensure the consistency of the model while avoiding most of the linearization error caused by a single
tangent space. The approach has been experimentally

We would like to thank the three anonymous reviewers
for their insights and comments that have significantly
contributed to improving this manuscript. This work
was partly funded by the Spanish MINECO project
RobInstruct TIN2014-58178-R and by the ERA-net CHISTERA project I-DRESS PCIN-2015-147.

17

MethodSubject 1
Train Test

Subject 5Train Test

Subject 6
Train Test

Subject 7
Train Test

Subject 8Train Test

Subject 9
Train Test

Subject 11
Train Test

Mean
Train Test

GD
LGD

5.433 5.435
6.409 6.410

5.453 5.451
6.421 6.421

5.452 5.451
6.453 6.452

5.509 5.502
6.516 6.511

5.408 5.413
6.375 6.380

5.442 5.442
6.443 6.444

5.410 5.415
6.409 6.411

5.444 5.444
6.432 6.433

Ours Block 9.360 11.925 10.170 7.449 10.503 5.039 11.319 2.129 8.943 12.188 9.854 9.597
9.495 11.793 9.949 8.589
Ours OAS 12.458 14.099 13.020 10.808 13.718 8.889 14.236 5.740 12.245 17.211 12.624 12.599 12.439 15.485 12.963 12.119

Table 8: Generalization of velocity estimation results. Log-likelihood evaluation of p(logxt (xt+1 )|V ∗ ) on asingle subject for testing and the rest for training. The subject indicated in the table is the subject used for testing.
We compare a Gaussian Diffusion (GD) strategy, a Local Gaussian Diffusion (LGD) against both our approaches.
For fairness we only consider the V1 manifold, as Gaussian Diffusion approaches do not use additional velocity
information.
References
1. Andriluka, M., Roth, S., Schiele, B.: Monocular 3D Pose Estimation and Tracking by Detection. In: IEEE Conference
on Computer Vision and Pattern Recognition (2010)
2. Archambeau, C., Verleysen, M.: Manifold constrained finite
gaussian mixtures. In: Computational Intelligence and Bioinspired Systems, pp. 820–828. Springer (2005)
3. Banerjee, A., Dhillon, I.S., Ghosh, J., Sra, S., Ridgeway, G.:
Clustering on the unit hypersphere using von mises-fisher
distributions. Journal of Machine Learning Research 6(9),1345–1382 (2005)
4. Boothby, W.M.: An Introduction to Differentiable Manifolds
and Riemannian Geometry. Revised 2nd Ed. Academic, NewYork (2003)
5. Brand, M.: Charting a Manifold. In: Neural Information Processing Systems, pp. 961–968 (2003)
6. Brubaker, M.A., Salzmann, M., Urtasun, R.: A Family of
MCMC Methods on Implicitly Defined Manifolds. Journalof Machine Learning Research 22, 161–172 (2012)
7. do Carmo, M.P.: Riemannian Geometry. Birkhäuser (1992)
8. Caseiro, R., Martins, P., Henriques, J.F., Batista, J.: A nonparametric riemannian framework on tensor field with application to foreground segmentation. Pattern Recognition
45(11), 3997–4017 (2012)
9. Caseiro, R., Martins, P., Henriques, J.F., Leite, F.S., Batista,
J.: Rolling riemannian manifolds to solve the multi-class classification problem. In: IEEE Conference on Computer Vision
and Pattern Recognition (2013)
10. Chang, J., Fisher III, J.W.: Parallel sampling of dp mixture
models using sub-cluster splits. In: Neural Information Processing Systems, pp. 620–628 (2013)
11. Chen, Y., Wiesel, A., Eldar, Y., Hero, A.: Shrinkage algorithms for mmse covariance estimation. IEEE Transactions
on Signal Processing 58(10), 5016–5029 (2010)
12. Darling, R.: Martingales on noncompact manifolds: maximal
inequalities and prescribed limits. Annales de l’IHP Probabilités et statistiques 32(4), 431–454 (1996)
13. Davis, B.C., Bullitt, E., Fletcher, P.T., Joshi, S.: Population
Shape Regression from Random Design Data. In: International Conference on Computer Vision (2007)
14. Dedieu, J.P., Nowicki, D.: Symplectic methods for the approximation of the exponential map and the newton iteration on riemannian submanifolds. Journal of Complexity
21(4), 487 – 501 (2005)
15. Dempster, A.P., Laird, N.M., Rubin, D.B.: Maximum likelihood from incomplete data via the em algorithm. Journalof the royal statistical society. Series B (methodological) pp.
1–38 (1977)

16. Deutscher, J., Reid, I.: Articulated body motion capture bystochastic search. International Journal of Computer Vision
61(2), 185–205 (2005)
17. Figueiredo, M., Jain, A.: Unsupervised Learning of Finite
Mixture Models. IEEE Transactions Pattern Analylis andMachine Intelligence 24(3), 381–396 (2002)
18. Fletcher, P., Lu, C., Pizer, S., Joshi, S.: Principal Geodesic
Analysis for the Study of Nonlinear Statistics of Shape. IEEE
Transactions on Medical Imaging 23(8), 995–1005 (2004)
19. Freifeld, O., Black, M.J.: Lie Bodies: A Manifold Representation of 3D Human Shape. In: European Conference on
Computer Vision (2012)
20. Gall, J., Rosenhahn, B., Brox, T., Seidel, H.P.: Optimizationand filtering for human motion capture. International Journal of Computer Vision 87, 75–92 (2010)
21. Harandi, M., Sanderson, C., Hartley, R., Lovell, B.: Sparse
coding and dictionary learning for symmetric positive definite matrices: A kernel approach. In: European Conference
on Computer Vision (2012)
22. Harandi, M.T., Salzmann, M., Hartley, R.: From manifold
to manifold: Geometry-aware dimensionality reduction for
spd matrices. In: European Conference on Computer Vision(2014)
23. Hauberg, S., Sommer, S., Pedersen, K.S.: Natural metrics
and least-committed priors for articulated tracking. Imageand Vision Computing 30(6), 453–461 (2012)24. Huckemann, S., Hotz, T., Munk, A.: Intrinsic Shape analysis:
Geodesic PCA for Riemannian Manifolds Modulo Isometric
Lie Group Actions. Statistica Sinica 20, 1–100 (2010)
25. Ionescu, C., Li, F., Sminchisescu, C.: Latent Structured Models for Human Pose Estimation. In: International Conference
on Computer Vision (2011)
26. Ionescu, C., Papava, D., Olaru, V., Sminchisescu, C.: Human3.6m: Large scale datasets and predictive methods for 3d
human sensing in natural environments. IEEE Transactions
Pattern Analylis and Machine Intelligence 36(7), 1325–1339(2014)
27. Jain, S., Govindu, V.: Efficient higher-order clustering on
the grassmann manifold. In: International Conference on
Computer Vision (2013)
28. Jayasumana, S., Hartley, R., Salzmann, M., Li, H., Harandi,
M.: Kernel methods on the riemannian manifold of symmetric positive definite matrices. In: IEEE Conference on Computer Vision and Pattern Recognition (2013)
29. Jayasumana, S., Hartley, R., Salzmann, M., Li, H., Harandi,
M.: Kernel Methods on Riemannian Manifolds with Gaussian RBF Kernels. IEEE Transactions Pattern Analylis andMachine Intelligence (2015)
30. Karcher, H.: Riemannian center of mass and mollifiersmoothing. Communications on Pure and Applied Mathematics 30(5), 509–541 (1977)

18
31. Lawrence, N.D.: Probabilistic Non-linear Principal Component Analysis with Gaussian Process Latent Variable Models.
Journal of Machine Learning Research 6, 1783–1816 (2005)
32. Lawrence, N.D., Moore, A.J.: Hierarchical Gaussian Process
Latent Variable Models. In: International Conference in Machine Learning (2007)
33. Ledoit, O., Wolf, M.: A well-conditioned estimator for largedimensional covariance matrices. Journal of Multivariate
Analysis 88(2), 365–411 (2004)
34. Ledoit, O., Wolf, M.: Nonlinear shrinkage estimation of largedimensional covariance matrices. Institute for Empirical Research in Economics University of Zurich Working Paper(515) (2011)
35. Lenglet, C., Rousson, M., Deriche, R., Faugeras, O.: Statistics on the manifold of multivariate normal distributions:
Theory and application to diffusion tensor mri processing.
Journal of Mathematical Imaging and Vision 25(3), 423–444(2006)
36. Li, R., Tian, T.P., Sclaroff, S., Yang, M.H.: 3d human motion tracking with a coordinated mixture ofăfactor analyzers.
International Journal of Computer Vision 87(1-2), 170–190(2010)
37. Moeslund, T.B., Granum, E.: A survey of computer visionbased human motion capture. Computer Vision and Image
Understanding 81(3), 231–268 (2001)
38. Moeslund, T.B., Hilton, A., Krüger, V.: A survey of advancesin vision-based human motion capture and analysis. Computer Vision and Image Understanding 104, 90–126 (2006)
39. Muralidharan, P., Fletcher, P.T.: Sasaki metrics for analysis
of longitudinal data on manifolds. In: IEEE Conference on
Computer Vision and Pattern Recognition (2012)
40. Ozakin, A., Gray, A.: Submanifold Density Estimation.
In: Neural Information Processing Systems, pp. 1375–1382(2009)
41. Pelletier, B.: Kernel Density Estimation on Riemannian
Manifolds. Statistics & Probability Letters 73(3), 297 – 304(2005)
42. Pennec, X.: Intrinsic Statistics on Riemannian Manifolds:
Basic Tools for Geometric Measurements. J. Math. Imaging Vis. 25(1), 127–154 (2006)
43. Pennec, X.: Statistical computing on manifolds: from riemannian geometry to computational anatomy. In: Emerging
Trends in Visual Computing, pp. 347–386. Springer (2009)44. Pennec, X., Fillard, P., Ayache, N.: A Riemannian framework
for tensor computing. International Journal of Computer
Vision 66(1), 41–66 (2006)
45. Quiñonero-candela, J., Rasmussen, C.E., Herbrich, R.: A
Unifying View of Sparse Approximate Gaussian Process Regression. Journal of Machine Learning Research 6, 1939–1959(2005)
46. Said, S., Courtry, N., Bihan, N.L., Sangwine, S.: Exact Principal Geodesic Analysis for data on SO(3). In: European
Signal Processing Conference (2007)
47. Sanin, A., Sanderson, C., Harandi, M., Lovell, B.: K-tangent
spaces on riemannian manifolds for improved pedestrian detection. In: International Conference on Image Processing(2012)
48. Sasaki, S.: On the differential geometry of tangent bundles ofriewannian manifolds. Tohoku Mathematical Journal, Second Series 10(3), 338–354 (1958)
49. Schäfer, J., Strimmer, K.: A shrinkage approach to large-scalecovariance matrix estimation and implications for functional
genomics. Statistical applications in genetics and molecularbiology 4(1) (2005)
50. Shirazi, S., Harandi, M., Sanderson, C., Alavi, A., Lovell,B.: Clustering on grassmann manifolds via kernel embedding
with application to action analysis. In: International Conference on Image Processing (2012)

51. Sigal, L., Bhatia, S., Roth, S., Black, M., Isard, M.: Tracking loose-limbed people. In: IEEE Conference on Computer
Vision and Pattern Recognition (2004)
52. Sigal, L., Isard, M., Haussecker, H.W., Black, M.J.: Looselimbed People: Estimating 3D Human Pose and Motion Using Non-parametric Belief Propagation. International Journal of Computer Vision 98(1), 15–48 (2012)
53. Simo-Serra, E., Quattoni, A., Torras, C., Moreno-Noguer, F.:
A Joint Model for 2D and 3D Pose Estimation from a Single
Image. In: IEEE Conference on Computer Vision and Pattern Recognition (2013)
54. Simo-Serra, E., Ramisa, A., Alenyà, G., Torras, C., MorenoNoguer, F.: Single Image 3D Human Pose Estimation from
Noisy Observations. In: IEEE Conference on Computer Vision and Pattern Recognition (2012)
55. Simo-Serra, E., Torras, C., Moreno-Noguer, F.: Geodesic Finite Mixture Models. In: British Machine Vision Conference(2014)
56. Simo-Serra, E., Torras, C., Moreno-Noguer, F.: Lie AlgebraBased Kinematic Prior for 3D Human Pose Tracking. In:
International Conference on Machine Vision Applications(2015)
57. Sivalingam, R., Boley, D., Morellas, V., Papanikolopoulos,
N.: Tensor sparse coding for region covariances. In: European
Conference on Computer Vision (2010)
58. Sminchisescu, C., Triggs, B.: Estimating Articulated Human
Motion With Covariance Scaled Sampling. International
Journal of Robotics Research 22(6), 371–391 (2003). Special issue on Visual Analysis of Human Movement
59. Sommer, S.: Anisotropic distributions on manifolds: Template estimation and most probable paths. Information Processing in Medical Imaging. Lecture Notes in Computer Science, Springer Berlin Heidelberg (2015)60. Sommer, S., Lauze, F., Hauberg, S., Nielsen, M.: Manifold
Valued Statistics, Exact Principal Geodesic Analysis and the
Effect of Linear Approximations. In: European Conference
on Computer Vision (2010)
61. Sommer, S., Lauze, F., Nielsen, M.: Optimization overgeodesics for exact principal geodesic analysis. Adv. Comput. Math. 40(2), 283–313 (2014)62. Straub, J., Chang, J., Freifeld, O., Fisher III, J.W.: A dirichlet process mixture model for spherical data. In: International
Conference on Artificial Intelligence and Statistics (2015)
63. Taylor, G., Sigal, L., Fleet, D., Hinton, G.: Dynamical binary latent variable models for 3d human pose tracking. In:
IEEE Conference on Computer Vision and Pattern Recognition (2010)
64. Tosato, D., Farenzena, M., Cristani, M., Spera, M., Murino,
V.: Multi-class classification on riemannian manifolds forvideo surveillance. In: European Conference on Computer
Vision, pp. 378–391 (2010)
65. Tosato, D., Spera, M., Cristani, M., Murino, V.: Characterizing humans on riemannian manifolds. IEEE Transactions
Pattern Analylis and Machine Intelligence 35(8), 1972–1984(2013)
66. Tournier, M., Wu, X., Courty, N., Arnaud, E., Reveret, L.:
Motion compression using principal geodesics analysis. Computer Graphics Forum 28(2), 355–364 (2009)
67. Turaga, P., Veeraraghavan, A., Srivastava, A., Chellappa, R.:
Statistical computations on grassmann and stiefel manifolds
for image and video-based recognition. IEEE Transactions
Pattern Analylis and Machine Intelligence 33(11), 2273–2286(2011)
68. Tuzel, O., Porikli, F., Meer, P.: Pedestrian Detection via
Classification on Riemannian Manifolds. IEEE Transactions
Pattern Analylis and Machine Intelligence 30(10), 1713–1727(2008)

19
69. Urtasun, R., Fleet, D.J., Fua, P.: 3D People Tracking with
Gaussian Process Dynamical Models. In: IEEE Conference
on Computer Vision and Pattern Recognition (2006)
70. Urtasun, R., Fleet, D.J., Lawrence, N.D.: Modeling humanlocomotion with topologically constrained latent variable
models. In: Proceedings of the 2nd Conference on Human
Motion: Understanding, Modeling, Capture and Animation(2007)
71. Varol, A., Salzmann, M., Fua, P., Urtasun, R.: A Constrained
Latent Variable Model. In: IEEE Conference on Computer
Vision and Pattern Recognition (2012)
72. Wallace, C.S., Freeman, P.R.: Estimation and inference bycompact coding. Journal of the Royal Statistical Society.
Series B (Methodological) pp. 240–265 (1987)
73. Wang, J., Fleet, D., Hertzmann, A.: Gaussian process dynamical models. In: Neural Information Processing Systems(2005)
74. Yao, A., Gall, J., Gool, L.V., Urtasun, R.: Learning probabilistic non-linear latent variable models for tracking complex
activities. In: Neural Information Processing Systems (2011)
75. Zhang, M., Fletcher, P.T.: Probabilistic principal geodesic
analysis. In: Neural Information Processing Systems, pp.
1178–1186 (2013)

A.2 M-stepIn this step, we have fixed w and are updating the other parameters θ = (µ, Σ) and α by solving:arg maxB(x, θ)θ,αKX

subject to(40)

αk = 1i=1αk ≥ 0, k = 1, . . . , K .
We shall follow the same approach as in the E-step and computethe partial derivatives to obtain the update rules. In particular,
both α and Σ are straight forward to compute, and do not significantly deviate from the standard formulation. Thus for α we
obtain:Nαk (t) =1 X (i)
wk
wk =
,
N
N

(41)

i

and for Σ:N
∂B(x, θ)1 X (i) =
wk logµk (x(i) ) logµk (x(i) )> − Σk ,
∂Σk
2

(42)

i=1

A Derivation of Mixture Models on
Riemannian Manifolds

and thus,

We follow the standard expectation-maximization approach tomaximize the log-likelihood of our model adapting it to Riemannian manifolds. For simplicity, we will not consider the Minimum Message Length criteria for model selection. We start out
by defining the log-likelihood of the model λ(x, θ) and bounding
it by Jensen’s equality:

λ(x, θ) =

N
X
i=1

≥

K
X

log

(i)
(i)
(i) >
i=1 wk logµk (x ) logµk (x ).
PN(i)i=1 wkPN

αk p(x(i) |θk )Σk (t) =For the mean µk we can follow the same approach, however, due
to the logarithmic map, it is slightly different to resolve. We startout by computing the partial derivative:
NX∂ logµk (x(i) )∂B(x, θ)(i)=wk Σk−1 logµk (x(i) )
.
∂µk
∂µk

αk
(i)
wk logi=1 k=1p(x(i) |θ
(i)

(44)

i=1

k=1

N X
K
X

(43)

k)

In general, there is no analytic solution to= B(x, θ) .

(37)

wk

(i)

∂ logµ (x∂ logµ (x(i) )k∂µk. How-)kever, under the assumption that= c where c 6= 0 is a∂µk
constant, and equating the partial derivative to 0, we can obtain:(i)with wk as auxiliary variables that represent membership probabilities. We can maximize over the lower bound B(x, θ) insteadof the untractable full likelihood.

A.1 E-step(i)

The E-step consists of maximizing the auxiliary terms wk which
are the membership probabilities of the samples. This is done bysolving:
arg maxNX(i)wk logµk (x(i) ) = 0 .

For simply connected and complete manifolds whose curvature is
non-positive (i.e., Hadamard manifolds) and bounded from below, there exists one and only one Riemannian center of mass
which is characterized by E[logµ (x)] = 0 [12]. Note that a compact and simply connected manifold with a non-positive andbounded from below curvature has no cut locus. In this case,as Eq. (45) is the discrete expectation of the weighted sum, we
can establish the update rule for the mean by:

B(x, θ)w

subject to

(45)

i=1

N
X

(i)

wk = 1, k = 1, . . . , Kµk (t) = arg min
(38)

i=1
(i)

wk ≥ 0, i = 1, . . . , N, k = 1, . . . , K .
This is straight forward to do by computing the derivative andequating it to 0 to obtain the update rule for step t:
αk (t − 1)p(xi |θk (t − 1))(i)wk (t) = PK.
k=1 αk (t − 1)p(xi |θk (t − 1))

p


N
X
N
d

i=1

wk

(i)

wk x(i) , p

2.

Note that this does not hold for the case in which there is a
cut locus, in which case there may not be only one Riemannian
center of mass. However, in practice, this approach will generally
converge to the center of mass. We will use Eq. (46) in all cases.
Finally, we perform a numerical analysis of the error for the∂ logµ (x(i) )

(39)

(46)

k
S 2 sphere by numerically computing
and visualizing∂µk
the results. In particular we visualize the change of Frobenius

20

8
6
4
2
0

Fig. 10: Plot of the change Frobenius norm of∂ logµ (x(i) )k. We compute the derivative numerically us∂µk

ing a first order approximation as there is no analytic
form. We can see for points near the center there issmall change in the derivative and thus little error in
the approximation we make by considering the derivative to be constant. For visualization purposes we only
display points with a change of under 10 units.

norm of the Jacobian in Fig. 10. We can see that points near
the origin have very little change in the derivative. Again, theuse of multiple tangent planes favors configurations in which thepoints are close to the center, and thus, keeps the error produced
by approximating the partial derivative to a constant within reasonable bounds.

B Clustering with von Mises-Fisher
distributions
Given a random vector x on the unit hypersphere of dimensionq−1, the probability density function of a von Mises-Fisher distribution with mean direction µ and concentration κ can be writtenas:T

fq (x|µ, κ) =cq (κ)eκµ x
cq (κ) =,κq/2−1(2π)q/2 Iq/2−1 (κ)

,

(47)

where kµk = 1, and Iq/2−1 is the modified Bessel function of firstkind and order q/2 − 1. Note that the concentration parameterκ is a single scalar that represents a uniform distribution on thesphere for κ = 0 and is unimodal for κ > 0.
The algorithm from [17] can be modified to use von MisesFisher distributions by adapting the way the distributions arerecalculated in the M-step. This can be computed by:
rk =

N
X

(i)

wk x(i)

,

(48)

i=1

µk (t) =rkkrk kκk (t) =krk k(q − krk k2 )1 − krk k2

,

(49)
.

(50)

As there exists no analytic form of Iq/2 (κk (t))/Iq/2−1 (κk (t)) = rk ,the computation of κk (t) is indeed an approximation [3].