Topological Constraints and Kernel-Based Density
Estimation
Carl Henrik EkHedvig KjellströmDanica Kragic ∗School of Computer Science and Communication
KTH Royal Institute of Technology, Stockholm, Sweden{fpokorny, chek, hedvig, danik}@csc.kth.se

Florian T. Pokorny

Abstract
This extended abstract1 explores the question of how to estimate a probability distribution from a finite number of samples when information about the topologyof the support region of an underlying density is known. This workshop contribution is a continuation of our recent work [1] combining persistent homology
and kernel-based density estimation for the first time and in which we explored an
approach capable of incorporating topological constraints in bandwidth selection.
We report on some recent experiments with high-dimensional motion capture data
which show that our method is applicable even in high dimensions and develop
our ideas for potential future applications of this framework.

1

Introduction

In recent years, novel topological techniques such as persistent homology have been developed to extract homologicalinformation from point-cloud data. This opens up a whole
new range of possibilities for exploiting topological methods
in a data-driven paradigm. Exciting theoretical results, providing guarantees on the reconstruction of homological inforFigure 1: A topologically interesting mation from point-cloud data for large sample sizes have been
point-cloud S ⊂ R2 .
obtained [2] and a growing community of researchers has recently found interesting applications of the persistence algorithm, such as the analysis of computer vision and medical data in [3]. However, much of themotivation for this research has been coming from a theoretical point of view or has been focussed
only on trying to extract topological information from point-cloud data without attempting to build
a probabilistic model for the data. In the work presented here, we ask the question of how else persistent homology might be used as a tool for machine learning and suggest motion capture data as apotential domain of application. We present some preliminary experiments showing that our topological density estimation framework [1], which combines classical kernel-based density estimation[4] with persistent homology, can be applied even to high-dimensional motion capture data. Insteadof trying to recover topological information from a finite number of samples, we are interested in thecase where we are given prior topological information about the support region of some otherwise
unknown probability distribution and where we would then like to utilize this information together
with the persistence algorithm to select an optimal bandwidth parameter.
∗This work was supported by the EU projects FLEXBOT (FP7-ERC-279933) and TOMSY (IST-FP7270436) and the Swedish Foundation for Strategic Research
1
Appearing in: Workshop on Algebraic Topology and Machine Learning, Advances in Neural Information
Processing Systems (NIPS) 25, 201212Background

Many probabilistic methods are based on Gaussian densities. Examples include Gaussian Mixture
Models (GMM) [5], Gaussian Processes (GP) [6] and Gaussian Mixture Regression (GMR). Whilethese methods are highly popular in application fields such as robotics, speech recognition and computer vision [5, 7, 8], a crucial and often overlooked property of such models is that they all generatedensities on Rd for which supp f = Rd ; i.e. such models will assign a non-zero probability to every subset of Rd with non-zero volume. This can cause problems in several real-world applications.
Consider for example a probabilistic model of human pose positions in joint space. The two human hands, for example, are often in close contact during manipulation and grasping, yet we should
assign zero probability to configurations corresponding to collisions with the environment or to selfintersections of the various body parts. Such constraints are impossible to enforce directly using
Gaussian densities. A better approach is to base such a model on densities with bounded support
which can be achieved using kernel-based density estimation [9, 10, 11]. There, one tries to reconstruct a probability density f : Rd → R of which only an i.i.d. sample
 S = {X1 , · · · , Xn } is known
Pnx−Xi1ˆˆusing an estimator fε,n defined by fε,n (x) = nεd i=1 K, and where the kernel function
εK : Rd → R is a suitably chosen probability density which is concentrated near the origin. We willin particular focus on spherical kernels that are symmetric functions of the norm kxk of their input
variable x ∈ Rd and which satisfy supp K = B1 (0), where Bε (p) = {x ∈ Rd : kx − pk 6 ε}.
SnThen supp fˆε,n = Yε (S), where Yε (S) = i=1 Bε (Xi ). By choosing ε small enough, we can hence‘design’ a support region that is as concentrated as we want around the data points and which does
not suffer from the Gaussians’ trivial support region which we alluded to earlier. In kernel-based
density estimation ε is called the bandwidth [4]. Theorem 2.1 provides reassurance that, for large
sample sizes, fˆε,n → f pointwise if the bandwidth is chosen correctly:
Theorem 2.1 (Theorem 3.1 and 3.2 [11]). Supposethat K is a Borel
measurable funcR
Rtion on Rd such that supx∈Rd |K(x)| < ∞, Rd |K(x)| dx < ∞, Rd K(x)dx = 1 andlimkxk→∞ kxkd K(x) = 0. Suppose that {εn }∞n=1 is a sequence of positive numbers such thatlimn→∞ εn = 0. Then limn→∞ E[fˆεn ,n (x)] = f (x) at every point of continuity of f . If furthermore limn→∞ nεdn = ∞, then limn→∞ E[(fˆεn ,n (x) − f (x))2 ] = 0 at every point of continuity off.
0.2
0.1

0.20.1102030102030Figure 2: Graph of a density function f (shaded, grey), a selection of 50 samples (displayed below the x-axis)
and a traditional Gaussian kernel-based estimate (σ 2 = 14 , black line) with a bandwidth that is too large (ε = 5,on the left) and too small (ε = 0.4, on the right) respectively.

The question of optimal bandwidth selection is an ongoing area of research [12] and optimality
is traditionally studied in an asymptotic L1 or L2 -error context [4, 13]. Given knowledge of theunderlying density f and a spherically symmetric kernel K, the asymptotic mean integrated squarederror is given by:
ZZ
1
ε4n22AM ISE(εn ) = dK(x) dx + µ2 (K){tr(Hess f (x))}2 dx,
nεn4R 2where µ2 (K) =xi K(x)dx is independent of j ∈ {1, . . . , d} by the spherical symmetry. AM ISE(εhn ) provides an asymptotici estimate for the natural mean integrated squared error
R2M ISE(ε) = E (fˆε,n (x) − f (x)) dx if limn→∞ εn = 0 and limn→∞ nεdn = ∞ [4]. A popularbranch of bandwidth selection methods attempts to determine bandwidths which minimize AMISE[12] instead of working with the more complicated MISE directly.
Considering the point cloud S displayed in Figure 1, a human observer – and especially anyone
familiar with persistent homology – will immediately notice that S intuitively has three connectedcomponents, one of which has a ‘hole’. We can reformulate the existence of the ‘hole’ in Figure
1 in a mathematically precise way using the machinery of persistent homology by inspecting the2corresponding barcode. We will assume here that the specialized audience of this workshop isfamiliar with the basic concepts from persistent homology ([3] provides an excellent introduction
to the field). Observe that, in persistent homology, we attempt to reconstruct information about theBetti numbers of an underlying space X from the samples S by studying the space Yε (S) as ε > 0varies – just like we vary ε in the case of bandwidth selection with spherical kernels. In practice,
the homology groups of Yε (S) are then approximated by computing the homology groups of theVietoris-Rips complex V2ε (S).
We consider three kernels with supp K = B1 (0) ⊂ Rd [1]. All these kernels are zero outsided(d+1)Γ( d)2
B1 (0) and defined by Ku = Vol(B1 (0))−1 (uniform), Kc (x) =(1 − kxk) (conic) anddKt (x) = (2πσ 2 )−d22π 2 −11kxk2,
Γ d−21 − Γ2 d2σe 2σ2 (truncated Gaussian) respectively for x ∈ B1 (0) (see(2)Figure 3 a-c).

(a) Ku(b) Kc(d) fˆ, Ku(c) Kt(e) fˆ, Kc(f) fˆ, Ku

Figure 3: (a-c): ε12 K( xε ), for ε = 4 and K = Ku , Kc , Kt respectively. (d-e): Density approximation for

K = Ku , Kc , Kt in dimension 2 and with ε = 4 and for a few data points in the plane. For Kt , σ 2 = 14 .

In [1], we introduced the idea of topological constraints to bandwidth selection. i.e., given a density
f such that some constraints on the Betti numbers of supp f and an i.i.d. sample S are known, we
compute the barcode for S and search for barcode intervals such that these constraints are satisfied.
We then chose the first such interval I = [εmin , εmax ] (which also coincided with the largest such
1min (n) − 4+dand were able
ninterval in our application) and defined εtop (n) = εmin (n) + εmax (n)−ε2
to show in examples that the empirical MISE for this bandwidth choice converges quickly to zero
as n → ∞ for Ku , Kt , Kc . We believe that this novel type of topological bandwidth selection, i.e.
choosing a bandwidth εn such that topological constraints are satisfied and such that the expectedsquared error is minimized for a given sample size n is a very interesting topic, which we believe
should be explored further both from an experimental and a theoretical perspective.

3Experiments

To gain an intuition for our approach, let us consider first a new 2D synthetic example of our topological density estimation framework. We consider the density f displayed Figure 4 a). Observe thatΩ = supp f has three connected components (b0 (Ω, Z2 )) = 3) and that two of these components
have a 2-dimensional hole each (b1 (Ω, Z2 ) = 2). Taking this observation as our prior topologicalknowledge about f during reconstruction, we are now able to compute a density estimator as displayed in 4 b) from only 300 samples and using εtop . We can also clearly see in that figure and in 4c) that the reconstructed density has the correct topological features. In this example, εmin u 1.36,εmax u 2.23 and εtop u 1.53. The reconstruction L1 and L2 errors for εtop were approximately
0.55 and 0.04 respectively. A very interesting topic that remains to be explored more fully andwhich we look forward to discussing during this workshop is the question of how else one could
choose εtop ∈ [εmin , εmax ] in order to minimize L1 or L2 errors. In [1], we determined an εtop
which ensures point-wise asymptotic convergence under some assumptions on εmin , εmax . In many
application domains, one might however prefer to choose εtop as small as possible in order not to’generalize’ the data too much. As outlined in the introduction, the probabilistic modelling of motion
is an interesting area of application for our density estimation approach which we would like to begin
to explore in this work. For this purpose, we generated a point-cloud S = {X1 , . . . , X2276 } ⊂ R60 ,where each vector Xi is a concatenation of 20 vectors in R3 . Each of these 20 vectors describes
the (x, y, z) positions of one of twenty key/joint positions of a person at a fixed point in time (seethe stick-person in Figure 5 for a visualization of one such Xi ). We start with a short 35 frame
motion-capture sequence of a person walking along a straight line taken from [14]. We then repeat
these frames and translate and rotate each frame to obtain semi-synthetic data of a person walking
3

(a)

(b)

(c)

Figure 4: a) a 2D density. b) a uniform kernel estimate based on 300 samples and with εtop = 1.53. c) thesupport region of the kernel estimate in b) with the correct number of ‘holes’ and components.

along the patterns in the plane displayed in Figure 5. Three of these paths are periodic closed-curve
walking patterns, two are straight line walking segments and the last is a single isolated static posture. As given topological data, we assume the knowledge that b0 (Ω, Z2 ) = 6 (i.e. there are 6
underlying motion patterns) and that b1 (Ω, Z2 ) = 3 (i.e. there are 3 ‘cyclic’ motions). This example
is typical in that we do not have knowledge of any underlying probability density and can hence notevaluate the error incurred by our estimator. We do know however that a probabilistic model should
not ’mix up’ the various motion patterns by creating a density with too few connected components
in its support region. We use JavaPlex [15] to determine the smallest bandwidth εmin such thatb0 (V2εmin (S), Z2 ) = 6 and b1 (V2εmin (S), Z2 ) = 3. In our experiment εmin u 12. In this experiment, we choose εtop = εmin as a conservative estimate for the bandwidth. This yields a densityfˆ with support equal to Y12 (S) from which we can generate samples by selecting a data point withuniform probability and then sampling from the rescaled kernel K at that data point. Here, we chose
K = Kc . Note that we can now continuously move from one data point to the next inside Y12 (S) ineach of the 6 connected components (e.g. to potentially obtain an interpolation of the whole walking
sequence from the samples). We observe also that, unlike the projection of the paths in Figure 5, thepaths in the full 60 dimensional space do not intersect.
Figure 5 shows a data point from S in 5 a) and 3 random samples from a conical kernel scaled by thebandwidth εtop = 12 and centred at that data point 5 b-d). We note that, even though our dataset isvery sparse in 60 dimensions, our framework enables us to generalize from the data while satisfying
the given topological constraints.
300
60

60

60

60

40

40

40

40

20

20

20

20

200
100
0
−100
0
−2060
80

00
−20

60
80

00
−2060

0
−2080

060
80

020 100

20 10020 10020 100

(a)

(b)

(c)

(d)

−200
−200

−150

−100

−50

050

100

150

200

250300(e)

Figure 5: Stick-person from the data (a) and three samples from a conical kernel with bandwidth 12
centred at that data-point (b-d). (e) displays the movement patterns in the xy-plane along which wegenerated our data points.

4Conclusion and Future WorkIn this work, we have begun to explore applications of our novel topological density estimation
framework [1] to a new synthetic data set in 2D and to motion capture. From our evaluation of
the estimation errors in [1] and the encouraging results with real motion data, we believe that thesenovel uses of persistent homology in kernel-based density estimation can be a valuable tool in realworld applications. We are intending to investigate various aspects of our approach in the comingmonths. Incorporating a time-dependence in the motion-capture data is an obvious starting point forexample.
4References[1] F. T. Pokorny, C. H. Ek, H. Kjellström, and D. Kragic, “Persistent homology for learning
densities with bounded support,” in Advances in Neural Information Processing Systems 25,pp. 1826–1834, 2012.
[2] 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.
[3] G. Carlsson, “Topology and data,” Bull. Amer. Math. Soc. (N.S.), vol. 46, no. 2, pp. 255–308,2009.
[4] M. P. Wand and M. C. Jones, Kernel Smoothing. Chapman and Hall/CRC, 1 ed., Dec. 1994.
[5] 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.
[6] C. E. Rasmussen and C. Williams, Gaussian Processes for Machine Learning. MIT Press,
2006.
[7] S. Calinon and A. Billard, “Incremental learning of gestures by imitation in a humanoid robot,”in ACM/IEEE International Conference on Human-Robot Interaction, 2007.
[8] D.-S. Lee, “Effective Gaussian mixture learning for video background subtraction,” PAMI,
vol. 27, no. 5, pp. 827–832, 2005.
[9] M. Rosenblatt, “Remarks on some nonparametric estimates of a density function,” The Annalsof Mathematical Statistics, vol. 27, no. 3, pp. 832–837, 1956.
[10] E. Parzen, “On estimation of a probability density function and mode,” Ann. Math. Statist.,
vol. 33, pp. 1065–1076, 1962.
[11] T. Cacoullos, “Estimation of a multivariate density,” Annals of the Institute of Statistical Mathematics, vol. 18, pp. 179–189, 1966.
[12] B. A. Turlach, “Bandwidth Selection in Kernel Density Estimation: A Review,” in CORE andInstitut de Statistique, pp. 23–493, 1993.
[13] L. Devroye, Nonparametric Density Estimation: The L1 View. Wiley, 1985.
[14] A. Agarwal and B. Triggs, “3d human pose from silhouettes by relevance vector regression,”
Computer Vision and Pattern Recognition, IEEE Computer Society Conference on, vol. 2,pp. 882–888, 2004.
[15] A. Tausz, M. Vejdemo-Johansson, and H. Adams, “JavaPlex: A software package for computing persistent topological invariants.” Software, 2011.

5