Persistent Homology for Learning Densities with
Bounded Support

Florian T. Pokorny
Carl Henrik EkHedvig Kjellström
Danica Kragic ∗Computer Vision and Active Perception Lab, Centre for Autonomous Systems
School of Computer Science and Communication
KTH Royal Institute of Technology, Stockholm, Sweden{fpokorny, chek, hedvig, danik}@csc.kth.se

Abstract
We present a novel method for learning densities with bounded support whichenables us to incorporate ‘hard’ topological constraints. In particular, we show
how emerging techniques from computational algebraic topology and the notionof persistent homology can be combined with kernel-based methods from machinelearning for the purpose of density estimation. The proposed formalism facilitateslearning of models with bounded support in a principled way, and – by incorporating persistent homology techniques in our approach – we are able to encode
algebraic-topological constraints which are not addressed in current state of theart probabilistic models. We study the behaviour of our method on two synthetic
examples for various sample sizes and exemplify the benefits of the proposed approach on a real-world dataset by learning a motion model for a race car. We show
how to learn a model which respects the underlying topological structure of the
racetrack, constraining the trajectories of the car.

1

Introduction

Probabilistic methods based on Gaussian densities have celebrated successes throughout machinelearning. They are the crucial ingredient in Gaussian mixture models (GMM) [1], Gaussian processes [2] and Gaussian mixture regression (GMR) [3] which have found applications in fields suchas robotics, speech recognition and computer vision [1, 4, 5] to name just a few. While Gaussian
distributions are convenient to work with for several theoretical and practical reasons (the central
limit theorem, easy computation of means and marginals, etc.) they do fall into the class of densities
on Rd for which supp f = Rd ; i.e. they assign a non-zero probability to every subset with non-zerovolume in Rd . This property of Gaussians can be problematic if an application dictates that certain
subsets of space should constitute a ‘forbidden’ region having zero probability mass. A simple example would be a probabilistic model of admissible positions of a robot in an indoor environment,
where one wants to assign zero – rather than just ‘low’ – probability to positions corresponding tocollisions with the environment. Encoding such constraints using e.g. a Gaussian mixture model is
not natural since it assigns potentially low, but non-zero probability mass to every portion of space.
In contrast to the above Gaussian models, we consider non-parametric density estimators based on
spherical kernels with bounded support. As we shall explain, this enables us to study topologicalproperties of the support region Ωε for such estimators. Kernel-based density estimators are wellestablished in the statistical literature [6] with the basic idea being that one should put a rescaled
version of a given model density over each observed data-point to obtain an estimate for the probability density from which the data was sampled. The choice of rescaling – or ‘bandwidth’ – ε
has been studied with respect to the standard L1 and L2 error and is still an active area of research[7]. We focus particularly on spherical truncated Gaussian kernels here which have been some∗This work was supported by the EU projects FLEXBOT (FP7-ERC-279933) and TOMSY (IST-FP7270436) and the Swedish Foundation for Strategic Research1what overlooked as a tool for probabilistic modelling. An important aspect of these kernels is that
their associated conditional and marginal distributions can be computed analytically, enabling us toefficiently work with them in the context of probabilistic inference.
A different interpretation of a density with support in an ε-ball can be given using the notion ofbounded noise. There, one assumes that observations are distorted by noise following a density withbounded support (instead of e.g. Gaussian noise). Bounded noise models are used in the signal
processing community for robust filtering and estimation [8, 9], but to our knowledge, we are thefirst to combine densities with bounded support and topology to model the underlying structureof data. Thinking of a set of observations S = {X1 , ..., XnS} ⊂ Rn as ‘fuzzy up to noise in an
ε-ball’ naturally leads one to consider the space Ωε (S) = i Bε (Xi ) of balls of size ε aroundthe data points. Persistent homology is a novel tool for studying topological properties of spaces
such as Ωε (S) which has emerged from the field of computational algebraic topology in recent
years [10, 11]. Using persistent homology, it becomes possible to study clustering, periodicity andmore generally the existence of ‘holes’ of various dimensions in Ωε (S) for ε lying in an interval.
Starting from the basic observation that one can construct a kernel-based density estimator fˆε whose
region of support is exactly Ωε (S), this paper investigates the interplay between the topologicalinformation contained in Ωε (S) and a corresponding density estimate. Specifically, we make thefollowing contributions:• Given prior topological information about supp f = Ω, we define a topologically admissible bandwidth interval [εmin , εmax ] and propose and evaluate a topological bandwidth
selector εtop ∈ [εmin , εmax ].
• Given no prior topological information, we explain how persistent homology can be of use
to determine a topologically admissible bandwidth interval.
• We describe how additional constraints defining a forbidden subset F ⊂ Rn of theparameter-space can be incorporated into our topological bandwidth estimation framework.
• We provide quantitative results on synthetic data in 1D and 2D evaluating the expected L2 errors for density estimators with topologically chosen bandwidth values ε ∈{εmin , εmid , εmax , εtop }. We carry out this evaluation for various spherical kernels andcompare our results to an asymptotically optimal bandwidth choice.
• We use our method in a learning by demonstration [12] context and compare our resultswith a current state of the art Gaussian mixture regression method.

2

Background2.1

Kernel-based density estimation

Let S = {X1 , ..., Xn } ⊂ Rd be an i.i.d. sample arising from a probability density f : Rd → R.
Kernel-based density estimation [13, 14,P
15] is an approachfor reconstructing f from the sample
iby means of an estimator fˆε,n (x) = nε1d ni=1 K x−X,
where
the kernel function K : Rd → R
εis a suitably chosen probability density. In this context, ε > 0 is called the bandwidth. If one isonly interested in an estimator that minimizes the expected L2 norm of fˆε,n − f , the choice of ε iscrucial, while the particular choice of kernel K is generally less important [7, 6]. Let {εn }∞n=1 bea sequence of positive bandwidth values depending on the sample size n. It follows from classicalresults [14, 15] that for any sufficiently well-behaved density K, limn→∞ E[(fˆεn ,n (x) − f (x))2 ] = 0
provided that limn→∞ εn = 0 and limn→∞ nεdn = ∞. Despite this encouraging result, the questionof determining the best bandwidth for a given sample is an ongoing research topic and the interested
reader is referred to the review [7] for an in-depth discussion.h One branch of methodsi [6] tries toR2ˆminimize the Mean Integrated Squared Error, M ISE(εn ) = E (fεn ,n (x) − f (x)) dx .
An asymptotic analysis reveals that, under mild conditions on K and f [6], M ISE(εn ) can be approximated asymptotically by AM ISE(εn ) as n → ∞ if limn→∞ εn = 0 and limn→∞ nεdn = ∞.
Here, AMISE denotes the Asymptotic Mean Integrated Squared Error. If we consider only spherical kernels that are symmetric functions of the norm kxk of their input variable x, an asymptotic
analysis [6] shows that, in dimension d,
Z
Z
1
ε4n
2
2
2
AM ISE(εn ) = d
K(x) dx + µ2 (K){tr(Hess f (x))} dx,
nεn
4
2

R
where µ2 (K) = x2j K(x)dx is independent of the choice of j ∈ {1, . . . , d} by the sphericalsymmetry and tr(Hess f (x)) denotes the trace of the Hessian of f at x. Due to the availability of a
relatively simple explicit formula for AMISE, a large class of bandwidth selection methods attempt
to estimate and minimize AMISE instead of working with MISE directly. One finds that AMISE is
minimized for
1
! 4+d
R
d K(x)2 dx
1εamise (n) =.
R
n µ2 (K)2 {tr(Hess f (x))}2 dx
Since f is assumed unknown in real world examples, so called plug-in methods can be used toapproximate εamise [7]. In this paper, we will work with two synthetic examples of densities forwhich we can compute εamise numerically in order to benchmark our topological bandwidth selection procedure. For our experiments, we choose three spherical kernels K : Rd → R that are
defined to be zero outside the unit ball B1 (0) and are defined by Ku = Vol(B1 (0))−1 (uniform),
Kc (x) =d(d+1)Γ( d
)
2
d

2 −d
2

(1 − kxk) (conic) and Kt (x) = (2πσ )

2π 2




 −1
kxk2
Γ d
, 1
−
2 2σ 21− Γ de 2σ2(2)(truncatedGaussian) respectively for kxk 6 1. These kernels can be defined in any dimension d > 0 and arespherical, i.e. they are functions of the radial distance to the origin only which enables us to efficiently evaluate them and to sample from the corresponding estimator fˆε,n even when the dimensiond

kxk2

d is very large. We will denote the standard spherical Gaussian by Ke (x) = (2πσ 2 )− 2 e− 2σ2 .

(a) Ku(c) Kt , σ 2 = 14(b) KcFigure 1: 412 K( x4 ) for the indicated kernels and a corresponding estimator fˆ4,3 for three sample points.
2.2

Persistent homology

Consider the point cloud S shown in Figure 2(a). For a human observer, it is noticeable that S looks‘circular’. One can reformulate the existence of the ‘hole’ in Figure 2(a) in a mathematically precise
way using persistent homology [16] which has recently gained increasing traction as a tool for the
analysis of structure in point-cloud data [10].

(a) Ω0(b) Ω0.25(c) Ω0.5(d) b0(e) b1Figure 2: Noisy data concentrated around a circle (a) and corresponding barcodes in dimension zero (d) andone (e). In (b) and (c), we display Ωε for ε = 0.25, 0.5 respectively together with the corresponding VietorisRips complex V2ε which we use for approximating the topology of Ωε . While the vertical axis in the ithbarcode has no special meaning, the horizontal axis displays the ε parameter of V2ε . At any fixed ε value,
the number of bars lying above and containing ε is equal to the ith Betti number of V2ε . The shaded region
highlights the ε-interval for which V2ε has one connected component (i.e. b0 (V2ε ) = 1) in (d) and for which asingle ‘circle’ (i.e. b1 (V2ε ) = 1) is detected in (e).

In the approach of [10], one starts with a subset Ω ⊂ Rd and assumes that there exists some probability density f on Rd that is concentrated near Ω. Given an i.i.d. sample S = {X1 , · · · , Xn } fromthe corresponding probability distribution, one of the aims of persistent homology in this setting is
to recover some of the topological structure of Ω – the homology groups Hi (Ω, Z2 ), for i = 1, . . . , d– from the sample S. Each Hi (Ω, Z2 ) is a vector space over Z2 and its dimension bi (Ω) is called3the ith Betti number. One of the properties of homology is that homology groups are invariant undera large class of deformations (i.e. homotopies) of the underlying topological space. A popular example of such a deformation is to consider a teacup that is continuously deformed into a doughnut.
One can think of b0 (Ω) as measuring the number of connected components while, roughly, bi (Ω),
for i > 0 describes the number i-dimensional holes of Ω. A closed curve in Rd that does not selfintersect can for example be classified by b0 = 1 (it has one connected component) and b1 = 1(it is topologically a circle). The reader is encouraged to consult [17] for a rigorous introduction tohomotopies and related concepts.
Sn
Given a discrete sample S and a distance parameter ε > 0, consider the set Ωε (S) = i=1 Bε (Xi ),for ε ∈ [0, ∞), where Bε (p) = {x ∈ Rd : kx − pk 6 ε}. In Figure 2(b) and 2(c) this set is displayed
for increasing ε values. Ωε (S) is a topological space and, in the case where Ω is a smooth compact
submanifold in Rd and f is in a very restrictive class of densities with support in a small tubularneighbourhood around Ω, [18, 11] have proven results showing that Ωε (S) is homotopy equivalentto Ω with high probability for certain large sample sizes. The key insight of persistent homology isthat we should study not just the homology of Ωε (S) for a fixed value of ε but for all ε ∈ [0, ∞)simultaneously. The idea is then to study how the homology groups Hi (Ωε (S), Z2 ) change with εand one records the changes in Betti number using a barcode [10] (see e.g. figure 2(d) and 2(e)).
Computing the barcode corresponding to Hi (Ωε (S), Z2 ) directly (via the Čech complex given by
our covering of balls Bε (X1 ), . . . , Bε (Xn ) [10]) is computationally very expensive and one hencecomputes the barcode corresponding to the homology groups of the Vietoris-Rips complex V2ε (S).
This complex is an abstract complex with vertices given by the elements of S and where we insert
a k-simplex for every set of k + 1 distinct elements of S such that any two are within distance lessthan 2ε of each other (see [10]). The homology groups of V2ε (S) are not necessarily isomorphic to
the homology groups of Ωε (S), but can serve as an approximation due to the interleaving propertyof the Vietoris-Rips and Čech complex, see e.g. Prop 2.6 [10]. For the computation of barcodes, we
use the javaPlex software [19]. The computed ith barcode then records the birth and death times of
topological features of V2ε in dimension i as we increase ε from zero to some maximal value M ,where M is called the maximal filtration value.

3

Our framework

Given a dataset S = {X1 , . . . , Xn } ⊂ Rd , sampled in an i.i.d. fashion from an underlying probability distribution with density f : Rd → R with bounded support Ω, we propose to recover fusing a kernel density estimator fˆε,n in a way that respects the algebraic topology of Ω. For this, weconsider only fˆε,n based on kernels K with supp K = B1 (0), and in particular, we experiment with
Kt , Ku and Kc . For such kernels, supp fˆε,n = Ωε (S) = ∪ni=1 Bε (Xi ) whose topological features
we can approximate by computing the barcodes for V2ε .
If no prior information on the topological features of Ω is given, we can then inspect these barcodesand search for large intervals in which the Betti numbers do not change. This approach is used in[10], who demonstrated that topological features of data can be discovered in this way. Alternatively,
one might be given prior information on the Betti numbers (e.g. using knowledge of periodicity,
number of clusters, inequalities involving Betti numbers) that one can incorporate by searching for εintervals on which such constraints are satisfied. Geometric constraints on the data can additionallybe incorporated by restricting to allowable ε-intervals to values for which Ωε (S) does not contain‘forbidden regions’. In the robotics setting, frequently encountered examples for such forbidden
regions are singular points in the joint space of a robot, or positions in space corresponding tocollisions with the environment.
Let us now assume that we are given constraints on some of the Betti numbers of Ω. For a given
sample S, we then compute the barcodes for V2ε in each dimension i ∈ {1, . . . , d} up to a large
maximal value M using javaPlex [19] and determine the set A of admissible ε values. If A is empty,
we consider the topological reconstruction to have failed. This will happen, for example, if ourassumptions about the data are incorrect, or if we do not have enough samples to reconstruct Ω.
If A is non-empty, we attempt to determine a finite union of disjoint intervals on which the Betti
numbers constraints are satisfied. Since, in our experiments, the interval I = [εmin (n), εmax (n)](determined up to some fixed precision) with smallest possible εmin (n) among those coincidedwith the largest such interval in most cases (indicating stable topological features), we decided to4

investigate this I ⊂ A for further analysis. For ε ∈ [εmin (n), εmax (n)], the resulting density fˆε,nthen has a support region Ωε (S) with the correct Betti numbers – as approximated by V2ε . We notethe following elementary observation:
Lemma 3.1. Let d ∈ N and εmin (n), εmax (n) ∈ R for all n ∈ N. Suppose that limn→∞ εmin (n) = 0
and that there exists a, b ∈ R such that 0 < a < εmax (n) < b and 0 6 εmin (n) < εmax (n) for all1
min (n) − 4+dn ∈ N. Then εtop (n) = εmin (n) + εmax (n)−εnsatisfies i) εtop (1) = εmid (1) and εtop (n) ∈2min (n)ii) limn→∞ εtop (n) = 0[εmin (n), εmid (n)] for all n ∈ N, where we define εmid (n) = εmax (n)+ε2dand iii) limn→∞ nεtop (n) = ∞.

It is our intuition that, for a large class of constraints on the Betti numbers and for tame densities
f : Rd → R (such as densities concentrated on a neighbourhood of a compact submanifold of Rd[11]), εmin (n) and εmax (n) exist for all large enough sample sizes n with high probability andthat the conditions of Lemma 3.1 are satisfied. In that case, Lemma 3.1 provides a motivation forchoosing {εtop (n)}∞n=1 as a topological bandwidth selector since – while it is difficult to analyse
εmin (n) asymptotically – at least the second summand of εtop (n) has the same asymptotics inn as the optimal AMISE solution. Furthermore, this choice of bandwidth then corresponds to asupport region Ωεtop (n) (S) with the correct Betti numbers (as approximated by the Vietoris-Rips
complex) since εtop (n) ∈ [εmin (n), εmax (n)]. Finally, ii) and iii) then imply that, point-wise,
limn→∞ E[(fˆεtop (n),n (x) − f (x))2 ] = 0 due to the results of [14, 15].
We note here that many different methods for choosing ε(n) ∈ [εmin (n), εmax (n)] can be considered. If the topologically admissible interval [εmin (n), εmax (n)] is for example determined by
the constraint of having three connected components of supp f as in 3(a), εmax (n) will increase
if we shift the connected components of supp f further apart. εtop (n) hence also increases andmight not yield good L2 error results for small sample sizes anymore. In that case, an estimator ε̂top (n) ∈ [εmin (n), εmax (n)] closer to εmin (n) might be a better choice. To give an initial
overview, we hence also display results for εmin (n), εmid (n), εmax (n) in our experiments. Notehowever also that the L2 error might not be the right quality measure for applications where thetopological features of supp f are most important – we illustrate an example of this situation in our
racetrack data experiment. We will show that – in the absence of further problem-specific knowledge – εtop (n) does yields a good bandwidth estimate with respect to the L2 error in our examples.

4

Experiments

Results in 1D We consider the probability density f : R → R displayed in grey in each of thegraphs in Figure 3. To benchmark the performance of our topological bandwidth estimators, we thencompute the AMISE-optimal bandwidth parameter εamise numerically from the analytic formula forf and for Kt , Ku , Kc and Ke . Here, we include the Gaussian kernel Ke for comparison purposes
only.
0.2

0.2

0.2

0.1

0.1

0.1

0

0

0

5

10

15

20

25

(a) fεtop ,10 using Kt .

30

0

5

10

15

20

25

(b) fεamise ,10 using Ke

30

0

0

5

10

15

20

25

30

(c) fεtop ,2500 using Kt

Figure 3: Density f (grey) and reconstructions (black) for the indicated sample size, bandwidth and kernel.
In order to topologically reconstruct f , we then assume only the knowledge of some points sampledfrom f and that b0 (supp f ) = 3 and no further information about f , i.e. we assume to know asample and that the support region of f has three components. We then find εtop (n) by computing
a topologically admissible interval [εmin (n), εmax (n)] from the barcode corresponding to the given
sample. To evaluate the quality of bandwidth parameters chosen inside [εmin (n), εmax (n)], we thensample at various sampling sizes and compute the mean L2 errors for the resulting density estimatorfε,n for ε = εtop , εmin , εmax and εmid = 21 (εmax + εmin ) for each of the spherical kernels that
we have described and compare our results to εamise . We set σ 2 = 14 for Ke and Kt . The results,
summarized in Figure 4, show that εtop performs at a level comparable to εamise in our experiments.
Note here that εamise can only be computed if the true density f is known, while, for εtop , we only
5

3.5

Εtop
Εmin

3.0

0.30

0.25

Εmid
2.5

Εmax

0.30

top
amise
min
mid
max

top
amise
min
mid
max

0.25

0.30

top
amise
min
mid
max

0.25

0.30

top
amise
min
mid
max

0.25

0.20

0.20

0.20

0.20

0.15

0.15

0.15

0.15

0.10

0.10

0.10

0.10

0.05

0.05

0.05

0.05

Εamise, Kt

2.0Εamise, KeΕamise, Ku

1.5Εamise, Kc
1.0
0.5
0.0

1

10

20

30

(a) bandwidth values

0.00

1

10

20

(b) Kt , σ

2

0.00

30

= 14

1

10

(c) Ke , σ

20

2

30

= 14

0.00

1

10

20

30

0.00

(d) Ku

1

10

20

30

(e) KcFigure 4: We generate samples from our 1D density using rejection sampling and consider sample sizes n from10 to 100 in increments of 10 (small scale) and from 250 to 5000 in increments of 250 (larger scale), resulting
in 30 increasing sample sizes n1 , . . . , n30 . In order to obtain stable results, we perform the sampling for eachsampling size 1000 times (small scale), 100 times (for 250, 500, 750, 1000) and 10 times (for n > 1000)
respectively. We then compute the corresponding kernel density estimators fˆε,n and the mean L2 norm off − fˆεn ,n . Figures (b)-(e) display these mean L2 errors (vertical axis) for the indicated kernel function andbandwidth selectors. Figure (a) displays the bandwidth values (vertical axis) for the given bandwidth selectors.
In all the above plots, a horizontal coordinate of i ∈ {1, . . . , 30} corresponds to a sample size of ni .

10

0

-10
-10

(a) density f010(b) 100 samples andΩεtop in grey.

0(c) fˆεtop ,100 using just
100 samples as in5(b)3(d) barcode
for b003(e) barcode
for b1Figure 5: 2D density, samples with inferred support region Ωεtop , topological reconstruction (using Kt ,σ 2 = 41 ) and barcodes with [εmin , εmax ] highlighted.

required the information that b0 (supp f ) = 3. In our experiments (sample sizes n > 10), we wereable to determine a valid interval [εmin (n), εmax (n)] in all cases and did not encounter a case where
the topological reconstruction was impossible.
Results in 2D Here, we consider the density f displayed in Figure 5(a). We chose this example to be representative for problems also arising in robotics, where the localization of a robot canbe modelled as depending on a probability prior which encodes space occupied by objects by zeroprobability. In such scenarios, we might be able to obtain topological information about the unobstructed space X, such as knowing the number of components or holes in X. Such information
could be particularly valuable in the case of deformable obstacles since their homology stays invariant under continuous deformations by homotopies. We set up the current experiment in a fashion
similar to our 1D experiments, i.e. we iterate sampling from the given density for various sample
sizes and compute the resulting mean L2 errors to evaluate our results. As we can see from Figure
6, our results indicate that bandwidths ε ∈ [εmin , εmax ] yield errors comparable with the AMISE
optimal bandwidth choice. While εtop does not perform as well as in the previous experiment, we
can observe that the corresponding L2 errors nonetheless follow a decreasing trend. Note also thatboth in 1D and 2D, εtop also yields good L2 error results for the standard spherical Gaussian kernel here. In applications such as probabilistic motion planning, the inferred structure of supp f ishowever of importance as well (e.g. since path-connectedness of supp f is important), making abounded support kernel a preferable choice (see also our racetrack example).
6

2.0

0.08

Εtop

Εmid

1.5

0.06

Εmax

0.08

top
amise
min
mid
max

Εmin

0.08

top
amise
min
mid
max

0.06

top
amise
min
mid
max

0.06

0.08

top
amise
min
mid
max

0.06

Εamise, KtΕamise, Ke

1.0

Εamise, Ku

0.04

0.04

0.04

0.04

0.02

0.02

0.02

0.02

Εamise, Kc
0.5

0.0

100

500

1000

0.00

1500

(a) bandwidth values

100

500

1000

(b) Kt , σ

2

0.00

1500

= 14

100

500

1000

(c) Ke , σ

2

0.00

1500

= 14

100

500

1000

(d) Ku

1500

0.00

100

500

1000

1500

(e) KcFigure 6: We generate samples from our 2D density using rejection sampling and consider sample sizes from100 to 1500 in increments of 100. We perform sampling 10 times for each sample size and compute thecorresponding kernel-based density estimator fˆε,n and the mean L2 norm of f − fˆεn ,n . Figures (b)-(e) display
these mean L2 errors (vertical axis) for the indicated sample size (horizontal axis) and kernel function. Figure(a) displays the indicated bandwidth values (vertical axis) and sample size (horizontal axis).
200

200

......
.....................
....
......
..
...
.............
..
......
.........
..
...
..
....
........
.........
.....
........
....
.......................
..........
...
............
...
.......
....
.......
........
.........
....
..............
.......
... . .
...................... ... . ...........
....
..............
......
................
.................................... .. . . . .................
.
.
.
.
.
.
.
.
.
...........
. . .............................
...........
......................
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
..........
....... .
.............. ...............
.
.
.
.
.
.
.
....
.
.
.
.
.
...
.
.
.
.
.
.
.
.
..
.
.
.........
..........
.....
.......
.....
......
.....
....
.................
.....
....
..
...
......
....
....
.......................
.....
.....
....... ...............
.........
... .
.........
.. .......
.....
.......
...........
....
.........
.......
........
.............
....... ..........
.....
.
.... ......
.
.
.
.
.
.........
..
.. . .....
....
.....
..........
....... ........
...........
........
........
........
......... ..............
......... ..................................
.. ..
.
......
.
.
.......
......... ................. ..............
....
........ ......... ...............
.
.
.
.
...
. .... . ..
......
...... ....... .................
.....
......
....
.....
............ ............
.............
.......
............
......
...
..... .......
......
......
...........
......
.....
... . .............
.... ......
..........
...........
.......
.........
... .....
........
..............
......
......... ............
.........
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.......... ...... ....
................................
. ...
.............................
. .....
.........
..................... ...
.... . ....
...........
.........
....
..............
.....
.....
.....
............ ...................
.......................
....

150

100

50

0

0

50

100

150

150

100

50

200

(a) Position component of ourracetrack data

0

0

50

100

150

200

(b) Projection of inferred support region, generated vector
field and sample trajectories(c) Inferred vector field, position
likelihood and sample trajectories using GMR.

Figure 7: Figure (a) shows the positions of a race car driving 10 laps around a racetrack. In (b), the resultsof our proposed method are displayed while Figure (c) shows the standard GMR approach. We exploit the
topological information that a racetrack should be connected and ‘circular’ when learning the density. As canbe seen, our model correctly infers the region of support as the track (grey). Using GMR, on the other hand, anon-zero probability is assigned to each location. We observe that the most probable regions are also lying over
the track (black being more likely). However, when sampling new trajectories using the learned density, we cansee that, whereas the trajectories using our method are confined to the track, the GMM results in undesirable
trajectories.

Application to regression We now consider how our framework can be applied to learn complex
dynamics given a topological constraint. We consider GPS/timestamp data from 10 laps of a race car
driving around a racetrack which was provided to us by [20]. For this dataset (see Figure 7(a)), weare given no information on what the boundaries of the racetrack are. One state of the art approach tomodelling data like this is to employ a learning by demonstration [12] technique which is prominent
especially in the context of robotics, where one attempts to learn motion patterns by observing a
few demonstrations. There, one uses data points S = {(Pk , Vk ) ∈ R2n , k = 1 . . . n}, where Pk
describes the position and Vk ∈ Rn the associated velocity at the given position. In order to modelthe dynamics, one can then employ a Gaussian mixture model [12] in R2n to learn a probability
density fˆ for the dataset (usually using the EM-algorithm). To every position x ∈ Rn , one canthen associate the velocity vector given by E(V |P = x) with respect to the learned density fˆ –this uses the idea of Gaussian mixture regression (GMR). The resulting vector field can then be
numerically integrated to yield new trajectories. Since E(V |P = x) for a Gaussian mixture model
can be computed easily, this method can be applied even in high-dimensional spaces. While it canbe considered as a strength of the GMR approach that it is able to infer – from just a few examples –7

a vector field that is non-zero on a dense subset of Rn , this can also be problematic since geometricand topological constraints are not naturally part of this approach and we cannot easily encode thefact that the vector-field should be non-zero only on the racetrack.
From our GPS/timestamp data, we now compute velocity vectors for each data-point and embedthe data in the manner just described in R4 . We then experimented with the software [21] to model
our racetrack data with a mixture of a varying number of Gaussians. While the model brakes down
completely for a low number of Gaussians, some interesting behaviour can be observed in the case ofa mixture model with 50 Gaussians displayed in Figure 7(c). We display the resulting velocity vector
field together with several newly synthesized trajectories. We observe both an undesired periodic
trajectory as well as a trajectory that almost completely traverses the racetrack before converging
towards an attractor. The likelihood of a given position is additionally displayed in 7(c) with black
being the most likely. While the most likely positions do occur over the racetrack, the mixture
model does not provide a natural way of determining where the boundaries of the track should lie.
The topmost trajectory in 7(c), for example, starts at a highly unlikely position.
Let us now consider how we can apply the density estimation techniques
we have described in this paper in this case. Given that we know thatthe racetrack is a closed curve, we assume that the data should be modelled by a probability density f : R4 → R whose support region Ω hasa single component (b0 (Ω) = 1) and Ω should topologically be a circle(b1 (Ω) = 1). In order for the velocities of differing laps around the track
not to lie too far apart , and so that the topology of the racetrack dominatesin R4 , we rescale the velocity components of our data to lie inside the in0
2.5
5 0
2.5
5
terval [−0.6, 0.6]. Figure 8 displays the barcode for our data. Using our(a) b0(b) b1procedure, we compute that [εmin , εmax ] u [3.25, 3.97] is the bandwidthinterval for which the topological constraints that we just defined are satFigure 8:
Barcodes isfied. Using the kernel K with σ 2 = 1 and the corresponding densityt4
in dimension zero (a)4ˆestimatorf,weobtainΩ
εtop
εtop ⊂ R with the correct topological propand one (b) and shaded[εmin , εmax ] interval for erties. Figure 7(b) displays the projection of Ωεtop onto R2 . As a next
our racetrack.
step, we suggest to follow the idea of the GMR approach to compute theposterior expectation E(V |P = x), but this time for our density fˆεtop . It
follows from the definition of our kernel-based estimator
that, for x such that (x, y) ∈ Ωεtop for


Pn

n

some y ∈ R , we have E(V |P = x) =i=1 Yi

Pn

R

i=1

x−X

i ,z dz
Kt ε
top


x−Xi
Kt ε
,z dz

R

. While we were not able to find a ref-top

erence for the use or computation of these marginals for spherical truncated Gaussians, a reasonably
simple calculation shows that these can in fact be computed analytically in arbitrary dimension:
Lemma 4.1. Consider d, k ∈ N, d > k and x ∈ Rk . Let Kt : Rd → R denote the spherical truncatedGaussian with parameter σ 2 > 0. ThenZ
Kt (x, y)dy =Rd−k2

P ( d−k, 1−kxk) − kxk22
1
2
2σ 2
e 2σ
d
1
2
k/2
(2πσ )P ( 2 , 2σ2 )denotes the normalized Gamma P function.
for kxk 6 1 and 0 otherwise. Here, P (a, b) = 1 − Γ(a,b)Γ(a)

For every point in the projection of Ωεtop onto the position coordinates, we can hence compute a
velocity E(V |P = x) to generate new motion trajectories. For points outside the support region, we
postulate zero velocity. Figure 7(c) displays the resulting vector-field and a few sample trajectories.
As we can see, these follow the trajectory of the data points in Figure 7(a) very well. At the sametime, the displayed support region looks like a sensible choice for the position of the racetrack.

5

Conclusion

In this paper, we have presented a novel method for learning density models with bounded support. The proposed topological bandwidth selection approach allows to incorporate topologicalconstraints within a probabilistic modelling framework by combining algebraic-topological information obtained in terms of persistent homology with tools from kernel-based density estimation.
We have provided a first thorough evaluation of the L2 errors for synthetic data and have exemplified
the practical use of our approach through application in a learning by demonstration scenario.
8

References[1] D. A. Reynolds, T. F. Quatieri, and R. B. Dunn, “Speaker verification using adapted Gaussian
mixture models,” Digital Signal Processing, vol. 10, no. 1–3, pp. 19–41, 2000.
[2] C. E. Rasmussen and C. Williams, Gaussian Processes for Machine Learning. MIT Press,
2006.
[3] D. A. Cohn, Z. Ghahramani, and M. I. Jordan, “Active learning with statistical models,” Journal of Artificial Intelligence Research, no. 4, pp. 129–145, 1996.
[4] S. Calinon and A. Billard, “Incremental learning of gestures by imitation in a humanoid robot,”in ACM/IEEE International Conference on Human-Robot Interaction, pp. 255–262, 2007.
[5] D.-S. Lee, “Effective Gaussian mixture learning for video background subtraction,” PAMI,
vol. 27, no. 5, pp. 827–832, 2005.
[6] M. P. Wand and M. C. Jones, Kernel Smoothing, vol. 60 of Monographs on Statistics andApplied Probability. Chapman and Hall/CRC, 1995.
[7] B. A. Turlach, “Bandwidth selection in kernel density estimation: A review,” in CORE andInstitut de Statistique, pp. 23–493, 1993.
[8] L. El Ghaoui and G. Calafiore, “Robust filtering for discrete-time systems with bounded noiseand parametric uncertainty,” IEEE Transactions on Automatic Control, vol. 46, no. 7, pp. 1084–
1089, 2001.
[9] Y. C. Eldar, A. Ben-Tal, and A. Nemirovski, “Linear minimax regret estimation of deterministic parameters with bounded data uncertainties,” IEEE Transactions on Signal Processing,vol. 52, no. 8, pp. 2177–2188, 2008.
[10] G. Carlsson, “Topology and data,” Bull. Amer. Math. Soc. (N.S.), vol. 46, no. 2, pp. 255–308,
2009.
[11] P. Niyogi, S. Smale, and S. Weinberger, “A topological view of unsupervised learning fromnoisy data,” SIAM Journal of Computing, vol. 40, no. 3, pp. 646–663, 2011.
[12] S. M. Khansari-Zadeh and A. Billard, “Learning stable non-linear dynamical systems with
Gaussian mixture models,” IEEE Transaction on Robotics, vol. 27, no. 5, pp. 943–957, 2011.
[13] M. Rosenblatt, “Remarks on some nonparametric estimates of a density function,” The Annalsof Mathematical Statistics, vol. 27, no. 3, pp. 832–837, 1956.
[14] E. Parzen, “On estimation of a probability density function and mode,” Annals of Mathematical
Statistics, vol. 33, pp. 1065–1076, 1962.
[15] T. Cacoullos, “Estimation of a multivariate density,” Annals of the Institute of Statistical Mathematics, vol. 18, pp. 179–189, 1966.
[16] H. Edelsbrunner, D. Letscher, and A. Zomorodian, “Topological persistence and simplification,” Discrete Comput. Geom., vol. 28, no. 4, pp. 511–533, 2002.
[17] A. Hatcher, Algebraic Topology. Cambridge University Press, 2002.
[18] P. Niyogi, S. Smale, and S. Weinberger, “Finding the homology of submanifolds with highconfidence from random samples,” Discrete Comput. Geom., vol. 39, no. 1-3, pp. 419–441,
2008.
[19] A. Tausz, M. Vejdemo-Johansson, and H. Adams, “JavaPlex: A software package for computing persistent topological invariants.” Software, 2011.
[20] KTH Racing, Formula Student Team, KTH Royal Institute of Technology, Stockholm, Sweden.
[21] A. Billard, “GMM/GMR 2.0.” Software.

9