Universität Karlsruhe (TH)Interactive Systems Laboratory
Prof. Dr. Alexander Waibel

Massachusetts Institute of Technology
Computer Science and Artificial Intelligence Lab
Prof. Dr. Trevor Darrell

Human Body Tracking
with Rank Priors for Continuous
Non-Linear Dimensionality Reduction
Diploma Thesis

Andreas Geiger
July 24th, 2008

Supervised by:Prof. Dr. Alexander Waibel
Dr. Rainer Stiefelhagen
Dipl.-Inform. Kai Nickel

Prof. Dr. Trevor Darrell
Dr. Raquel Urtasun

Universität Karlsruhe (TH)MIT CSAIL, UC Berkeley EECSErklärung

Hiermit erkläre ich, dass ich die vorliegende Arbeit selbständig und ohne unerlaubte fremde Hilfe angefertigt, andere als die angegebenen Quellen und Hilfsmittel nicht benutzt und die den benutzten Quellen wörtlich oder inhaltlich entnommenen Stellen als solche kenntlich gemacht habe.

Karlsruhe, den 24. Juli 2008.................................................................

3

4

Abstract

Non-linear dimensionality reduction methods are powerful techniques to deal withhigh-dimensional datasets. However, they often are susceptible to local minima
and perform poorly when initialized far from the global optimum, even when
the intrinsic dimensionality is known a priori. In this work a prior over the dimensionality of the latent space is introduced, which allows for simultaneously
optimizing both the latent space and its intrinsic dimensionality. Ad-hoc initialization schemes are unnecessary with this approach; the latent space is initialized
to the observation space and the latent dimensionality is automatically inferredusing an optimization scheme that drops dimensions in a continuous fashion. The
prior is applied to various tasks involving probabilistic non-linear dimensionality
reduction, and it is shown that this method can outperform graph-based dimensionality reduction techniques as well as previously suggested ad-hoc initialization
strategies. Results are reported for artificial datasets which cause problems tonearest neighbor-based spectral methods and for a real-world example, where thetask consists of tracking a person performing activities in a kitchen scenario.

56

Zusammenfassung

Eine Vielzahl unterschiedlicher Sensoren erlaubt heute die Registrierung von Datenauf modernen Rechenmaschinen. Allerdings ist die Dimension D der Daten (zumBeispiel D = 1310720 im Falle einer 1280 × 1024-Videokamera) oft zu hoch um
sie direkt weiter verarbeiten zu können. Daher sind Methoden für nichtlineare Dimensionsreduktion unerlässlich. Solche Methoden versuchen, durch eine geeignet
gewählte Abbildung (in einen Raum niederer Dimension), die Datenmenge zu reduzieren ohne dabei die für die Aufgabe relevante Beziehungen in den Daten zu
zerstören. Dabei wird angenommen dass die Daten auf oder nahe einer Mannigfaltigkeit niederer Dimension verteilt sind, welche in einen hochdimensionalen
Raum eingebettet ist. Je nach Methodik kann dabei auch a-priori Vorwissen über
den Prozess mit eingearbeitet werden, sofern dieses zur Verfügung steht.
Verschiedene Verfahren zur Reduktion der Dimension sind in den letzten Jahren
entwickelt worden. Ein besonderes Augenmerk fällt dabei auf die nicht-linearen
Verfahren, da diese oft eine kompaktere Repräsentation erzeugen können als lineare Verfahren [Pea01]. Dies ist im Besonderen dann der Fall, wenn die Daten
Abhängigkeiten aufweisen, die sich nicht in Korrelationen äussern. Einige dieser
Verfahren bieten geschlossene Lösungen an [RS00, TSL00, TB99, WS06, BN03,
ZZ03], reagieren aber oft empfindlich auf die Wahl der Parameter oder die Datenverteilung.
Andere basieren auf der Optimierung nicht-konvexer Funktionen [Law05, LQC06]
und sind daher anfällig für lokale Minima im Optimierungsprozess.
Wieder andere setzen a-priori Wissen ein [WFH08, UD07, UFG+ 08], welches
allerdings, je nach Aufgabe, nicht immer zur Verfügung steht. Zudem wird eine
geeignete Wahl der Dimension bei den meisten Verfahren als bekannt angenommen. Diese muss zunächst, zum Beispiel durch aufwändige Vergleichsprüfung(Cross-Validation) oder andere Verfahren [mFSA07, LB05] ermittelt werden.
In dieser Diplomarbeit wird eine neue Methode vorgestellt, welche auf dem Gaussian Process Latent Variable Model (GPLVM) aufbaut. Durch einen Rank Priorüber die Anzahl der Dimensionen wird die Dimensionalität kontinuierlich (während
des Optimierungsprozesses) reduziert. So lassen sich lokale Minima des GPLVM
vermeiden und die Anzahl der Dimensionen wird automatisch ermittelt. Die resultierende kompakte Repräsentation bietet zudem eine probabilistische Abbildung
vom niederdimensionalen (latenten) Raum in den hochdimensionalen (Daten-)Raum. Dadurch lassen sich Aussagen über die Unsicherheit eines abgefragten
Datums treffen.

7

Im Abschnitt Experiments wird das entwickelte Verfahren auf 3 Fälle angewandt:
Zunächst wird illustriert, wie sich mit der Methode eine Schweizer Rolle trotz
spärlicher Abtastung noch zuverlässig entfalten lässt. Danach wird die Methodeauf Beispiele angewandt, in denen es die geeignete Dimension zu bestimmen gilt.
Abschliessend werden Bewegungen eines artikularen menschlichen Körpers eingelernt. Diese werden dann verwendet um eine Person, welche Aktivitäten in derKüche (wie beispielsweise Auswellen, Mahlen, Schneiden) durchführt, erfolgreich
zu verfolgen und die Bewegungen korrekt zu klassifizieren.

8Contents
1 Introduction

11

2 Dimensionality Reduction
2.1 Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Linear Dimensionality Reduction . . . . . . . . . . . . . . . . . .
2.3 Non-linear Dimensionality Reduction . . . . . . . . . . . . . . . .
2.3.1 Graph Based Methods . . . . . . . . . . . . . . . . . . . .
2.3.2 Latent Variable Models . . . . . . . . . . . . . . . . . . . .
2.3.3 Sparsification . . . . . . . . . . . . . . . . . . . . . . . . .

13
13
14
16
16
18
26

3 Priors for Latent Variable Model
3.1 Dynamical Priors . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Topological Priors . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Rank Priors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Rank Prior Formulation . . . . . . . . . . . . . . . . . . .
3.3.2 Derivatives . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.3 Putting it all together: the Algorithm . . . . . . . . . . . .

31
31
32
33
33
36
37

4 Experimental results
4.1 A sparsely sampled swiss roll . . . . . . . . . . . . . . . . . . . . .
4.2 Discovering the correct dimensionality . . . . . . . . . . . . . . .
4.3 Tracking and classification in the kitchen domain . . . . . . . . .
4.3.1 Learning . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.2 Feature extraction and image likelihood . . . . . . . . . . .
4.3.3 Particle filter . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

41
41
47
55
57
59
60
63

5 Implementation

67

6 Acknowledgments

71

9

Contents

10

1 Introduction
Many real-world problems involve high dimensional datasets that are computationally challenging to handle. In such cases it is desirable to reduce the dimensionality of the data while preserving the original information in the data distribution,
allowing for more efficient learning and inference. Linear dimensionality reduction methods (e.g., PCA) are efficient but miss important structure in non-linear
data. Graph-based techniques, e.g., LLE [RS00] and ISOMAP [TSL00], capture
non-linear dependencies but require highly dense and homogeneously sampled
manifolds for accurate modeling.
Non-linear dimensionality reduction techniques can be applied to more complex
data, but generally suffer from local minima. Choosing the dimensionality of thelatent space is non-trivial, and existing methods typically rely on cross-validation.
Even when given the correct latent dimensionality, these techniques often do not
succeed in practice when initialized far from the global minimum [UFG+ 08]. Factors which contribute to this include the distortion introduced by the initializationand the non-convexity of the optimization: when optimization is performed in a
low dimensional space the model may not have the requisite degrees of freedom
to avoid local minima.
In this work we develop a prior on the dimensionality of the set of latent coordinates, encouraging low dimensional representations. Our Rank Prior enforces apenalty on the non-sparsity of the singular values of the matrix of latent variables,
and automatically discovers the latent space and its dimensionality using a continuous optimization that drops dimensions on the fly. By initializing the latentcoordinates to the original space no distorsion is introduced; since we decrease thedimensionality slowly (starting from a high-dimensional space) extra flexibility isgained which allows the method to avoid local minima during optimization. To
our knowledge, ours is the first non-linear dimensionality reduction technique thatpenalizes the latent space rank and simultaneously optimizes the structure of thelatent space as well as its intrinsic dimensionality.
We demonstrate the effectiveness of our approach when learning probabilistic latent variable models with the Gaussian Process Latent Variable Model (GPLVM)[Law05], a generalization of Probabilistic PCA to the non-linear case that models
the mapping from the latent space to the data space as a Gaussian process. The
GPLVM has proven successful in many applications, but initialization, knowledgeof the latent dimensionality, and/or additional prior knowledge were assumed (e.g.,11CHAPTER 1. INTRODUCTION[GMHP04, SUF08, UFG+ 08]). Incorporating our prior with the GPLVM objective results in an optimization problem that allows us to discover the latent space
and the intrinsic dimensionality of artificial and real datasets.
We chose articulated body tracking as an application of our technique. Since afull human body model consists of many degrees of freedom (DOFs), reducing
the dimensionality is key to performing robust and reliable tracking. Since theobservation likelihood for high-dimensional models is highly non-convex, tracking
in the original space is a hard problem, especially when using a small number ofcameras or a simple observation likelihood. Furthermore the need for many evaluations of the likelihood function inhibits a reasonable tracking speed. Consider
for example a particle filter where the number of particles needed grows exponentionally with the number of dimensions. Using a prior model created via motion
capturing data helps in two ways: While making the algorithm computationallytractable it helps disambiguating ambiguous observations. We thus demonstrate
the effectiveness of our approach in tracking and classifying complex articulated
human body motions from video by using such a prior model which is explainedin more details in this thesis.

12

2 Dimensionality Reduction
A fundamental challenge of artificial intelligence is to develop useful internal representations of the external world, ideally by introducing as little prior knowledge
as possible. The human brain excels in tasks like memorizing objects or faces byextracting features from a vast amount of input information perceived by sensors
like the eyes, for instance. Since visual information is usually high-dimensional inraw form (consider the image of a 10 Megapixels consumer camera), reducing thedimensionality is key for further using the data by processing it on a computer.
Only with a good representation and relevant features is higher-level decision making possible (a high-level decision in an image could be a single binary feature,such as the answer to the question if a given object is present in the image ornot).
While features for visual input can be designed in problem-specific fashion (e.g.
by building bag-of-word features or by looking for strong gradients, corners andusing methods like SIFT [Low99]) signal data may be on hand in many different
forms and an automatic feature retrieval method is eligible.
From a statistical view, dimensionality reduction (also called dimension reduction)is the process of reducing the number of random variables under consideration.
This can be further divided into feature selection and feature extraction. While
feature selection tries to find a subset of the original variables which is mostrepresentative, feature extraction relates the high dimensional data to a lower
dimensional space while trying to maintain the information contained in the input
data. The projection can be either linear or non-linear, both approaches arepresented within this thesis.

2.1 Notation
In order to avoid confusion the following notation has been retained throughout
the paper:• N : Number of training points• D: Dimensionality of the high dimensional input (=data) space Y• Q: Dimensionality of the low dimensional latent space X

13

CHAPTER 2. DIMENSIONALITY REDUCTION• yi : input (=data) point i, yi ∈ <D
• xi : latent point i, xi ∈ <Q
• Y = (y1 , ..., yN )T : Collection of input (=data) points, Y ∈ <N ×D• X = (x1 , ..., xN )T : Collection of latent points, X ∈ <N ×Q• Yij is the entry of matrix Y at row i and column j• Xij is the entry of matrix X at row i and column j• yij = Yji is the entry of the transposed matrix YT at row i and column j• xij = Xji is the entry of the transposed matrix XT at row i and column j• x ∼ N (µ, Σ): x follows a Gaussian distribution• f (x) ∼ GP(m(x), k(x, x0 )): f (x) follows a Gaussian Process distributionP• x: mean of x: x = N1 N
i=1 xiP• y: mean of y: y = N1 N
i=1 yi

2.2 Linear Dimensionality Reduction
Linear dimensionality reduction finds a linear projection from a low dimensional
subspace X to a higher dimensional space Y. For example Principal Component
Analysis (PCA) [Pea01, DHS00], sometimes also called the Karhunen-Loève transform, seeks an optimal affine subspace X of Y which maximizes variance (along
the dimensions of X) and minimizes the reprojection error when projecting backfrom X to Y. More formally let the projection from the latent space X to thedata space Y beŷi = y +QX

xij ej(2.1)

j=1

where {ej , ...eQ } spans a Q-dimensional orthonormal basis in the data space Y.
Then the least squares reprojection error is defined as

E=

N
X

kŷi − yi k2 =

i=1

N
X
i=1

ky +

Q
X

xij ej − yi k2 .

j=1

Considering that {ej , ...eQ } forms a orthonormal basis, we can expand E to

E=

" Q
N
X
X
i=1

#
x2ij + 2xij eTj (y − yi )



j=1

14

+ ky − yi k

2

.

(2.2)

CHAPTER 2. DIMENSIONALITY REDUCTION
Minimizing E with respect to xij can be done in closed form and leads to∂E!
= 2xij + 2eTj (y − yi ) = 0∂xij⇒ xij = −eTj (y − yi ).
Applying this result to equation 2.2 we have

E =

N
X

"
−

#
x2ij + ky − yi k2

(2.3)

j=1

i=1

= −

Q
X

Q
X

eTj Sej +

j=1

N
X

ky − yi k2

(2.4)

i=1

where S is the scatter matrix of Y, defined as

S=

N
X

(y − yi )(y − yi )T .

i=1

The constraint kei k = 1 is enforced by introducing Lagrange Multipliers λ into
the target equation:

0

E =−

Q
X

eTj Sej +

j=1

N
X

2

ky − yi k +

i=1

Q
X
j=1

Minimizing E 0 with respect to ej is straightforward∂E 0!
= −2Sej + 2λej = 0∂ejleading to the eigenvalue problem
Sej = λj ej .
Since
eTj Sej = λj eTj ej

15

λj (eTj ej − 1)

(2.5)

CHAPTER 2. DIMENSIONALITY REDUCTION
E 0 is minimized by selecting ej as the top Q eigenvectors of S. Once the basis{ej , ...eQ } has been calculated, projections in both directions (from X to Y andvice versa) are readily given (since we assumed linearity in equation 2.1) as

xi = ET (yi − y)andyi = Exi + ywhere E is a D × Q matrix comprising the basis vectors ej as columns.

2.3 Non-linear Dimensionality Reduction
Linear Dimensionality Reduction (LDR) is often used due to its simplicity andclosed form solution. When data dependencies are mostly expressed in terms ofcorrelations this leads to reasonable results. Considering real world data (and also
very simple artificial examples), linearity cannot be assumed in general. Consider
for example a one dimensional manifold1 coiled up to a spiral when embeddedin two dimensions (see figure 3.2). Applying PCA in this case leads to a badsolution since points which are far away in the two dimensional input space Y areprojected to nearby points on the one dimensional subspace X (in the figure, thismeans that red points (circles) and yellow points (crosses) are not separated after
applying 1D PCA).
To handle non-linear dependencies one has to resort to Non-Linear Dimensionality
Reduction (NLDR) techniques which are discussed in this section of the thesis.

2.3.1 Graph Based Methods
Graph based methods allow for non-linear dimensionality reduction by taking into
consideration a k-nearest-neighbor graph. They often offer closed-form solutions
and have proven successful in real-world applications like classification of digits,
for example. A short overview of such techniques is given next.
ISOMAP [TSL00] finds a projection that preserves the global, nonlinear geometry of the data by preserving the geodesic manifold distances between points. It
first approximates the geodesic distances (by means of a nearest neighbor graph)and then runs multidimensional scaling to find a projection that preserves these
distances.
1

A manifold is an abstract mathematical space in which every point has a neighborhood which
resembles Euclidian space, but in which the global structure may be more complicated.

16

CHAPTER 2. DIMENSIONALITY REDUCTION
Laplacian Eigenmaps [BN03] calculates the embedding map by constructing aweighted graph (based on nearest neighbors or -balls) and calculating the eigenvectors of the graph Laplacian, which is an approximation to the Laplace Beltrami
operator
Lf = −divO(f )The weights are chosen either by the heat kernelkyi − yj k2
wij = exp −tor by setting wij = 1 if an edge exists and wij = 0 otherwise.
Locally Linear Embeddings (LLE) [RS00] assume that the data points Y canbe locally reconstructed in the latent space X. First, weights are calculated byminimizing the cost function
E=

N
X

kyi −

i=1

X

wij yj k2

j∈ηiwhere ηi denotes the set of indices of nodes in the neighborhood of yi . The
P i|invariance constraint |η
j=1 wij ensures that the weights sum up to one. Finallythe latent embedding is found by minimizing
0

E =

N
X

kxi −

i=1

X

wij xj k2

j∈ηiwith respect to x while keeping the weights fixed. It turns out that - by removing
the translational and rotational degree of freedom - this leads to a well posed
problem and can be solved by a standard eigenproblem solver.
Local Tangent Space Alignment (LTSA) [ZZ03] represents the local geometry of the manifold by fitting an affine subspace in a neighborhood of each data
point. The tangent spaces are then aligned to give the internal global coordinatesof the data points with respect to the underlying manifold by way of a partial
eigendecomposition of the neighborhood connection matrix. To do so, a smoothfunction
yi = f (xi ) + 
which generates the data points yi from latent points xi is assumed. Thus a firstorder Taylor expansion
f (x) = f (x) + Jf (x)(x − x) + O(kx − xk2 )is used to reconstruct the data points locally by means of the tangent spacespanned by the columns of the Jacobi matrix Jf (x). Here x is taken as each17

CHAPTER 2. DIMENSIONALITY REDUCTION
local neighborhood’s mean, respectively. [ZZ03] shows that finding the best Qdimensional affine subspace approximation leads to a singular value problem.
Global coordinates are then extracted by minimizing the overall reconstruction
error, leading to an eigenvalue problem.
Maximum Variance Unfolding (MVU) [WS06] assumes that nearby inputs
match the distances between nearby outputs. It places rigid connections between
neighboring nodes by enforcing hard constraints. The optimizing problem is formulated as a quadratic program:
MaximizeXkxi − xj k2

ij

subject to(1) kxi − xj k2 = kyi − yj k2 for all (i, j) with wij = 1
PN(2)
i=1 xi = 0

(2.6)
(2.7)

where wij = 1 if an edge in the neighborhood graph exists and wij = 0 otherwise.
The problem is made tractable by reformulating the optimization as a semidefinite
program (SDP) resulting in a linear program.
All spectral NLDR methods described above rely on a local neighborhood structure where the neighborhood size k is a parameter and the dimensionality Q is
assumed to be known. Predicting a new data point yi often requires to re-run
the algorithm again. Furthermore they do not in general provide a measure of
uncertainty about a given point, i.e. they are non-probabilistic.
In the following section probabilistic Latent Variable Models are discussed with afocus on the Gaussian Process Latent Variable Model, to which the prior proposedin this thesis will be applied in section 3.3.

2.3.2 Latent Variable Models
Latent Variable Models (LVMs), e.g., MDS or Probabilistic PCA [TB99] (a probabilistic extension of PCA and a special case of the Gaussian Process Latent Variable Model), assume that the data Y has been generated by some latent (whichmeans unobserved) variables X that lie on or close to a low-dimensional manifold.
Probabilistic LVMs relate the latent variables to a set of observed variables via aprobabilistic mapping. Since we focus on the GPLVM, the next section introduces
Gaussian Processes first.

18

CHAPTER 2. DIMENSIONALITY REDUCTION(b) Posterior p(f |X, Y)(a) Prior p(f )Figure 2.1: Sampling from a Gaussian Process. (a) shows 3 functions drawn
at random from a GP prior. (b) depicts 3 functions drawn at random
from the posterior, i.e. the prior conditioned on the five noise freeobservations indicated. The gray areas show the uncertainty for eachpoint x respectively. This figure was taken from [RW06].
Gaussian Processes
Gaussian Processes (GPs) [RW06] underwent a revival in the Machine Learning
community over the past ten years. Since the implementation is simple they
became useful tools for probabilistic regression and classification. Here GPs are
introduced in the context of regression and we give a short illustrative example.
A Gaussian Process is a stochastic process and can be seen as an infinite dimensional Gaussian distribution.
Definition 1 A Gaussian Process is a collection of random variables, any finite
number of which have a joint Gaussian distribution. It is completely specified by its
mean function m(x) and its covariance function k(x, x’) and we write a GaussianProcess as
f (x) ∼ GP(m(x), k(x, x’))
wherem(x) = E(f (x))
k(x, x’) = E [(f (x) − m(x))(f (x’) − m(x’))]

(2.8)
(2.9)

The specification of the covariance function implies a distribution over functions.
To see this, one can draw samples from this distribution, as done in figure 2.1.

19

CHAPTER 2. DIMENSIONALITY REDUCTION
Note that the definition of GPs automatically implies a consistency requirement,
also known as the marginalization property, since the covariance function k(x, x’)
specifies entries of the covariance matrix K. The covariance function is also calleda kernel function. This property means that if the GP (e.g.) specifies (y1 , y2 ) ∼N (µ, Σ), then it must also specify y1 ∼ N (µ1 , Σ11 ) and y2 ∼ N (µ2 , Σ22 ).
As for section 2.3.1, let Y = [y1 , · · · , yN ]T be the set of observations yi ∈ <D ,and let X = [x1 , · · · , xN ]T be the set of latent variables xi ∈ <Q , with Q  D.
We further assume a mapping from X to Y as following:y = f (x) + η(2.10)

where y differs from the function values f (x) by additive noise, which is assumed
to be independent, identically and Gaussian distributed:η ∼ N (0, σn2 )By putting a GP prior over the function space, marginalization of f leads to themarginal likelihoodZp(y|X) = p(y|f, X)p(f|X)df(2.11)

where p(f|X) ∼ N (0, K) and f is a vector containing the function values of X.
Since the likelihood p(y|f, X) and the prior p(f|X) are Gaussian, the marginal
likelihood is again Gaussian distributed. More insight can be gained by comparing
the weight space view to the function space view. Consider the decomposition of
f into non-linear basis functions ϕif (x) =∞X

θi ϕi (x)i=1where θi are normally distributed weights θi ∼ N (0, σθ2 ). Marginalizing the weights
in the weight space view corresponds to marginalizing the functions in the function space view and is straightforward, due to the conjugate property of the prior:
Since f (x) is a linear combination of Gaussian distributed weights θi , the resulting
distribution over f is again Gaussian. Considering an instance of the GP for a finite number of training points N , the mean and the covariance are straightforward
to calculate:20CHAPTER 2. DIMENSIONALITY REDUCTION∞Xµk = E!
θi ϕi (xk )

i=1
∞
X

Σkl = E=!

θi ϕi (xk )

i=1
∞
X

= E

∞
X

E (θi ) ϕi (xk ) = 0

i=1
∞
X

(2.12)!!
θi ϕi (xl )(2.13)i=1!

θi ϕi (xk )θi ϕi (xl )(2.14)E (θi ϕi (xk )θi ϕi (xl ))(2.15)

i=1

=

∞
X
i=1

=

σθ2

∞
X

ϕi (xk )ϕi (xl )(2.16)i=1= σθ2 k(xk , xl )(2.17)

Since we are operating in an inner product space, we make use of the kernel
trick (in the last equality): A non-degenerate kernel like the radial basis function(RBF)
kxi − xj k22
k(xi , xj ) = θ1 exp −2θ22!

allows for implicitly making use of an infinite dimensional space (i = 1, ..., ∞)without evaluating it explicitly (which would be intractable). Here the hyperparameters Θ = {θ1 , θ2 } stand for the RBF lengthscale (or kernel height) θ1 and thekernel width θ2 .
Putting it all together we arrive at a simple expression for equation 2.10. If we
restrict the number of training points to be finite, f (x) is Gaussian distributed.
Thus we have a conjugate prior p(f|X) and y is normally distributed, too. This
results inp(y|X) ∼ N (0|K)
with K the covariance matrix defined by a RBF + noise kernel:kxi − xj k22
Kij = k(xi , xj ) = θ1 exp −
2θ22

21

!
+ θ3 δij .

(2.18)

CHAPTER 2. DIMENSIONALITY REDUCTION
Here Θ = {θ1 , θ2 , θ3 } are the hyperparameters, where θ1 is the RBF lengthscale,θ2 the kernel width, and θ3 the noise. The latter not only models the observationnoise, but also contributes to a numerically stable inversion of the covariance K.
The GP key predictive equations provide the mean µ(f∗ ) and the covariance σ 2 (f∗ )for a given test point x∗ given the training points {X, Y}. They are readily derived
from the Gaussian conditional distribution asp(f∗ |x∗ , X, y) ∼ N (µ(f∗ ), σ 2 (f∗ ))
where

−1
µ(f∗ ) = K∗,f Kf,f + σn2 Iy
−1
σ 2 (f∗ ) = K∗,∗ − K∗,f Kf,f + σn2 I
Kf,∗(2.19)(2.20)

Here Kf,f denotes the covariance matrix of the training points, K∗,f is the Kernel
Density Estimate (KDE) with respect to the test point f∗ and K∗,∗ is the varianceof the test point. Note that Gaussian Processes always represent a projection
from <Q to <1 , whereas this limitation is overcome by the GPLVM, introduced
in section 2.3.2.
To illustrate GP regression we might first consider a collection of only 2 randomvariables which gives rise to a 2 dimensional Gaussian distribution as an instanceof the full GP. To illustrate the shape of a Gaussian, figure 2.2 shows face plotsof 2 dimensional Gaussian distributions generated from

N (0,

1 Υ
Υ 1


)

with different off-diagonal elements Υ. Υ = K12 = K21 = 0.9 is assumed for figure
2.3. Given the training point y1 , the mean and variance for y2 can be calculated
easily using equation 2.19 and 2.20.
A slightly more complex example is shown in figure 2.4. 3 input points are given
along with the covariance matrix where Kij = k(xi , xj ) measures the distancealong the abscissa between two points. This figure shows clearly that the uncertainty increases with the distance of the test points from the training points.
Free GP toolkits for visualization of GPs, regression and classification (MATLAB
and C) are available online2 .
2

http://www.gaussianprocess.org/, http://www.rainsoft.de/projects/gausspro.html

22

CHAPTER 2. DIMENSIONALITY REDUCTION

2

1.5

1.5
2

1
y2

y2

1

y2

y2

1

1
0.5

–1.5

–1

–0.5

0.5

0.5

1

1.5

–1.5

–1

–0.5

0.5

y1
–0.5

1

–2

1.5

–1

1

2

–2

–1

1

2

y1

y1

y1
–0.5
–1

–1

–1

–1

–2
–1.5

–1.5
–2

(a)

(b)

(c)

(d)

Figure 2.2: Face plots of 2D Gaussians. Here face plots of 2 dimensional
Gaussians are shown for different parameters of the covariance matrix.
While (a) is uncorrelated (Υ = 0), (b) is correlated with Υ = 0.5. (c)exhibits almost perfect linear coherence (Υ = 0.99) whereas (d) is
negatively correlated (Υ = −0.8).

Figure 2.3: 2D illustration of Gaussian Process regression. Here y1 is given
and the correlation between y1 and y2 is assumed to be known bymeans of a kernel function k(xi , xj ). Since y1 and y2 are correlated,the mean of the conditional distribution p(y2 |y1 ) is similar to the meany1 itself while the small variance indicates only little uncertainty.

23

CHAPTER 2. DIMENSIONALITY REDUCTION

Figure 2.4: 5D illustration of Gaussian Process regression. Here 3 training
points are given compared to only one in figure 2.3. Since we have2 test points, the regression results in a 5 dimensional Gaussian distribution. Note that uncertainty increases with the distance of a test
point from the training data.
Gaussian Process Latent Variable Model
While the training points for GP regression include Y as well as X, the GPLVM
assumes that the latent variables X are unknown and optimizes the model likelihood with respect to X and the hyperparameters Θ in an unsupervised fashion inorder to learn the latent space. Furthermore it allows for a D-dimensional data
space Y by multiplying D single Gaussian Processes, one for each data dimension,
resulting in the final GPLVM likelihood. As it turns out, starting from a goodinitial guess is key to avoid local minima, since the target function is non-convex.
More formally, lety (d) = f (d) (x) + ηwith y (d) the d-th coordinate of y, and η ∼ N (0, θ3 ) iid Gaussian noise. The Gaussian Process Latent Variable Model (GPLVM) [Law05] places a Gaussian process
prior (see section 2.3.2) over the space of mapping functions f . Marginalizing overthe functions f and assuming conditional independence of the output dimensions
given the latent variables results in the GPLVM likelihood

24

CHAPTER 2. DIMENSIONALITY REDUCTION

p(Y|X) =DY(d)p(y |X) =

d=1

D
Y

N (y(d) |0, K)(2.21)

d=1

where y(d) is the d-th column in Y, and K is the covariance matrix, definedin terms of a kernel function. Here we use the RBF + noise kernel, defined inequation 2.18, since it allows for a variety of smooth, non-linear mappings using
only a limited number of hyperparameters (|Θ| = 3).
The log-likelihood of equation 2.21 is given asln p(Y|X) = −D
1DN
ln 2π − ln |K| − trace(K−1 YYT ).
222(2.22)

While the first term is a constant, the second term can be interpreted as a regularizer, which prevents K from getting too big (for example by means of highnoise or a large kernel height) and the third term is the data-fitting term, whichtries to maintain the relationship between points in the input space.
Learning the GPLVM is performed by maximizing the posteriorp(X|Y) ∝ p(Y|X)p(X)with respect to the latent variables X and the kernel hyperparameters Θ. Herep(X) encodes prior knowledge about the latent space X. In [Law05] the redundancy between the overall scale of the matrix X and the kernel
QN width is removed
by chosing a product of isotropic Gaussian priors p(X) = i=1 N (xi |0, I) whichleads to a Wishart distributed scatter matrix XT X. In practice however, thiscan lead to degenerate solutions (like X = 0). To prevent from this behaviour,
we applied constrained optimization instead and used another prior, which penalizes the dimension, as described in section 3.3. While [Law05] uses SCG (Scaled
Conjugate Gradients) [Mol93] we made use of SNOPT (Software for Large-Scale
Non-Linear Programming), a constrained non-linear optimizer [GMS02] based onsequential quadratic programming.
PCA and graph-based techniques are commonly used to initialize the latent spacein GPLVM-based dimensionality reduction; both offer closed-form solutions. However, PCA [Pea01] cannot capture non-linear dependencies, LLE [RS00] gives agood initialization only if the data points are uniformly sampled in the manifold,and ISOMAP [TSL00] has difficulties with non-convex datasets [Har07]. Generally, when initialized far from the true minimum, the GPLVM optimization can
get stuck in local minima [Law05, UFG+ 08].
To avoid this problem, different priors over the latent space have been developed(see chapter 3).

25

CHAPTER 2. DIMENSIONALITY REDUCTION
A further issue evolves from the fact that GPLVMs rely on the calculation of theinverse of the covariance matrix, which prohibits the use of a large number oftraining points and will be addressed in the next section.

2.3.3 Sparsification
One of the main limitations of using Gaussian Processes in practical applications is
the computational complexity of learning a model - like for many kernel methods ingeneral. Since learning involves the inversion of a N × N -matrix, computational
cost3 is approximately O(N 3 ). This restricts applicability to cases with N <5000.
While several methods for Gaussian Process sparsification have been developed[SNS06, YDD05, UD08, SW02, QCR05], it is a bit more subtle to overcome thetime limitiations for GPLVMs. This is due to the fact that in each iteration the
inverse of the covariance matrix K−1 is needed for calculating the gradient ofequation 2.22 with respect to X. One way of increasing the speed is to simply
remove points from the training set by subsampling or applying information criterions like the one proposed in the Informative Vector Machine (IVM) [LSH02].
If further speed-up is needed without rejecting training points from the dataset,
one has to resort to more sophisticated techniques, like the ones described in thefollowing.
Sparse Gaussian Process Approximations
A unifying view, which includes many existing probabilistic sparse approximations
for Gaussian Process regression has been developed in [QCR05]. The approach relies on expressing several methods in terms of the so-called effective prior. Rather
than the interpretation ”approximate inference with exact prior” the methods arereinterpreted as ”exact inference with an approximated prior”.
Let us first consider inference for GP regression. Sincep(f, f∗ |y) =p(f, f∗ )p(y|f∗ , f)
p(y)

and

p(y|f∗ , f) = p(y|f)

inference can be done by marginalizing out the unwanted training variablesZp(f∗ |y) =3

1
p(f, f∗ |y)df =p(y)Zp(y|f)p(f, f∗ )df(2.23)The Strassen algorithm can solve the matrix inversion problem in O(N 2.807 ). The fastest
known algorithm - derived by Coppersmith and Winograd - lies in O(N 2.376 ), but has hugeconstants hidden in the O notation and is therefore not used in practice.

26

CHAPTER 2. DIMENSIONALITY REDUCTION
where p(y|f) = N (f, ση2 I) is the likelihood, p(f, f∗ ) is the conjugate joint prior and1
is a normalization constant. The predictive distribution p(f∗ |y) is Gaussianp(y)with the first and second central moment given by equation 2.19 and equation
2.20.
In order to reduce computational requirements from 2.23, several approximations
to the joint prior p(f∗ , f) are reviewed in [QCR05]. The prior has been rewritten
by introducing another set of M latent variables u = [u1 , ..., uM ]T which are called
inducing variables. They are values of the Gaussian process, corresponding to aset of input locations Xu , which are called inducing inputs. Since GPs exhibit
the marginalization property, the joint GP prior p(f∗ , f) can be rewritten as themarginalZp(f∗ , f) =Z
p(f∗ , f, u)du =p(f∗ , f|u)p(u)duwhere p(u) = N (0, Ku,u ). The approximation is then done by assuming that f∗
and f are conditionally independent given u, such thatZp(f∗ , f) 'p(f∗ |u)p(f|u)p(u)du.

The different computationally efficient algorithms proposed in the literature correspond to different additional assumptions about the two approximate inducing
conditionals p(f∗ |u) and p(f|u). Applying sparse Gaussian Process approximation
reduces the computational complexity from O(N 3 ) to O(N M 2 ) where N is thenumber of training points and M is the number of inducing variables. Further details as well as the individual approximation schemes are discussed in [QCR05].
In [Law07] the method reviewed in [QCR05] was transferred to the Gaussian Process Latent Variable Model and the optimization was performed with respect toX, Xu and Θ. While being an elegant way to reduce computational complexity,
introducing additional variables (Xu ) into the model also makes the optimization
prone to overfitting and introduces local minima. In practice, the inducing variables also tend to be placed at non-intuitive locations after optimization which
has to be further investigated.

Nyström Method
The Nyström Method [WS01] computes an approximation to the eigendecomposition of the Gram matrix K. It is done by carrying out an eigendecomposition
on a smaller system of size M < N and then expanding the result back to N
dimensions. The computational complexity therefore reduces to O(N M 2 ).

27

CHAPTER 2. DIMENSIONALITY REDUCTION
Considering the kernel function

k(x, y) =RXλi φi (x)φi (y)i=1where R ≤ ∞ and λi are the eigenvalues to the eigenfunctions φi , such thatZk(x, y)φi (x)p(x)dx = λi φi (y)where p(x) denotes the probability density of the input vector x. The integral overp(x) is replaced by an empirical average (which is called the Nyström Method ) toobtainq

1Xk(y, xk )φi (xk ) ≈ λi φi (y).
q k=1This motivates the matrix eigenproblem
KU = UΛ.
Computational costs of computing K can now be cut down by using the eigendecomposition K = UΛUT where U is orthonormal and Λ a diagonal matrix
with the eigenvalues of K. Williams et al approximate K + σI by UΛUT + σI.
Furthermore they show that the eigenvalues and eigenvectors which are needed
can be approximated by

(N)

λi

(N)

ui

N (M)
λ
M i
r
M 1
(M)
≈
KN,M uiN λ(M)i≈(2.24)(2.25)

n
oN
(N)
(N)
where (λi , ui )are eigenvalues and eigenvectors of the full kernel matrix
i=1
n
o
M
(M)
(M)
K and (λi , ui )are eigenvalues and eigenvectors of a M × M submatrix
i=1
of K.
If the eigenvalues of K decay rapidly, good approximations can be obtained byusing M  N . However there is no clear answer to the choice of the subset
of points which give rise to the M × M submatrix of the N × N matrix K.

28

CHAPTER 2. DIMENSIONALITY REDUCTION
Furthermore, using our datasets, it turned out that the eigenvalues of K didn’t
decay as quickly as needed for using the Nyström method, mainly because of alarge kernel width.
Therefore we employed the simple method described next to speed up learning oflarger motion datasets.
Product of local and global Gaussian Processes
The method proposed here is simple to implement and based on the assumptionthat topological knowledge about the latent points is given. More precisely, it issufficient to know local neighborhoods. Thus this method can be applied when
dealing with motion capture data, for example.
Instead of mapping the full dataset from the latent space X to the data space
Y, different subsets of X are mapped to different subsets of Y. This results in a
product of products of Gaussian Processes, or a product of GPLVMs (see equation
2.21). The likelihood of the sparse GP can be written asp(Y|X) ≈DY Y

N (y(d)s |0, Ks )(2.26)

s∈L∪G d=1

where L is a set of indices capturing overlapping local neighborhoods along the
motion trajectory and G is a set of indices which puts points into global relation
to each other. The idea is illustrated in figure 2.5. The resulting computational
complexity is O(M 3 ) with M = maxs∈L∪G |s|. Note that this is not a proper
sparsification technique in a probabilistic sense, but rather serves as a reasonablead-hoc approximation to the gradients needed for the GPLVM.

29

CHAPTER 2. DIMENSIONALITY REDUCTION

Figure 2.5: Sparsification: Global and local GPs. This figure illustrates thesparsification method used for motion datasets with a large numberof training points. The colored points are training points, lying on a
motion trajectory (black) in the data space. Instead of inverting thefull Gaussian Process covariance matrix (which takes O(N 3 ), with N
the number of training points), we make use of prior knowledge about
the nature of motion data, i.e. we assume smoothness over time (orframes). Thus only a set of smaller local GPs (which overlap withneighboring sets, shown in gray) is calculated, reducing the computation time dramatically. To also keep track of global relationships, eachpoint in a local GP set is linked to one point in each other local GP set
via a global GP (which thus again consists of a very limited numberof training points, contributing to a speed-up in learning the model).
In this example 3 global GPs are depicted (red, green and blue).

30

3 Priors for Latent Variable Model
To overcome the problem of local minima and to learn smooth latent spaces forhuman motion, different priors over the latent space have been developed. In[WFH08] a prior was introduced in the form of a Gaussian process over the dynamics in the latent space (see section 3.1). This results in smoother manifolds
but performs poorly when learning stylistic variations of a motion or multiple motions. Urtasun et al. [UFG+ 08] proposed a prior over the latent space, inspired by
the LLE cost function, that encourages smoothness and allows the introductionof prior knowledge, e.g., topological information about the manifold (see section
3.2). However, such prior knowledge is not commonly available, reducing considerably the applicability of their technique. In contrast, the method developed in
section 3.3 below introduces a generic prior that requires no specific prior knowledge, directly penalizing the dimensionality of the latent space to learn effective
low-dimensional representations.

3.1 Dynamical Priors
Dynamical Priors, as introduced in the Gaussian Process Dynamical Model (GPDM)[WFH08], assume knowledge about the dynamics in the latent space, i.e. that data
points have been recorded consecutively with a constant frame rate. They thus
employ a model where an additional mapping within the latent space captures afirst-order Markov dynamics
x(t+1) = g(x(t) ) + η (t)xy(t) = f (x(t) ) + η (t)
y

(3.1)
(3.2)

where t represents a discrete time index. The prior p(X) is then given as
Z Y
N
(1)
p(X) = p(x )p(x(t) |x(t−1) , g)p(g)dg.
t=2

In [WFH08] a GP prior and a linear + RBF kernel was used for modeling the
dynamics mapping. They show that incorporating the dynamics prior into the
GPLVM likelihood results in smoother latent trajectories which can be used forbetter motion generation or articulated body tracking. However, when the training
data contains large stylistic variations and multiple motions, the generic GPDM
and the back-constrained GPDM [UFG+ 08] do not produce useful models.

31

CHAPTER 3. PRIORS FOR LATENT VARIABLE MODEL

3.2 Topological Priors
In [UFG+ 08] we introduced Topologically-Constrained Latent Variable Models which
incorporate a prior based on knowledge about walk cycle phases to align motions
and separate style (e.g. weak vs. exaggerated arm movement) from content (thewalking cycle itself). Furthermore when switching between styles (e.g. walking
vs. running) a transition is restricted to be possible only when at least one foottouches the ground. The application is graphics animation, where an animator
specifies foot constraints, for example, and a motion is generated from a specificmodel with respect to the specified constraints.
In order to do so we inspired ourselves by the Locally Linear Embedding (LLE)
NLDR method (see section 2.3.1). Our proposed prior is given by
!
N
D
X (j)
1 XX
kXij −
wik Xij k2
p(X|W) ∝ exp − 2
σ i=1 j=1k∈ηiwhere W is the weight matrix, σ 2 represents a global scaling of the prior and ηiis the local neighborhood of data point i. In LLE the weights are computed bysolving ∀j∀k ∈ ηi
X
(j)
Clk wik = 1lsim
where Clk
= (yi − yl )T (yi − yk ) if l, k ∈ ηi , and 0 otherwise.

To incorporate transitions, the covariance matrix C can be modified as
Cijtrans = 1 − δij exp −ζ(ti − tj )2
 Mwith ζ a constant, and δij = 1 if ti and tj are in the same set t̂k k=1 and otherwise M
δij = 0. Here t̂k k=1 is defined as a subset of frames where transitions are
possible. This covariance encourages the points at which transitions are physicallypossible to be close together in the latent space.
Cyclic topologies can be enforced by introducing covariances based on the phase,
specifying different covariances for each latent dimension:cos
Cjk= (cos(φi ) − cos(φj ))(cos(φi ) − cos(φk ))

(3.3)

sin
Cjk

(3.4)

= (sin(φi ) − sin(φj ))(sin(φi ) − sin(φk ))We also show, that it is possible to use the back-constrained GPLVM [LQC06]
to incorporate prior knowledge by modifying the kernel function. Experiments
show, that the learned models are smooth and exhibit desirable properties for32

CHAPTER 3. PRIORS FOR LATENT VARIABLE MODEL
character animation, for example. However, the prior knowledge needed is not always available. While the cyclic assumption is reasonable for walking or runningmotions, motions in general do not exhibit this structure, thus violating the assumption (consider for example arbitrary motions performed in a kitchen scenario
as discussed in section 4.3 of this thesis, see figure 4.14). The Rank Prior introduced in the next section seeks for a viable embedding without making additional
assumptions about the dynamics or the topology of the dataset.

3.3 Rank Priors
Here we introduce a novel method for probabilistic non-linear manifold learningwhich avoids the initial distortion induced by ad-hoc initialization. The idea is
illustrated in the figures 3.1 and 3.2. While figure 3.1 describes the general idea
motivated by a spring-force system, figure 3.2 shows the actual algorithm whenapplied to a spiral manifold embedded in two dimensions.
In our method we initialize the latent space to the high-dimensional observation
space and define a Rank Prior which favors latent spaces with low dimensionality.
Minimizing the dimensionality of the latent space is equivalent to minimizing R,
the rank of the unbiased covariance estimate (along the dimensions) of X̃
R = rank

1
T
X̃ X̃
N −1


(3.5)

where X̃ denotes the mean subtracted latent variables X̃ = X − X and Xij =PN
1
n=1 Xnj .
N

3.3.1 Rank Prior Formulation
The cost function in Eq. (3.5) is discrete and thus difficult to minimize. However,it can be expressed in terms of singular values and singular vectors of the covariance matrix along the dimensions. Let the singular value decomposition (SVD) ofX̃ be UΣVT , where U and V are matrices containing the left and right singular
vectors, and Σ comprises the singular values {σ1 (X̃), ..., σD (X̃)} on its diagonal.
T
Then X̃ X̃ = VΣ2 VT , and R can be minimized by minimizing the number of nonzero singular values of X̃. This connection can also be drawn
 by comparing
 T theTnull spaces of X̃ and X̃ X̃, leading to the insight that rank X̃ = rank X̃ X̃ .
Reducing the dimensionality of X̃ is equivalent to reducing the rank of X̃ or thenumber of non-zero singular values.

33

CHAPTER 3. PRIORS FOR LATENT VARIABLE MODEL

Figure 3.1: Illustration using a spring-force system. From a simplified viewpoint the method can be seen as a system of observation points (red)in the high-dimensional data space, connected by springs (green andblue) of strengths defined according to the distance of the two points
they are connnecting, respectively. Thus the spring strength can be
defined by a radial basis function, for example. Here, e.g., the blue
spring is weaker than the two green springs since the data points it
connects are further apart from each other. Before applying the algorithm, the springs are in their home position, imposing no force.
The spring system stands for the GPLVM, which preserves dissimilarities. On the other hand, the pressure plates surrounding the systemapply forces towards the covariance directions. They represent the
Rank Prior which favors low-dimensional spaces. Here the pressure
along the smaller variance is strongest since we want to get rid ofthe ”smaller” dimensions first. This is achieved by using a non-linear
penalty function, as described in section 3.3.3.

34

CHAPTER 3. PRIORS FOR LATENT VARIABLE MODEL

Figure 3.2: Illustration of our Rank Prior with a GPLVM. The goal is torecover the 1D manifold embedded in 2D. The GPLVM gets stuck in
local minima very early (upper row) since PCA initialization does not
capture non-linear dependencies, whereas our method decreases dimensionality gradually and recovers the correct manifold (lower row).
We transform the discrete optimization criteria in Eq. (3.5) into a continuous one
by introducing a sparsity penalty on the singular values. In particular we introduce
a prior of the form!
DX
1
ϕ(σ̂m (X̃))(3.6)
p(X) = exp −Z
m=1
where σ̂m (X̃) = √N1−1 σm (X̃) is the normalized singular value of X̃, ϕ is a spar
 P
R
ϕ(σ̂
(
X̃))
dX̃ a normalizationsity penalty function, and Z = exp − Dmm=1
constant1 .
One might consider a variety of sparsity penalty functions. The identity function
ϕ(σ̂) = σ̂ results in the L1 -norm since the singular values are always positive.
Of particular interest to us are functions that drive small singular values fastertowards 0 than larger ones. Examples of such functions are the logarithmic ϕ(σ̂) =α ln(1 + β σ̂ 2 ) and the sigmoid ϕ(σ̂) = α(1 + exp(−β(σ̂ − γ)))−1 functions, with α,β and γ constant parameters.
Minimizing the negative log posterior − ln p(Y|X) − ln p(X) results in an optimization that reduces the dimensionality in a continuous fashion:
"
#
D
X
D
D
min
ln |K(X, Θ)| + tr(K(X, Θ)−1 YYT ) +ϕ(σ̂m (X̃))
s.t. ∆E = 0
X
2
2
m=1
(3.7)
1

Please note that the fact that this is an improper prior has no impact in the optimization
since it acts as a constant when minimizing the negative log posterior.

35

CHAPTER 3. PRIORS FOR LATENT VARIABLE MODEL
where∆E = |E(X̃) − E(Ỹ)|
is the difference between the energies of the spectrum of the mean subtracted
latent coordinates X̃ and the mean subtracted observations, Ỹ = Y − Ȳ, withȲ the mean of the observations. This constraint keeps the overall energy of thesystem constant during optimization and thus removes the redundancy between
the overall scale and the kernel width. The energy of the spectrum at the current
iteration is defined as
D
X
2
(X̃)
E(X̃) =
σ̂mm=1while the energy of the initial spectrum equals the energy of the data points’
spectrum:DX
2
E(Ỹ) =
σ̂m
(Ỹ).
m=1

Equation 3.7 is straightforward to optimize, the required derivatives are given
below.

3.3.2 Derivatives
The derivative of the GPLVM log likelihood (equation 2.22) with respect to thekernel matrix K is given by∂ ln p(Y|X)
1
D= K−1 YYT K−1 − K−1∂K22
and can be composed with∂k(xi , xj )∂xi
∂k(xi , xj )∂θ1∂k(xi , xj )∂θ2∂k(xi , xj )∂θ3
kxi − xj k2
θ1= − 2 (xi − xj ) exp −
θ2
2θ22


kxi − xj k2= exp −
2θ22
θ1 kxi − xj k2
kxi − xj k2
=
exp −
θ23
2θ22
= δij

(3.8)
(3.9)
(3.10)(3.11)

via the chain rule36

CHAPTER 3. PRIORS FOR LATENT VARIABLE MODEL∂f (K)= trace∂Xij"∂f (K)
∂K

T

∂K
∂Xij

#

The derivative of the rank prior is also easy to calculate. Since X̃ = UΣVT wehave UT X̃V = Σ. Thus, since Σ is diagonal, the derivative of the rank prior withrespect to the latent coordinates follows as
D
D
X
∂ X
1
∂ϕ(σ̂m )ϕ(σ̂m ) = √Uim Vjm ,∂Xij m=1
N − 1 m=1 ∂ σ̂m(3.12)

m)depends on the sparsity function and Uim and Vjm are elements of
where ∂ϕ(σ̂∂ σ̂mthe respective singular vectors.

We use the SNOPT [GMS02] non-linear constraint optimizer to minimize Eq.
(3.7). After optimization, we select the latent dimension as
Q = argmaxm

σ̂mσ̂m+1 + 

where   1, and σ̂1 ≥ σ̂2 · · · ≥ σ̂D . Alternatively one could also threshold thesingular value spectrum.

3.3.3 Putting it all together: the Algorithm
The final steps consist of applying PCA in the optimized Q-dimensional spaceand optimizing p(Y|X, Θ) with respect to the hyperparameters Θ. Note that the
mapping is still non-linear since PCA is performed in the latent space, not in
the observation space, and simply rotates the data to produce the most compact
Q-dimensional representation.
An overview of the complete algorithm is given in the following.
Algorithm 1 Continuous Non-Linear Dimensionality Reduction
1: Normalize Y and set X = Y
2: Select parameters (Θ, α, β, γ)3: Optimize equation 3.7 with respect to X
4: Apply PCA to X for final rotation (this reduces Q)
5: Optimize equation 2.22 with respect to Θ while keeping X fixed
Figure 3.2 compares the GPLVM (initialized with PCA) with the result of optimizing Eq. (3.7) on a toy example where a 1D manifold is embedded in 2D space.
PCA provides a non-optimal initialization, and the GPLVM gets trapped in local
minima whereas our method recovers the correct structure. Note that our final

37

CHAPTER 3. PRIORS FOR LATENT VARIABLE MODEL
PCA projection rotates the latent space and results in a 1D manifold. In thisexample, using spectral methods could lead to a successful initialization for the
GPLVM. However, for more complex datasets this is not necessarily the case ingeneral, as shown in the examples section in Figs. 4.1 and 4.15.
Figure 3.3 (a),(b) and (c) depicts the evolution of the first ten singular values
when optimizing Eq. (3.7) with a linear, sigmoid and logarithmic sparsity penalty
function respectively (for a motion database composed of 30D observations). Note
how our method drops dimensions as the optimization evolves (i.e., the smallest
singular values drop to zero within the first few iterations). A comparison of thespectrum of different sparsity penalty functions is shown in figure 3.3 (d). The L1 norm results in a poor estimation of the dimensionality, while the more aggressive
sigmoid and logarithmic functions are able to recover the correct dimensionality inthis example. In the remainder of the paper we use the logarithmic function since
it converges faster than the L1 -norm and has fewer parameters than the sigmoid.
Since the final result often depends strongly on the sparsity penalty function andits parameters further investigation will be needed concerning this choice.

38

2

2

1.5

1.5

1

1

σ

σ

CHAPTER 3. PRIORS FOR LATENT VARIABLE MODEL

0.5

0.5

0
0

20

40
60
iteration number

80

0
0

100

(a) Linear penalty function

20

40
60
iteration number

80

100

(b) Sigmoid penalty function

2.5

2observation space
linear sparsity function
logarithmic sparsity function
sigmoid sparsity function

2
1.5
σ

σm

1.5
1

1
0.5

0.5
0
0

20

40
60
iteration number

80

0

100

2

4

6

8

10

m

(c) Logarithmic penalty function(d) Comparing the resulting spectrum

Figure 3.3: Spectrum of a 30D motion database. (a),(b) and (c): Evolutionof the first ten singular values as a function of the optimization iteration number for a linear, sigmoid and logarithmic sparsity penalty
functions. (d): Spectrums learned after convergence by different sparsity penalty functions compared to the observation space spectrum(red).

39

CHAPTER 3. PRIORS FOR LATENT VARIABLE MODEL

40

4 Experimental results
In this section we demonstrate our approach in three different scenarios. In section 4.1 we compare our method to graph-based techniques and the GPLVM withdifferent initializations in artificial data (a sparsely sampled swiss roll). We further illustrate our method’s ability to estimate the latent space dimensionality incomplex synthetic data in section 4.2. Finally we present an application of ourtechnique to the challenging problem of tracking and classifying 3D articulated
human body motion in the kitchen domain (see section 4.3).

4.1 A sparsely sampled swiss roll
The swiss roll is a widely used example of a 2D manifold which is embedded ina 3D space. Since the manifold is embedded in the data space in a non-linear
fashion, PCA [Pea01] fails to recover it. However, many state-of-the-art graphbased techniques (like ISOMAP [TSL00], Laplacian Eigenmaps [BN03], LocallyLinear Embedding [RS00], Local Tangent Space Alignment [ZZ03] and Maximum
Variance Unfolding [WS06]), which rely on local neighborhoods and have a closedform solution, can be used to unravel it correctly if the data is homogenously
sampled, the noise is small, and the neighborhood size is selected appropriately.
However, real data often violates these assumptions resulting in poor performance.
Imagine for example a motion dataset: Using a high-frequency motion capturingsystem, many training points can be obtained from a motion along its trajectoryin the data space, but the sampling is very sparse in the direction of different
styles. This explains the need for more sophisticated methods than the spectral
methods mentioned above when it comes to handle complex and sparse datasets
like human motion.
Therefore we have created an example where only a very sparse and noisy subsetof the swiss roll is assumed to be known to the algorithm (see figure 4.1a, points
in black). The goal is to recover the correct underlying manifold (depicted on theleft side). We then compare the reconstruction of the test points in 2D and theirmean reprojection in 3D to the ground truth (figure 4.1a).
The first row in figure 4.1 (b) shows the result of applying PCA, Isomap, Laplacian Eigenmaps, LLE, LTSA and MVU (see [vdMPvdH07] for a review on thesetechniques). The second row depicts our technique and the result of optimizing

41CHAPTER 4. EXPERIMENTAL RESULTS

Figure 4.1: Finding a 2D manifold in 3D space on a sparsely sampledswiss roll. Only a sparse, noisy subset (depicted in black) of the fullmanifold is assumed to be known (a). (b) shows the initialization (withneighborhood size k=6), GPLVM result and 2D/3D reconstruction of
the full manifold (from top to bottom).

42

CHAPTER 4. EXPERIMENTAL RESULTS
the GPLVM with different initializations. Finally the last two rows of figure 4.1(b) show the test data (i.e., colored samples) reconstructed in the latent space andin the original space.
The results of our method are shown in the last column of figure 4.1b. To visualize
our algorithm, figure 4.2 shows how the latent space (which in our case is initialized
to the data space) evolves over time.
Note that our method, unlike PCA, graph-based techniques and the GPLVM withany of the initializations, is able to recover the correct manifold. The bending
effect at the boundaries of the manifold (which causes the manifold to spread at
the red and yellow end) is explained by the fact that the GPLVM preserves global
distances (the inner points are close to the outer points) rather than only local
relationships between data points next to each other.
We evaluate the performance of the different algorithms on this example computing a global and a local measure of accuracy. The reconstruction error is a
global measure of the ability to generalize, and was obtained by first finding thelatent coordinates x∗ of the test data y∗ by maximizing p(x∗ |y∗ , X, Y), and thencomputing the average mean prediction error
1 X
kµ(x∗i ) − y∗i k2
Nt iwith Nt the number of test data. The relationship error, Rerror , measures how
well local neighborhoods are preserved and is defined as
Rerror =Nt X
X

Γij − Γ̄ij

2

i=1 j∈ηikx −x kwhere ηi is the set of neighbors of the i-th test data, Γij = kyii −yjj k22 is the ratio
between the distance in the latent space and the distance in the observation space
for two neighbors, and Γ̄ij is the mean ratio in the local neighborhood. figure 4.3
and figure 4.4 depict these two error measures when performing the experimentin figure 4.1 averaged over 20 random partitions of the data. We use a local
neighborhood of size 4 to compute the relationship error in all experiments, anda logarithmic sparsity function with α = 10, β = 10, and Θ = {0.5, 1.5, 0.01} .
The hyperparameters were optimized for the GPLVM baselines. Note that our
method outperforms the baselines independent of the initialization used for the
GPLVM.

43

CHAPTER 4. EXPERIMENTAL RESULTS(a) Iteration 0(b) Iteration 10(c) Iteration 29(d) Iteration 48(e) Iteration 60(f) Iteration 91

Figure 4.2: Learning a sparsely sampled Swiss Roll. Here snapshots of thelearning process are shown at different iterations during optimization
using our method. One can clearly see how the roll is unraveled correctly. The final manifold has 2 dimensions.

44CHAPTER 4. EXPERIMENTAL RESULTS(a) Reconstruction Error(b) Relationship Error

Figure 4.3: Example: A sparsely sampled swiss roll. Reconstruction Error(a) and Relationship Error (b) for the experiment in figure 4.1 averaged over 20 random partitions of the data for the 8 best performing
dimensionality reduction techniques. A more detailed result is given
in figure 4.4.

45

CHAPTER 4. EXPERIMENTAL RESULTS

our method
PCA init
ISOMAP init (k=6)
ISOMAP init (k=9)
ISOMAP init (k=12)
ISOMAP init (k=15)
Laplacian init (k=6)
Laplacian init (k=9)
Laplacian init (k=12)
Laplacian init (k=15)
LLE init (k=6)
LLE init (k=9)
LLE init (k=12)
LLE init (k=15)
LTSA init (k=6)
LTSA init (k=9)
LTSA init (k=12)
LTSA init (k=15)
MVU init (k=6)
MVU init (k=9)
MVU init (k=12)
MVU init (k=15)

Reconstruction Errormeanstddev
0.00410.01070.08450.18600.38130.17940.17200.84530.35830.42210.23500.23263.70272.65532.31361.31501.60380.28780.55612.03213.55142.18712.61752.51411.6235
1.13862.73004.54882.63770.72733.51493.88352.17343.20994.09503.66810.37830.5381
0.33830.36650.32280.34770.07060.2560Relationship Errormeanstddev
0.00080.00130.42320.19150.00380.01200.00110.00760.01300.00350.07310.03149.27352.55400.6426
8.1716
7.97841.50292.24000.05087.55911.48078.54851.48810.10580.1064
1.78208.92151.26361.64982.05752.40382.37003.07182.18194.51500.12380.00550.0491
0.04650.06720.01640.58562.8057

Figure 4.4: Example: A sparsely sampled swiss roll. Reconstruction and
Relationship Error for the experiment in figure 4.1 averaged over 20
random partitions of the data. Here more detailed results are shown,including PCA and graph-based methods with different neighborhood
sizes. The 8 best performing techniques are depicted in 4.3.

46

CHAPTER 4. EXPERIMENTAL RESULTS

4.2 Discovering the correct dimensionality
Here we illustrate our method’s ability to discover the intrinsic dimensionalityof the underlying manifold in 5 complex synthetic examples, which all are embedded in 3D space (D = 3), but exhibit a different intrinsic dimensionality(Q = {1, 2, 3}). In figure 4.5a a spiral with a wide separation between rings isreduced to a 1D manifold. When the distance between the different rings decreases, the intrinsic manifold dimensionality changes from 1D to 2D (see figure4.5b), since relationships between points that have the same phase are considered.
In figure 4.5c the underlying 2D manifold from a cut-off sphere sampled along
longitudinal lines is discovered. The manifold in figure 4.5d is intrinsically 3Dand thus cannot be reduced (as discovered by our method), while the manifold infigure 4.5e can be reduced to 2D. The learning process is illustrated in figure 4.6,
4.7, 4.8, 4.9, 4.10, 4.11 and 4.12.

Figure 4.5: Dimensionality estimation. (Top) Five 2D manifolds embedded in
3D. (Bottom) Latent spaces and intrinsic dimensionalities Q learned
using our continuous dimensionality reduction method.

47

CHAPTER 4. EXPERIMENTAL RESULTS(a) Iteration 0(b) Iteration 19(c) Iteration 23(d) Iteration 27(e) Iteration 33(f) Iteration 50

Figure 4.6: Learning a hat shaped manifold generated by sine functions.
Snapshots of the learning process are shown at different iterations
during optimization using our method. The discovered dimensionality
Q is 2.

48

CHAPTER 4. EXPERIMENTAL RESULTS(a) Iteration 0(b) Iteration 10(c) Iteration 13(d) Iteration 16(e) Iteration 19(f) Iteration 50

Figure 4.7: Learning a hat shaped manifold with a hole. Snapshots of thelearning process are shown at different iterations during optimization
using our method. The discovered dimensionality Q is 2.

49

CHAPTER 4. EXPERIMENTAL RESULTS(a) Iteration 0(b) Iteration 21(c) Iteration 30(d) Iteration 50(e) Iteration 65(f) Iteration 100

Figure 4.8: Learning a cut-off bowl manifold sampled along longitudal
lines. Snapshots of the learning process are shown at different iterations during optimization using our method. The discovered dimensionality Q is 2.

50

CHAPTER 4. EXPERIMENTAL RESULTS(a) Iteration 0(b) Iteration 6(c) Iteration 20(d) Iteration 25(e) Iteration 31(f) Iteration 60

Figure 4.9: Learning a 2D wave shaped manifold. Snapshots of the learning
process are shown at different iterations during optimization using our
method. The discovered dimensionality Q is 2.

51

CHAPTER 4. EXPERIMENTAL RESULTS(a) Iteration 0(b) Iteration 11(c) Iteration 16(d) Iteration 30(e) Iteration 50(f) Iteration 100

Figure 4.10: Learning a spiral manifold embedded in 3D. Snapshots of thelearning process are shown at different iterations during optimization
using our method. The discovered dimensionality Q is 1.

52

CHAPTER 4. EXPERIMENTAL RESULTS(a) Iteration 0(b) Iteration 14(c) Iteration 21(d) Iteration 35(e) Iteration 47(f) Iteration 100

Figure 4.11: Learning a spiral manifold embedded in 3D. Snapshots of thelearning process are shown at different iterations during optimization
using our method. The discovered dimensionality Q is 1.

53

CHAPTER 4. EXPERIMENTAL RESULTS(a) Iteration 5(b) Iteration 10(c) Iteration 14(d) Iteration 18(e) Iteration 24(f) Iteration 50

Figure 4.12: Learning a narrow spiral manifold embedded in 3D. Snapshots of the learning process are shown at different iterations during
optimization using our method. The discovered dimensionality Q is2.

54

CHAPTER 4. EXPERIMENTAL RESULTS

4.3 Tracking and classification in the kitchendomainAn interesting real-world application of discovering low-dimensional structure ina high-dimensional space is tracking and classifying human motion from video
sequences. Tracking consists of inferring the 3D locations of body joints fromimages. Since tracking requires good generalization abilities and latent trajectories
to be smooth this is a difficult task.
The goal of this experiment is to evaluate a model learned from different motions
by tracking and classifying human movements performed in a kitchen scenario.
This is done by learning the model from previously captured motion capturing data
and tracking the person’s movements using a particle filter approach (see figure4.17). Once tracking can be performed robustly, classification is straightforward,
e.g. by assigning the class label of the nearest neighbor in the latent space to thecurrent frame.

rolling
millingbroomingØlength
18 s
15 s
11 s

#total6
42#learning2
2
2#tracking
4
2
0

Figure 4.13: Dataset used in the tracking experiment. Here an overview ofthe dataset used in the experiment is given. First, the average length
per sequence is given in seconds. Next the total number of utilizable
sequences, the number of sequences used for learning the model andthe number of sequences used for tracking is shown.
We show how learning multiple motions in a single latent space can be doneusing our method since this is a difficult task in general. We use the kitchen
dataset of [KPFW08], that consists of images from 2D cameras and ground truth
joint angles of 3 motion types (i.e., rolling, milling, brooming), performed several
times by the same subject. Unfortunately, only a subset of the entire dataset was
found to be not corrupted by mislabelings and could be used for the experiments.
More specifically those were 6 rolling motions, 4 milling motions and 2 brooming
motions. For our experiments we learned a common latent space using the firsttwo motions of each type and tested our algorithm on the remaining motions.
Therefore no tests on the brooming motions were performed although they were
included in the training set. The motion capture data was recorded using a VICON
system (at 100 Hz) while the video streams were captured by four Sony VWF 500
color cameras (at 30 Hz). An overview of the available sequences is given in table
4.13 and in figure 4.14.

55

CHAPTER 4. EXPERIMENTAL RESULTS(a) Camera 1: Rolling(b) Camera 2: Rolling(c) Camera 1: Brooming(d) Camera 2: Brooming(e) Camera 1: Milling(f) Camera 2: Milling

Figure 4.14: Tracking sequences. For the tracking and classification experiments the database from [KPFW08] was used. It contains several
types of motion, captured with a VICON motion capturing systemand four Sony VWF 500 color cameras (resolution 640x480). In thisfigure snapshots of 3 motions (rolling, brooming and milling) are depicted from 2 views each.

56

CHAPTER 4. EXPERIMENTAL RESULTS

Figure 4.15: Learning different types of motion into one single 3D latent
space. Each type is visualized with a different color (red = rolling,green = brooming, blue = milling). (a) depicts a 3D-GPLVM initialized to LLE, discontinuities are emphasized in black. (b,c) show the
manifold learned by our method along with the variance.

4.3.1 Learning
We learned a single 3D latent space from 30D joint angle observations of 2 trials
for each of the 3 different motions. Since learning involves the inversion of thecovariance matrix K and thus scales as O(N 3 ) where N represents the number oftraining examples, only a limited number of samples can be used for training themodel.
Due to variation in dynamics, learning a latent space using a simple subsampling(like 1 out of 4) of the dataset leads to suboptimal results, since similar points inthe latent space do not contribute to a richer model. Thus a small threshold has
been employed which determines the minimum distance between 2 observation
points. Observations were then added incrementally to an initially empty vector
if the condition is fulfilled. Alternatively, one could also apply the IVM criteria,
described in [LSH02]. The number of training examples which we used finally was
N = 1120, comprising samples from 3 motions with 2 trials each.
For a further speedup the sparsification technique, as described in section 2.3.3
was employed. Altogether we were able to learn the model in 4 hours on a 2-Ghz
single CPU machine.
Since articulated body motion is non-linear we used LLE initialization for the
GPLVM as a baseline. Figure 4.15 shows the result of learning the latent spaceof the full model using a GPLVM initialized with LLE (a) and our method (b).
Note that our method, unlike the GPLVM initialized with LLE, leads to a smooth(i.e., consecutive frames in time are close in the latent space) result, and separates
well the different classes. As shown in the figure, smoothness implies a lower

57

CHAPTER 4. EXPERIMENTAL RESULTS

Figure 4.16: Learning 2 motions (blue = milling, green = rolling) using
our method. The upper-body skeleton (in red) depicts mean-poses
calculated at the corresponding points in the latent space. One cansee that ”grabbing the mill” (using the left hand, poses to the left) is
well separated from ”grabbing the rolling pin” (using the right hand,poses to the right).
relationship error. To quantify class separation, we also compute the Fisher score[UD07] defined asFscore = tr Sw−1 Sb
where Sw is the within class matrix
Sw =

C
X
Ni
i=1

N

"

N

i 

T
1 X
(i)
(i)
(i)
(i)
x −x
xj − xNi j=1 j#

and Sb is the between class matrix
Sb =

C
X
Ni
i=1

N


T
x(i) − x x(i) − x .

Here C denotes the number of classes (e.g. C = 3 in figure 4.15), Ni is the number(i)of points in class i and N is the total number of training points. Furthermore xj
stands for the j’th point in class i, x(i) is the mean of class i and x represents
the overall mean. Note that our method performs significantly better than the
GPLVM with LLE initialization in terms of both, the relationship error and thefisher score.
To gain more insight into the latent space, figure 4.16 shows the model after
optimizing training data from 2 types of motion. Here the mean pose is depict fordifferent points in the latent space.

58

CHAPTER 4. EXPERIMENTAL RESULTS

4.3.2 Feature extraction and image likelihood
Tracking was performed using silhouette features. In contrast to [UFF06], wherepoints in the video sequence (for example knees, hands, elbows, etc.) had to be
tracked by a separate 2D tracker, this appearance-based approach has the advantage of no need to being reinitialized manually (after occlusions) or relying on
a separate tracking algorithm. Disadvantages include that gradient-based refinements (for accurate tracking) are harder to incorporate and a more sophisticatedlikelihood had to be used than the one we used here. However, we found that with
our test data reliable tracking was possible with an easy-to-implement likelihood
described in this section.
To retrieve silhouettes from the input image, a small number of frames (we used
4 frames of the sequence) was labeled manually in each view (see figure 5.2).
Departing from the labeled frames, a background and a color representation has
been created. The background representation consists of storing all pixel values
belonging to the background for each labeled frame. The color representation
was created by storing a foreground and a background color histogram, where 16
bins per color (since each color channel provides 256 values, 16 values have beencondensed together) have been used. Thus, in total 16 × 16 × 16 = 4096 bins havebeen used.
The background segmentation was done by comparing all pixels of a 640 × 480
input image to all pixels (up to 4) in the background representation for the corresponding pixel coordinate. A pixel was then considered belonging to the foreground (or to the silhouette) if its Euclidian distance to all background pixels
at this location (in the RGB color space) was larger than 18 (which was foundemperically).
For the (skin) color segmentation we used the Bayesian classifier described in[PBC05] which employs a minimum cost decision rule: Classify a pixel with colorc as skin (or foreground silhouette), iff: P (skin|c) > P (skin|c). Since
P (skin|c) =P (c|skin)P (skin)
P (c)

and

P (skin|c) =P (c|skin)P (skin)P (c)we have
P (c|skin)
P (skin)>.
P (skin)P (c|skin)
Here P (c|skin) and P (c|skin) are given as counts by the previously collected color
histograms:
P (c|skin) =number of foreground pixels with same color
total number of foreground pixels

P (c|skin) =number of background pixels with same color
total number of background pixels

59

CHAPTER 4. EXPERIMENTAL RESULTS

Figure 4.17: Tracking setup. First, some (3-5) frames of the sequence werelabeled manually (a). Then a background and a color representation
has been generated automatically in order to segment the completesequence (b). The likelihood is based on a matching (green) between
the generative model and the segmentation in each frame (c). Finallya particle filter was used to track the sequence (d).
(skin)
was set to = 0.3, which was
For the experiments reported in this thesis PP (skin)found empirically.

We finally combined the background segmented image and the color segmented
image using a logical AND operator with respect to the foreground region. For
calculating the image likelihood p(z|y) (here z denotes an observation, which in our case - comes frome one or two color cameras) given a (30 dimensional)
data point y, the silhouette of the articulated body model (see figure 4.17) was
generated from y and rendered into the back buffer of an ATI graphics card
via OpenGL. By reading the graphics buffer back into memory, p(z|y) was finallyobtained as p(z|y) ∝ %, where % denotes the number of pixels (in the projected andthe segmented image) belonging to the same class (foreground or background).

4.3.3 Particle filter
Tracking itself (see figure 5.3) was performed using a simplified particle filter withsecond order Marcov dynamics in the latent space. We assumed that all observed
poses are captured by our pose model. The probabilistic framework of the tracker
is described next.
Let st be the state of our system at time t and let zt be the observation of
the cameras at time t. To shorten notation we abbreviate s1:t = {s1 , ..., st } andz1:t = {z1 , ..., zt }.
Following [IB98] we assume that the object dynamics forms a temporal Markov
chain, i.e.
p(st |s1:t−1 ) = p(st |st−1 ).

60

CHAPTER 4. EXPERIMENTAL RESULTS
Furthermore the observations zt are assumed to be conditional independent andindependent of the underlying dynamical processp(z1:t−1 , st |s1:t−1 ) = p(st |s1:t−1 )t−1Yp(zi |si )

i=1

which leads to the mutual conditional independence (by marginalizing st )p(z1:t |s1:t ) =tYp(zi |si ).

i=1

Using these assumptions, the rule for propagation of state density over time canbe formulated asp(st |z1:t ) = kt p(zt |st )p(st |z1:t−1 )whereZp(st |z1:t−1 ) =p(st |st−1 )p(st−1 |z1:t−1 )dst−1

and kt is a normalisation constant that does not depend on st . The rule canbe interpreted as the Bayes’ rule for inferring the posterior state density fromobservations for the time-varying case. The prove can be found in [IB98].
Since the observation density is non-Gaussian, the evolving state density p(st |z1:t )is also non-Gaussian. Applying a non-linear filter to evaluate the state density
cannot be done analytically. Thus we are using an approximation by evaluating
a finite set of particles.

Our state vector st consists of st = x(t) , x(t−1) , the estimation of the positionin the latent space at time t and time t − 1. Including the latent position attime t − 1 is needed for drawing samples from the second order Marcov dynamics.
Given the state vector st at time t we estimate the state vector st+1 at time t + 1
by creating P samples x̂i with i = {1, ..., P } in the latent space and calculating
the image likelihood of the mean projection p (zt+1 |µ(x̂i )) as described in section
4.3.2. Then st+1 is chosen as (t+1)  xargmaxx̂i [p (zt+1 |µ(x̂i ))]st+1 ==.
x(t−1)x(t)The samples x̂i were generated from the (predictive) Mixture-of-Gaussians distribution (see figure 4.18)





2
1
1
2
(t) 1
(t)
(t−1)
(t)
(t−1)
,
I + pKDE x(t)p(x̂i |x , x) = N x , I + N 2x − x
5
25
5
625
5
where the first term is a Gaussian distribution centered at the position of x attime step t, the second term represents a Gaussian distribution centered on the61CHAPTER 4. EXPERIMENTAL RESULTS

Figure 4.18: Predictive Distribution. This figure illustrates the dynamical
model. The predictive distribution is taken as a Mixture-of-Gaussians
with 3 components: (1) A random walk sampled directly around xt
to marginalize the danger of local minima, (2) a second order Marcov
model which extrapolates linearly from xt and xt−1 (Gaussian aroundblue point), and (3) Gaussians around training points, weighted bytheir corresponding kernel density estimate using a RBF kernel. Here2 trajectories with their training points are shown in the latent space.
Gaussians are depict by a circle. The line width of the circles indicate
their respective mixture weight.
second order Marcov estimate and the third term is a Kernel Density Estimate(KDE) of x at time t, defined as
pKDE (x(t) ) = PN

N
X

1

(t)
i=1 k(x , x̂i ) i=1

k(x(t) , x̂i )N (x̂i ,

1
I)
625

where N is the number of training points in the model. Here the inference kernelk(xi , xj ) is defined as a radial basis function:!
kxi − xj k22
k(xi , xj ) = θ1 exp −+ θ3 δij
2θ22
In our experiments we set the hyperparameters of the inference kernel to θ1 = 1,
θ2 = 0.2 and θ3 = 0.0001.
Since we are dealing with Gaussian Processes and large datasets, speed becomes
an issue. Calculating the mean pose µ(x̂i ) of a sample point x̂i using a fullGP involves a one time inversion of the covariance matrix K which scales asO(N 3 ). Furthermore, for each sample x̂i , a full KDE and N -dimensional matrixvector multiplication has to be performed which scales as O(N 2 ). To speed up
inference several approximation techniques could be used [SNS06, YDD05, UD08,
SW02, QCR05]. We chose one of the most intuitive ones, which makes use of thesparseness of K. As in [UD08] we infer the mean pose µ(x(i) ) using only a subset

62

CHAPTER 4. EXPERIMENTAL RESULTS(we chose the 50 nearest neighbors) of the training points, which in practice leadto a sufficiently good approximation.
Improvement in tracking accuracy could be gained by further optimizing the posterior p(s|z) as done in [UFF06], which allows for tracking poses beyond the ones
captured by the model. In this scenario, however, we found direct sampling fromthe model sufficient to obtain good tracking results (see figure 4.19).

4.3.4 Results
figure 4.19 depicts tracking and classification performance for milling and rolling
motions1 when using the models depicted in figure 4.15. We used the particle filter
described in section 4.3.3, that operates in the low dimensional latent space andmodels the dynamics with a second order Markov model. Our image likelihood
is based on low-level silhouette features (see section 4.3.2). The reason for using
a simple likelihood is that it was found sufficient for the experiments shown inthis thesis. This is mainly because most of the ambiguities are resolved by thestrong pose prior we use. For more complicated sequences however, designing amore sophisticated likelihood can improve results. This could be done by using
features based on optical flow, edges, appearance or incorporating an additional
discriminative model into the tracking process [UD08].
In practice we could achieve a sample evaluation rate of 100 frames per second on a
2 Ghz single CPU machine with hardware graphics acceleration. This means thatusing 25 particles and 2 cameras (which was found sufficient for reliable tracking
in our scenario, see figure 4.19) results in a tracking frame rate of 2 frames per
second, for example.
The tracking error was measured by averaging the positions of 6 joints (shoulder
left/right, elbow left/right, hand left/right). It is given as
6

1X
Tracking Error =kp̂ − pi k2
6 i=1 iwhere p̂i denotes the 3D position of the i’th joint as estimated by the particle
filter and pi stands for the ground truth.
For classification purposes, we labeled the dataset using the following 7 classes:• rest (both hands rest on the table next to the board)• grasp pin (reaching for the pin with the right hand)• rolling (perform a cyclic rolling motion with the rolling pin)
• grasp broom (grasp the broom and the dustpan with both hands)
1

No results are shown for brooming since no test data was available.

63

CHAPTER 4. EXPERIMENTAL RESULTS
Tracking milling sequences

Tracking rolling sequences

14

14
GPLVM+LLE, 1 camera
GPLVM+LLE, 2 cameras
our method, 1 camera
our method, 2 cameras

10
8
6
4
2
0
0

GPLVM+LLE, 1 camera
GPLVM+LLE, 2 cameras
our method, 1 camera
our method, 2 cameras

12
Tracking Error (cm)Tracking Error (cm)

12

10
8
6
4
2

20

40
60
80
Number of particles

0
0

100

20

(a)

Classifying rolling sequences
100

80
60
40
GPLVM+LLE, 1 camera
GPLVM+LLE, 2 cameras
our method, 1 camera
our method, 2 cameras

20

20

40
60
80
Number of particles

Classification Accuracy (%)

100
Classification Accuracy (%)

100

(b)

Classifying milling sequences

0
0

40
60
80
Number of particles

80
60
40

0
0

100

GPLVM+LLE, 1 camera
GPLVM+LLE, 2 cameras
our method, 1 camera
our method, 2 cameras

20

20

(c)

40
60
80
Number of particles

100

(d)

Figure 4.19: Tracking and nearest neighbor classification performance formilling and rolling motions using our method (red) and GPLVM initialized to LLE (blue) as a function of the number of particles usedin the particle filter.
• brooming (perform a brooming motion using the broom and the dustpan)• grasp mill (grasp the mill with the left hand)• milling (perform a cyclic milling motion)
Classification was done by computing the nearest neighbor of the current pose
estimate in the latent space, and assigning its class label to the current frame.
The classification accuracy is then given byClassification Accuracy =number of correctly labeled frames
× 100%.
total number of frames

Our method significantly outperforms the GPLVM with LLE initialization in bothtracking accuracy and NN classification performance (see figure 4.19) when using
both, one or two cameras. We observed that the milling motion is more ambiguousthan the rolling motion and thus can be tracked reliably only when using two
cameras, whereas one camera proved sufficient for the rolling motion.

64

CHAPTER 4. EXPERIMENTAL RESULTS
We furthermore observed that although the GPLVM model initialized to LLE
exhibits many jumps and discontinuities in the latent space, a particle filter canstill perform relatively well, depending on the neighborhood size which was used
for generating the Locally Linear Embedding. This is due to the fact that inorder to disambiguate motions with any of the models, the standard deviation
used for creating samples can not drop below a certain value. Thus jumps in thelatent space were found to be not as critical as initially assumed. In contrast, if
one wants to generate motion by sampling, discontinuities start to play a much
more important rule, leading to unrealistic motions when used for animation. We
also expect a tracking system which further optimizes the observation likelihood
together with the pose prior to behave significantly worse with a model like theone shown in figure 4.15a. Thus we think that a method like ours, which can
reduce dimensionality in a non-linear way without generating jumps in the latent
space can be worth further investigations.

65CHAPTER 4. EXPERIMENTAL RESULTS

66

5 Implementation
This project was implemented in C++ using Eclipse CDT and Ubuntu 8.04 - theHardy Heron. The following libraries were used:• FLTK (Fast Light Toolkit, cross-platform C++ GUI toolkit)• OpenGL (visualizing the latent space, rendering the generative model)• OpenCV (import, export and undistortion of images, 2D drawing)• SNOPT (constrained non-linear optimization for learning the model)• BLAS, LAPACK (efficient algorithms for linear algebra, e.g. SVD)
Figures 5.1, 5.2 and 5.3 show the GUI of the program created during this thesis.
First a script has been used to undistort images using the OpenCV function
cvUndistort2 with the intrinsic camera parameters found by the MATLAB camera
calibration toolbox 1 . Furthermore errors in the motion datasets were corrected bylooking for discontinuities along the latent trajectories. The errors occured during
registration of the motion capture dataset due to mislabeled markers. If a small
number of consecutive defective frames were detected, a linear interpolation wasused to interpolate in between. If a large number of consecutive frames was
defective the trial was rejected from being used for learning or tracking.
In the second step the motions were loaded into the program depict in figure 5.1
and a latent space was learned using the method described in section 3.3. This
program also allows for showing the quality of the model by placing a line into
the latent space and reconstructing the mean pose y from samples along this line.
For samples close to the training points those artificially generated motions shouldlook realistic, which indicates good interpolation and extrapolation abilities anda useful model for tracking or motion generation.
The third step consists of labeling 4 dissimilar frames of an image sequence which
is about to be tracked manually, as can be seen in 5.2. This program further
allows for creating the image statistics necessary for tracking and for checking theresult of a segmented test image.
In the final step, another program allows for tracking the scene using the particle
filter described in section 4.3.3 (see figure 5.3). The program parameters comprise the number of particles which are used for tracking as well as the number1

http://www.vision.caltech.edu/

67

CHAPTER 5. IMPLEMENTATION

Figure 5.1: Learning the latent space. This application allows for loading
a motion dataset and learning a latent space using different methods like GPLVM or GPLVM with Rank Prior. Initializations withgraph based techniques were obtained via the MATLAB toolkit from[vdMPvdH07]. The main window shows the current latent space, projected to the first three dimensions via PCA. It also allows for displaying the covariance plot and the particle sample distribution in
the latent space after a tracking experiment. The upper-right window
displays the evolution of the singular values over time during an optimization. The lower-right window depicts the current state of thearticulated human body model by applying the mean projection fromthe latent to the data space of a given test point.
of cameras. Since our test sequences behaved well, we found one or two cameras sufficient for this task. Tracking error and classification accuracy are stored
in text files. Furthermore images which show the tracking and classification result from different views are created. Videos of those results can be found at
http://www.rainsoft.de.
After tracking a sequence the particle distribution over time can be displayed in
the latent space by loading the samples and the model back again into the program
depict in figure 5.1.

68

CHAPTER 5. IMPLEMENTATION

Figure 5.2: Image likelihood and features. Prior to tracking, this application
asks the user to label the silhouette of the subject in 4 frames (upperleft). Using this information, a background- and a color-representationis created. The lower images show the result of applying those representations to a new frame (upper-right), resulting in a foregroundsegmentation (lower-left) and edges (lower-right). For tracking, theforeground representation was matched with a generative model ofthe articulated body.

69CHAPTER 5. IMPLEMENTATION

Figure 5.3: Tracking using a simple particle filter. After learning the modeland creating the feature representations this application performs the
tracking itself. A particle filter can be started, which is initialized bycomputing the training point with the highest likelihood. The tracking
step itself involves the generation of 5 to 100 particles (depending on
the experiment) using a second order Markov model. After computing
the mean prediction, the generative articulated body model silhuetteis read from the OpenGL back buffer and matched against the observation silhouette (for 1 or 2 cameras, depending on the experiment).
The sample with the highest likelihood is kept and the mean prediction
is displayed to the OpenGL front buffer.

70

6 Acknowledgments
I want to thank everybody who made my diploma thesis possible. First of all
many thanks to my supervisors1 who agreed on having me jointly working atUniversität Karlsruhe (TH) and Massachusetts Institute of Technology. They
helped me in every situation and answered all of my questions. Many thanks also
go to the DAAD (German Academic Exchange Service) who supported me with ascholarship for my 5-month stay in Boston, USA. Furthermore I want to thank myparents, friends and my girlfriend for their steady support and encouragement.

1Prof. Dr. Trevor Darrell studied at University of Pennsylvania and at MIT
where he received his Phd in 1996. After visiting Stanford University he firstbecame assistant professor (2000) and then associate professor (2003) at MIT
where he lead the Vision Interface Group. In August 2008 he moved to the
University of California in Berkeley, EECS & ICSI.
Dr. Raquel Urtasun studied telecommunication engineering at Universityof Navarra (UPNA) in Pamplona, Spain. She hold positions at UPNA, EURECOM and ENST (France) and received her Phd in 2006 from EPFL (Switzerland). After being postdoctoral associate with Trevor Darrell at CSAIL, MIT
she moved to the University of California in Berkeley, EECS & ICSI.
Prof. Dr. Alexander Waibel studied electrical engineering and computerscience at MIT and received his PhD in 1986 from CMU at Pittsburgh. He currently is professor at CMU and Universität Karlsruhe (TH) where he leads the”international center for Advanced Communication Technologies” (interACT).
His research interests include speech technologies and machine learning.
Dr. Rainer Stiefelhagen is assistant professor with Alexander Waibel atUniversität Karlsruhe (TH). He leads the visual perception group of interACT.
He studied at Universität Karlsruhe (TH) where he also received his Phd in2002 with a work about ”focus of attention”. His research is about visual
human perception, intelligent environments and machine learning.
Dipl.-Inform. Kai Nickel graduated from Universität Karlsruhe (TH) and
now works as a research assistant with Rainer Stiefelhagen. His research interests include visual person tracking and guesture recognition. he is involved in
the CHIL-project (computers in the human interaction loop) and in the SFB588(german collaborative research center 588) ”humanoid robots”.

71

CHAPTER 6. ACKNOWLEDGMENTS

72

List of Figures

2.1

2.2

2.3

2.4

Sampling from a Gaussian Process. (a) shows 3 functions
drawn at random from a GP prior. (b) depicts 3 functions drawn
at random from the posterior, i.e. the prior conditioned on thefive noise free observations indicated. The gray areas show theuncertainty for each point x respectively. This figure was takenfrom [RW06]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

Face plots of 2D Gaussians. Here face plots of 2 dimensional
Gaussians are shown for different parameters of the covariance matrix. While (a) is uncorrelated (Υ = 0), (b) is correlated with
Υ = 0.5. (c) exhibits almost perfect linear coherence (Υ = 0.99)
whereas (d) is negatively correlated (Υ = −0.8). . . . . . . . . . .

23

2D illustration of Gaussian Process regression. Here y1 isgiven and the correlation between y1 and y2 is assumed to be knownby means of a kernel function k(xi , xj ). Since y1 and y2 are correlated, the mean of the conditional distribution p(y2 |y1 ) is similar
to the mean y1 itself while the small variance indicates only little
uncertainty. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

5D illustration of Gaussian Process regression. Here 3 training points are given compared to only one in figure 2.3. Since wehave 2 test points, the regression results in a 5 dimensional Gaussian distribution. Note that uncertainty increases with the distanceof a test point from the training data. . . . . . . . . . . . . . . . .

24

73

List of Figures
2.5

3.1

3.2

Sparsification: Global and local GPs. This figure illustrates
the sparsification method used for motion datasets with a large
number of training points. The colored points are training points,
lying on a motion trajectory (black) in the data space. Instead of
inverting the full Gaussian Process covariance matrix (which takesO(N 3 ), with N the number of training points), we make use of priorknowledge about the nature of motion data, i.e. we assume smoothness over time (or frames). Thus only a set of smaller local GPs(which overlap with neighboring sets, shown in gray) is calculated,
reducing the computation time dramatically. To also keep track ofglobal relationships, each point in a local GP set is linked to one
point in each other local GP set via a global GP (which thus again
consists of a very limited number of training points, contributing
to a speed-up in learning the model). In this example 3 global GPs
are depicted (red, green and blue). . . . . . . . . . . . . . . . . .

30

Illustration using a spring-force system. From a simplified
viewpoint the method can be seen as a system of observation points(red) in the high-dimensional data space, connected by springs(green and blue) of strengths defined according to the distance ofthe two points they are connnecting, respectively. Thus the spring
strength can be defined by a radial basis function, for example.
Here, e.g., the blue spring is weaker than the two green springs since
the data points it connects are further apart from each other. Before
applying the algorithm, the springs are in their home position, imposing no force. The spring system stands for the GPLVM, which
preserves dissimilarities. On the other hand, the pressure plates
surrounding the system apply forces towards the covariance directions. They represent the Rank Prior which favors low-dimensional
spaces. Here the pressure along the smaller variance is strongest
since we want to get rid of the ”smaller” dimensions first. Thisis achieved by using a non-linear penalty function, as described in
section 3.3.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

Illustration of our Rank Prior with a GPLVM. The goal is torecover the 1D manifold embedded in 2D. The GPLVM gets stuck
in local minima very early (upper row) since PCA initialization does
not capture non-linear dependencies, whereas our method decreases
dimensionality gradually and recovers the correct manifold (lowerrow). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

74

List of Figures
3.3

4.1

4.2

4.3

4.4

4.5

4.6

4.7

4.8

Spectrum of a 30D motion database. (a),(b) and (c): Evolution of the first ten singular values as a function of the optimization iteration number for a linear, sigmoid and logarithmic sparsity
penalty functions. (d): Spectrums learned after convergence by different sparsity penalty functions compared to the observation space
spectrum (red). . . . . . . . . . . . . . . . . . . . . . . . . . . . .

39

Finding a 2D manifold in 3D space on a sparsely sampledswiss roll. Only a sparse, noisy subset (depicted in black) of thefull manifold is assumed to be known (a). (b) shows the initialization (with neighborhood size k=6), GPLVM result and 2D/3D
reconstruction of the full manifold (from top to bottom). . . . . .

42

Learning a sparsely sampled Swiss Roll. Here snapshots of
the learning process are shown at different iterations during optimization using our method. One can clearly see how the roll isunraveled correctly. The final manifold has 2 dimensions. . . . . .

44

Example: A sparsely sampled swiss roll. Reconstruction Error (a) and Relationship Error (b) for the experiment in figure 4.1
averaged over 20 random partitions of the data for the 8 best performing dimensionality reduction techniques. A more detailed result is given in figure 4.4. . . . . . . . . . . . . . . . . . . . . . . .

45

Example: A sparsely sampled swiss roll. Reconstruction
and Relationship Error for the experiment in figure 4.1 averaged
over 20 random partitions of the data. Here more detailed results
are shown, including PCA and graph-based methods with different
neighborhood sizes. The 8 best performing techniques are depicted
in 4.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

Dimensionality estimation. (Top) Five 2D manifolds embedded in 3D. (Bottom) Latent spaces and intrinsic dimensionalities Q
learned using our continuous dimensionality reduction method. . .

47

Learning a hat shaped manifold generated by sine functions. Snapshots of the learning process are shown at different
iterations during optimization using our method. The discovered
dimensionality Q is 2. . . . . . . . . . . . . . . . . . . . . . . . . .

48

Learning a hat shaped manifold with a hole. Snapshots of
the learning process are shown at different iterations during optimization using our method. The discovered dimensionality Q is
2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49

Learning a cut-off bowl manifold sampled along longitudal
lines. Snapshots of the learning process are shown at different
iterations during optimization using our method. The discovered
dimensionality Q is 2. . . . . . . . . . . . . . . . . . . . . . . . . .

50

75

List of Figures
4.9

Learning a 2D wave shaped manifold. Snapshots of the learning process are shown at different iterations during optimization
using our method. The discovered dimensionality Q is 2. . . . . .
4.10 Learning a spiral manifold embedded in 3D. Snapshots of
the learning process are shown at different iterations during optimization using our method. The discovered dimensionality Q is
1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.11 Learning a spiral manifold embedded in 3D. Snapshots of
the learning process are shown at different iterations during optimization using our method. The discovered dimensionality Q is
1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.12 Learning a narrow spiral manifold embedded in 3D. Snapshots of the learning process are shown at different iterations during
optimization using our method. The discovered dimensionality Q
is 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.13 Dataset used in the tracking experiment. Here an overview
of the dataset used in the experiment is given. First, the average
length per sequence is given in seconds. Next the total number ofutilizable sequences, the number of sequences used for learning themodel and the number of sequences used for tracking is shown. . .
4.14 Tracking sequences. For the tracking and classification experiments the database from [KPFW08] was used. It contains several
types of motion, captured with a VICON motion capturing systemand four Sony VWF 500 color cameras (resolution 640x480). In
this figure snapshots of 3 motions (rolling, brooming and milling)
are depicted from 2 views each. . . . . . . . . . . . . . . . . . . .
4.15 Learning different types of motion into one single 3D latent
space. Each type is visualized with a different color (red = rolling,green = brooming, blue = milling). (a) depicts a 3D-GPLVM initialized to LLE, discontinuities are emphasized in black. (b,c) show
the manifold learned by our method along with the variance. . . .
4.16 Learning 2 motions (blue = milling, green = rolling) using
our method. The upper-body skeleton (in red) depicts meanposes calculated at the corresponding points in the latent space.
One can see that ”grabbing the mill” (using the left hand, poses
to the left) is well separated from ”grabbing the rolling pin” (using
the right hand, poses to the right). . . . . . . . . . . . . . . . . .
4.17 Tracking setup. First, some (3-5) frames of the sequence were labeled manually (a). Then a background and a color representation
has been generated automatically in order to segment the completesequence (b). The likelihood is based on a matching (green) between the generative model and the segmentation in each frame(c). Finally a particle filter was used to track the sequence (d). . .

76

51

52

53

54

55

56

57

58

60

List of Figures
4.18 Predictive Distribution. This figure illustrates the dynamical model. The predictive distribution is taken as a Mixture-ofGaussians with 3 components: (1) A random walk sampled directly
around xt to marginalize the danger of local minima, (2) a second
order Marcov model which extrapolates linearly from xt and xt−1
(Gaussian around blue point), and (3) Gaussians around training
points, weighted by their corresponding kernel density estimate using a RBF kernel. Here 2 trajectories with their training points areshown in the latent space. Gaussians are depict by a circle. Theline width of the circles indicate their respective mixture weight.

62

4.19 Tracking and nearest neighbor classification performancefor milling and rolling motions using our method (red) and GPLVM
initialized to LLE (blue) as a function of the number of particles
used in the particle filter. . . . . . . . . . . . . . . . . . . . . . .

64

5.1

5.2

Learning the latent space. This application allows for loading
a motion dataset and learning a latent space using different methods like GPLVM or GPLVM with Rank Prior. Initializations withgraph based techniques were obtained via the MATLAB toolkit
from [vdMPvdH07]. The main window shows the current latent
space, projected to the first three dimensions via PCA. It also
allows for displaying the covariance plot and the particle sample
distribution in the latent space after a tracking experiment. The
upper-right window displays the evolution of the singular values
over time during an optimization. The lower-right window depicts
the current state of the articulated human body model by applying
the mean projection from the latent to the data space of a given
test point. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

68

Image likelihood and features. Prior to tracking, this application asks the user to label the silhouette of the subject in 4 frames(upper-left). Using this information, a background- and a colorrepresentation is created. The lower images show the result of applying those representations to a new frame (upper-right), resulting
in a foreground segmentation (lower-left) and edges (lower-right).
For tracking, the foreground representation was matched with agenerative model of the articulated body. . . . . . . . . . . . . . .

69

77

List of Figures
5.3

Tracking using a simple particle filter. After learning themodel and creating the feature representations this application performs the tracking itself. A particle filter can be started, which is
initialized by computing the training point with the highest likelihood. The tracking step itself involves the generation of 5 to
100 particles (depending on the experiment) using a second orderMarkov model. After computing the mean prediction, the generative articulated body model silhuette is read from the OpenGL
back buffer and matched against the observation silhouette (for 1
or 2 cameras, depending on the experiment). The sample with thehighest likelihood is kept and the mean prediction is displayed to
the OpenGL front buffer. . . . . . . . . . . . . . . . . . . . . . . .

78

70

Index
L1 -norm, 35

Derivative of the Rank Prior, 36
Dimensionality estimation, 47
Dimensionality Reduction, 13
Discontinuity, 64
Dynamical prior, 31
Dynamics, 57

Acknowledgments, 71
Algorithm (GPLVM + Rank Prior),
37
Ambiguity, 63
Approximated prior, 26

Eclipse CDT, 67
Eigendecomposition, 27
Eigenvalue, 28
Eigenvalue problem, 15, 27, 28
Eigenvalue spectrum, 28
Eigenvector, 28
Embedding, 16, 41
Energy of the spectrum, 36
Experimental results, 41

back-constrained GPLVM, 32
Background representation, 59
Basis funtion, 20
BLAS, 67
C++, 67
Camera calibration, 67
Chain rule for matrices, 36Class label, 63Classification, 63
Classification accuracy, 63Classifying human motion, 55
Color histogram, 59
Color representation, 59
Complexity, 26Computational complexity, 26, 27
Conditionally independence, 27, 60
Conjugate prior, 21
Conjugate property, 20
Consistency assumption, 19
Continuous dimensionality reduction,
33
Coppersmith-Winograd-algorithm, 26
Cost function, 33
Covariance function, 19
Covariance Matrix, 32
Cyclic topologies, 32

Feature extraction, 59
First order markov dynamics, 31
Fisher Score, 57
FLTK, 67
Function space view, 20
Gaussian distribution, 19
Gaussian Process, 19, 20
Gaussian Process Dynamical Model,
31
Gaussian Process Latent Variable Model,
24
Global neighborhood, 29
GP, 19GP prior, 20
GP regression, 22
GPDM, 31
GPLVM, 24, 25, 64
Graph based methods, 16

Darrell, Prof. Dr. Trevor, 71

79

IndexHyperparameter, 61
Hyperparameters, 21, 25

Log-likelihood, 25
Logarithmic sparsity penalty function,
35
Low dimensionality, 33
LTSA, 17, 41
LVM, 18

iid Gaussian noise, 24
Image likelihood, 59
Improper Prior, 33
Inducing inputs, 27
Inducing variables, 27
Inference kernel, 61
Infinite dimensional space, 21
Informative Vector Machine, 26
Initial guess, 24
Initialization, 37
Inner product space, 21
ISOMAP, 16, 25, 41
Isotropic Gaussian prior, 25
IVM, 26, 57

Manifold, 16, 41
Manual labeling, 67
Marginal distribution, 27
Marginal likelihood, 20, 26
Marginalization, 26
Marginalization property, 19, 27
Markov chain, 60
MATLAB, 67
Maximum Variance Unfolding, 18
Mean pose, 62
Mean projection, 61
Mislabeled marker, 67
Mixture-of-Gaussian, 61
MVU, 18

K nearest neighbors, 16
Karhunen-Loève transform, 14
KDE, 61, 62
Kernel Density Estimate, 61
Kernel function, 19, 25, 27
Kernel height, 21
Kernel trick, 21
Kernel width, 21
Kitchen, 55

Nearest neighbor classifier, 64
Nickel, Dipl.-Inform. Kai, 71
NLDR, 16
Non-degenerate kernel, 21
Non-linear Dimensionality Reduction,
16
Notation, 13
Nyström Method, 27, 28

LAPACK, 67
Laplace Beltrami operator, 16
Laplacian Eigenmaps, 16, 41
Latent space, 33
Latent Variable, 18
Latent Variable Model, 18
Lengthscale, 21
Linear + RBF kernel, 31
Linear Dimensionality Reduction, 14
Linear program, 18
Linear sparsity penalty function, 35
LLE, 17, 25, 32, 41, 57, 64
Local minimum, 24, 25, 27, 31
Local neighborhood, 16, 29
Local Tangent Space Alignment, 17
Locally Linear Embedding, 17, 32
Log posterior, 35

Observation density, 61
Observations, 60
OpenCV, 67
OpenGL, 67
Overfitting, 27
Overlapping neighborhood, 29
Particle distribution, 68
Particle filter, 60, 67
PCA, 14, 25, 37, 41
Posterior, 19, 25
PPCA, 18
Predictive distribution, 27, 61
Predictive equation, 22

80

IndexPrincipal Component Analysis, 14
Prior, 19, 25, 31
Prior over the dimensionality of thelatent space, 33
Probabilistic Principal Component Analysis, 18
Product of local and global GPs, 29
Radial Basis Function, 21
Rank Prior, 33, 57
RBF, 21
RBF + noise kernel, 21, 25
Reconstruction Error, 43
Regression, 22
Regularizer, 25
Relationship Error, 43
Rotation, 37

Swiss roll, 41
Tayler expansion, 17
Topological Prior, 32
Tracking, 67
Tracking accuracy, 63
Tracking error, 63
Tracking human motion, 55
Transition, 32
Ubuntu 8.04, 67
Unbiased covariance estimate, 33
Undistort, 67
Urtasun, Dr. Raquel, 71
VICON, 55
Video, 67
Waibel, Prof. Dr. Alexander, 71
Weight space view, 20
Wishart distribution, 25

Sample, 61
Sample evaluation rate, 63
Scaled Conjugate Gradients, 25
SCG, 25
Semidefinite program, 18
SIFT, 13
Sigmoid sparsity penalty function, 35
Silhouette features, 59
Singular value, 33
Singular Value Decomposition, 33
Singular vector, 33
Skin segmentation, 59
Smooth latent space, 31
SNOPT, 25, 37, 67
Sparse approximation, 26
Sparsification, 26, 57
Sparsity Penalty function, 33
Spectral methods, 16State density, 61
State vector, 61
Stiefelhagen, Dr. Rainer, 71
Stochastic Process, 19
Subsampling, 57
Subspace, 16Supervisor, 71
SVD, 3381

Index82

Bibliography[BN03]

Mikhail Belkin and Partha Niyogi. Laplacian eigenmaps for dimensionality reduction and data representation. Neural Computation,
15(6):1373–1396, June 2003.

[DHS00]R. O. Duda, P. E. Hart, and D. G. Stork. Pattern Classification,
2nd edition. Wiley-Interscience, 2000.

[GMHP04]

K. Grochow, S. Martin, A. Hertzmann, and Z. Popovic. Stylebased inverse kinematics. In Proceedings of SIGGRAPH 2004.
ACM Press / ACM SIGGRAPH, 2004.

[GMS02]

Philip E. Gill, Walter Murray, and Michael A. Saunders. Snopt:
An sqp algorithm for large-scale constrained optimization. SIAM
Journal on Optimization, 12(4):979–1006, 2002.

[Har07]

Stefan Harmeling. Exploring model selection techniques for nonlinear dimensionality reduction. Technical report, School of Informatics, Edinburgh University, Scotland, March 2007.

[IB98]

Michael Isard and Andrew Blake. Condensation – conditional density propagation for visual tracking. International Journal of Computer Vision, 29(1):5–28, 1998.

[KPFW08]

H. Koehler, M. Pruzinec, T. Feldmann, and A. Woerner. Automatic human model parametrization from 3d marker data formotion recognition. In Proceedings of International Conferences
in Central Europe on Computer Graphics, Visualization and Computer Vision, 2008.

[Law05]

Neil Lawrence. Probabilistic non-linear principal component analysis with gaussian process latent variable models. Journal of Machine Learning Research, 6:1783–1816, 2005.

[Law07]

N. D. Lawrence. Learning for larger datasets with the gaussian
process latent variable model. In M. Meila and X. Shen (eds)
Proceedings of the Eleventh International Workshop on Artificial
Intelligence and Statistics, 2007.

[LB05]

Elizaveta Levina and Peter J. Bickel. Maximum likelihood estimation of intrinsic dimension. In Lawrence K. Saul, Yair Weiss, andLéon Bottou, editors, Advances in Neural Information Processing
Systems 17, pages 777–784. MIT Press, Cambridge, MA, 2005.

83

Bibliography[Low99]David G. Lowe. Object recognition from local scale-invariant features. In Proc. of the International Conference on Computer Vision
ICCV, Corfu, pages 1150–1157, 1999.

[LQC06]

Neil D. Lawrence and Joaquin Quinonero-Candela. Local distancepreservation in the gp-lvm through back constraints. In ICML
’06: Proceedings of the 23rd international conference on Machinelearning, pages 513–520, New York, NY, USA, 2006. ACM.

[LSH02]

Neil D. Lawrence, Matthias Seeger, and Ralf Herbrich. Fast sparse
Gaussian process methods: The informative vector machine. In
Advances in Neural Information Processing Systems 15, pages 609–
616. MIT Press, 2002.

[mFSA07]

Amir massoud Farahmand, Csaba Szepesvári, and Jean-Yves Audibert. Manifold-adaptive dimension estimation. In ICML ’07:
Proceedings of the 24th international conference on Machine learning, pages 265–272, New York, NY, USA, 2007. ACM.

[Mol93]

Martin Fodslette Moller. A scaled conjugate gradient algorithm forfast supervised learning. Neural Networks, 6(4):525–533, 1993.

[PBC05]

S.L. Phung, Sr. Bouzerdoum, A., and Sr. Chai, D. Skin segmentation using color pixel classification: analysis and comparison. Pattern Analysis and Machine Intelligence, IEEE Transactions on,
27(1):148–154, Jan 2005.

[Pea01]

K. Pearson. On lines and planes of closest fit to systems of points
in space. Philosophical Magazine, 2(6):559–572, 1901.

[QCR05]

Joaquin Quinonero-Candela and Carl Edward Rasmussen. A unifying view of sparse approximate gaussian process regression. Journalof Machine Learning Research, 6:1939–1959, 2005.

[RS00]Sam T. Roweis and Lawrence K. Saul. Nonlinear Dimensionality
Reduction by Locally Linear Embedding. Science, 290(5500):2323–
2326, 2000.

[RW06]

Carl Edward Rasmussen and Christopher K. I. Williams. GaussianProcesses for Machine Learning. MIT Press, 2006.

[SNS06]

Yirong Shen, Andrew Ng, and Matthias Seeger. Fast gaussian
process regression using kd-trees. In Y. Weiss, B. Schölkopf, and
J. Platt, editors, Advances in Neural Information Processing Systems 18, pages 1225–1232. MIT Press, Cambridge, MA, 2006.

[SUF08]

M. Salzmann, R. Urtasun, and P. Fua. Local deformation models
for monocular 3d shape recovery. In Proceedings of the Conferencein Computer Vision and Pattern Recognition (CVPR) 2008. IEEE
Computer Society, 2008.

84

Bibliography[SW02]

Matthias Seeger and Christopher K. I. Williams. Fast forward
selection to speed up sparse gaussian process regression. In Advances in Neural Information Processing Systems 14. Cambridge,
MA, 2002.

[TB99]

Michael E. Tipping and Christopher M. Bishop. Probabilistic principal component analysis. Journal Of The Royal Statistical Society
Series B, 61(3):611–622, 1999.

[TSL00]

Joshua B. Tenenbaum, Vin de Silva, and John C. Langford. A
Global Geometric Framework for Nonlinear Dimensionality Reduction. Science, 290(5500):2319–2323, 2000.

[UD07]

Raquel Urtasun and T. Darrell. Discriminative gaussian processlatent variable models for classification. In ICML ’07: Proceedingsof the 24th international conference on Machine learning, June
2007.

[UD08]

Raquel Urtasun and Trevor Darrell. Local probabilistic regression for activity-independent human pose inference. In Proceedingsof the Conference in Computer Vision and Pattern Recognition(CVPR) 2008, June 2008.

[UFF06]

Raquel Urtasun, David J. Fleet, and Pascal Fua. 3d people tracking with gaussian process dynamical models. In Proceedings of theConference in Computer Vision and Pattern Recognition (CVPR)
2006, pages 238–245, Washington, DC, USA, 2006. IEEE Computer Society.

[UFG+ 08]

R. Urtasun, D. J. Fleet, A. Geiger, J. Popovic, T. Darrell, andLawrence N.D. Topologically-constrained latent variable models.
In ICML ’08: Proceedings of the 25th international conference on
Machine learning, 2008.

[vdMPvdH07] L. J. P. van der Maaten, E. O. Postma, and H. J. van den Herik. Dimensionality reduction: A comparative review, unpublished. Matlab implementation available for all methods., 2007.
[WFH08]

Jack M. Wang, David J. Fleet, and Aaron Hertzmann. Gaussian process dynamical models for human motion. IEEE Transactions on Pattern Analysis and Machine Intelligence, 30(2):283–298,
2008.

[WS01]

Christopher K. I. Williams and Matthias Seeger. Using the
Nyström method to speed up kernel machines. In T. Leen, T. Dietterich, and V. Tresp, editors, Neural Information Processing Systems 13, pages 682–688. MIT Press, 2001.

[WS06]

Kilian Q. Weinberger and Lawrence K. Saul. An introduction to
nonlinear dimensionality reduction by maximum variance unfold-85

Bibliography
ing. In Proceedings of the Twenty First National Conference on
Artificial Intelligence (AAAI), 2006.
[YDD05]

Changjiang Yang, Ramani Duraiswami, and Larry Davis. Efficient kernel machines using the improved fast gauss transform. In
Lawrence K. Saul, Yair Weiss, and Léon Bottou, editors, Advancesin Neural Information Processing Systems 17, pages 1561–1568.
MIT Press, Cambridge, MA, 2005.

[ZZ03]

Zhenyue Zhang and Hongyuan Zha. Nonlinear dimension reduction
via local tangent space alignment. In IDEAL 2003: International
conference on intelligent data engineering and automated learning,
pages 477–481, 2003.

86