Topological Trajectory Clustering with Relative Persistent Homology
Florian T. Pokorny, Ken Goldberg and Danica Kragic

Abstract— Cloud Robotics techniques based on Learningfrom Demonstrations suggest promising alternatives to manualprogramming of robots and autonomous vehicles. One challenge
is that demonstrated trajectories may vary dramatically: it canbe very difficult, if not impossible, for a system to learn control
policies unless the trajectories are clustered into meaningfulconsistent subsets. Metric clustering methods, based on a
distance measure, require quadratic time to compute a pairwise
distance matrix and do not naturally distinguish topologicallydistinct trajectories. This paper presents an algorithm for
topological clustering based on relative persistent homology,
which, for a fixed underlying simplicial representation anddiscretization of trajectories, requires only linear time in thenumber of trajectories. The algorithm incorporates global
constraints formalized in terms of the topology of sublevel orsuperlevel sets of a function and can be extended to incorporate
probabilistic motion models. In experiments with real automobile and ship GPS trajectories as well as pedestrian trajectories
extracted from video, the algorithm clusters trajectories intomeaningful consistent subsets and, as we show in an experimentwith ship trajectories, results in a faster and more efficient
clustering than a metric clustering by Fréchet distance.

I. I NTRODUCTION
Very large datasets of vehicle, robot, and human trajectories from GPS and other sensors are increasingly available
in the Cloud [26], [9]. Such Big Data has great potential for
Cloud Robotics [19] techniques such as automated Learningfrom Demonstrations to provide an alternative to manualprogramming of robots and autonomous vehicles. This paper
presents a new algorithm using results from topology to efficiently and effectively partition trajectories into meaningfulclusters.
Some of the key difficulties in applying machine learning
methods to large databases of trajectories arise from thefact that trajectories are naturally of varying length andhence not easily representable in a vector space of fixed
dimension – which is the natural domain for popular methods
such as Support Vector Machines. Furthermore, probabilistic
approaches to motion analysis such as Gaussian Processes[20] are smooth in nature, often making the formulation of
discrete constraints, such as whether a trajectory passes an
obstacle to the right or to the left, difficult to encode. However, for an autonomous system to learn motion primitives
and control policies, consistent subsets of training trajectory data that do not intermingle global behavior patterns
Florian T. Pokorny and Ken Goldberg are with the Departmentof Computer Science and Electrical Engineering, University of California, Berkeley. Ken Goldberg is also with the Department of Industrial Engineering and Operations Research, University of California, Berkeley. Danica Kragic is with CAS/CVAP, KTH Royal Instituteof Technology ftpokorny@berkeley.edu, dani@kth.se,

goldberg@berkeley.edu

Fig. 1. Distinguishing topologically distinct trajectories can facilitate learning of control policies for automated driving, flying, and robot motion. As an
example, this figure illustrates 177 automobile trajectories based on recorded
GPS coordinates at a highway intersection (data from OpenStreetMap [26]).
One can observe many distinct trajectories travelled between the lower left
shaded entry zone and upper right shaded exit zone. The presented algorithm
for topological trajectory clustering with relative persistent homology efficiently identified 24 distinct trajectory classes (colored accordingly) using
the complex DCr (X, Y ) at r = 0.000109. Eight of these trajectory classes
are shown in the lower figures.

are beneficial. The topological clustering method proposed
here (see Fig. 1) holds promise in this respect because
we are able to respect global topological properties of the
trajectory data and environment. Following the discovery of
persistent homology [12], the incorporation of topological
techniques with Machine Learning is now beginning toreceive increased attention, for example at recent workshops
at the flagship machine learning conferences ICML and NIPS[17], [25]. While trajectories can be clustered geometricallyusing several distance measures such as the Hausdorff orFréchet metric, or using techniques such as string kernels [6],[31] and Dynamic Time Warping, the complexity of typicaldistance based approaches scales quadratically in the size ofthe trajectory dataset, making these methods challenging toapply with very large databases of trajectories. Furthermore,distance based clustering is rather sensitive to the chosen distance measure, with each such measure having benefits anddrawbacks. The Hausdorff distance, for example, is highly
sensitive to point-wise differences between trajectories, while
L2 based averaging approaches are on the other hand rather
insensitive to point-wise outliers, leading to ongoing research
into trajectory distance measures [33]. In this work, we
propose a topological rather than purely geometric approach
to clustering trajectories into consistent subsets with potential future applications to anomaly detection and policy
learning for robotics and autonomous driving. We show inexperiments with GPS trajectories that the proposed method
can yield fundamentally different clusters compared to aclustering by Fréchet distance and show how our approach
can be combined with probabilistic reasoning by considering
sublevel sets of a probability density function for a pedestrian
motion dataset.
II. BACKGROUND AND R ELATED W ORK
The processing, classification and clustering of trajectorydata is a sizable research area. Some of the current methods
are reviewed in [33]. In the robotics domain, trajectories,
recorded for example as sequences of joint-angles of arobotic arm, play an important role in the learning fromdemonstration framework [1], [5], [27]. There, trajectory
sequences are recorded during a demonstration phase, wherea human instructor ‘teaches’ a robot how to perform acertain task. The robot then uses the trajectory data to model
motion primitives that can be adapted based on environment
conditions. In [14], the automated extraction of clusters of
trajectories in order to obtain vocabularies of motion wasstudied using a system comprising filtering, segmentation andclustering using K-means. Jenkins et al. [18], focused on theextraction of behaviors for humanoid motion in particular.
Knepper et al. [21] studied classes of path segments toenable a robot to reason about motion alternatives and toplan trajectories efficiently. Trajectories have been clusteredusing local sub-trajectories in the work of Lee [22] and also
Buchin [7], who focused on extracting commuting patterns
using the Fréchet distance.
Generally, previous clustering methods are either geometric in nature such as [22] or are based on probabilistic models, as in [15], who applied a mixture model for trajectory
clustering and [24], who applied a hidden Markov model todetect activities from trajectories. Our work is distinct fromthese approaches, applying instead the topological techniquesof persistent homology [12] to arrive at a trajectory classification, which relies on global topological information extracted
either from the trajectory data itself, or using a simplicial
model of the environment containing these trajectories. Ourcurrent work extends our efforts [29] that introduced the useof persistent homology as a trajectory clustering technique,and showing how filtrations of Delaunay-Čech complexes
can be utilized to cluster trajectories with fixed common
start and end points. This work was recently extended toa construction in [30] to accommodate trajectories withvarying start and end points. This construction however reliedon an enlarged simplicial complex construction, substantially
increasing the size of the required simplicial data structure.

Our present work generalizes [29] instead by means ofrelative persistent homology, allowing us to cluster trajectories with varying end-points based on the same simplicial
complexes used for the classification of trajectories with fixedstart and end points. The present work studies large datasetsof real world GPS traces as a particular source of trajectorydata, which, as well as trajectories extracted from video [32],forms a common real data-source studied for example inanomalous event detection problems [28].
A. Mathematical BackgroundThe key tool used in this work is the machinery of relative
persistent homology with coefficients in a field F, a variantof persistent homology [12]. This requires a few definitions
from Algebraic Topology, requiring mathematical techniques
which are difficult to fully explain and do justice in a shortpaper. The book [13] provides an excellent introduction – weshall only review the necessary key notation and ideas here.
We focus on persistence with binary coefficient field Z2 ={0, 1}, but more generally fields such as Zp for some prime
p or the field of rational numbers Q could be considered.
The binary field Z2 has the particular advantage of being
simple and efficiently implementable on a computer. In 2
dimensions, the field Zp allows us to distinguish, e.g. up top − 1-fold winding of trajectories around obstacles – Z2 inparticular only enables us to detect whether we move to theleft or right of voids/obstacles [16], [13].
1) Simplicial Complexes: A geometric k-simplex σ =[v0 , . . . , vk ] in Rd is a convex hull of k + 1 affinely independent ordered points v0 , . . . , vk ∈ Rd . We call k thedimension of a k-simplex. If τ ⊆ σ, τ, σ ∈ K, then τ is calleda face of σ. In the special case of Z2 homology, the ordering
can in fact be ignored. A geometric simplicial complex K
is a non-empty set of simplices such that if σ ∈ K and∅=6 τ ⊆ σ ∈ K, then τ ∈ K and the intersection of any two
simplices σ, τ ∈ K is a face of both σ and τ . We write |K|for set of points in Rd contained in the union of all simplices
in K. The set |K| is a topological space with the subspace
topology from Rd . A subset of simplices A ⊂ K that isitself a simplicial complex is called a subcomplex of K. Notethat 0-simplices just correspond to points, 1-simplices are
finite line segments and 2-simplices are triangles in Rd . The
conditions in the definition of a simplicial complex ensurethat the simplices are assembled in a natural manner and thenotion of a simplicial complex generalizes both the notionof a geometric graph and a triangulation.
2) Relative Homology:
For a field F, a p-chain c is aPkkλσformal sum c =i=1 i i of p-simplices {σi }i=1 ⊂ Kwith λi ∈ F and Cp (K) denotes the F-vector space of all pchains. In particular, for finite simplicial complexes, 1-chainsare finite linear combinations of edges and 2-chains are
finite linear combinations of triangles. When no confusion
arises, we write Cp for Cp (K) to simplify notation. For
every geometric p-simplex
Pp σ = i[v0 , . . . , vp ] let ∂p σ be thep − 1-chain ∂p σ =i=0 (−1) [v0 , . . . , vi−1 , vi+1 , . . . , vp ]consisting of a signed sum of faces of σ. For each p ∈{0, . . . , d}, ∂p extends to a linear map ∂ : Cp → Cp−1 ,called the boundary operator.
A p-chain c such that c = ∂p+1 ω for some ω ∈ Cp+1
is called a p-boundary, and a p-chain c such that ∂p c = 0
is called a p-cycle. The vector spaces of p-boundaries andp-cycles are denoted Bp , and Zp , respectively. For a 1-chainc corresponding to an oriented path from s to t, ∂c = t − s.
We have ∂c = 0, so that c is a 1-cycle, for any closedoriented 1-chain c. Similarly, the boundary ∂w of a 2-chain
w corresponding to an oriented collection of triangles (2simplices) corresponds to the oriented geometric boundaryof these triangles (however only modulo F) - the nameboundary operator is hence also geometrically natural.
The p-th homology group of K is defined by Hp =Hp (K) = Zp /Bp . For a cycle c ∈ Zp , we denote by [c] ∈ Hp
the resulting element in homology. Note that each p-cycle
c yields an element in Hp , but this representative is only
unique up to elements in Bp . We are interested in H1 inparticular, which consists of equivalence classes of closed
1-cycles up 1-cycles that are boundaries of 2-cycles. In thecase of F = Z2 , we can visualize 1-chains as a collectionof edges in K which have non-zero coefficients in the chain.
See Fig. 2 for an example. The importance of homology inmathematics arises from the fact that it captures global topological properties about the topological space defined by |K|.
In particular homology remains invariant under continuous
deformations of the space |K| (homotopy equivalences of|K|). In particular, bp = dim(Hp (K)) is called the pth Betti
number and counts the number of connected components(b0 ), tunnels (b1 ), and higher dimensional voids in |K|. The
left part of Fig. 2 illustrates an exemplary 1-cycle c lying in a
simplicial complex K and forming a basis of H1 (K), whichis in this case 1-dimensional and where we pick F = Z2
coefficients. There, |K| is in fact homotopy equivalent to acircle and dim(H1 (K)) = 1.
In our prior work [29], we used a basis for H1 (K) totopologically cluster trajectories α0 , . . . , αn , represented asedge paths in K. There, it was initially necessary to assume
that all trajectories had the same fixed start and end verticess, t ∈ K, respectively, in order to form Z2 -cycles ci =α0 + αi that could then be classified in homology. We were
only able to consider more general trajectories with general
initial and terminal regions S, T ⊆ K, by introducing a cone
construction [30], which increased the size and complexityof the approach.
Here we solve this problem with relative homology, which
provides a generalization of standard homology: For a subcomplex A ⊂ K, we define the quotient vector space ofrelative p-cycles Cp (K, A) = Cp (K)/Cp (A) The boundary
operator descends to a linear operator on relative chains.
We denote it by ∂ˆp [16]. We define the set of relative
p-boundaries by Bp (K, A) = im(∂ˆp+1 : Cp+1 (K, A) →Cp (K, A)). These correspond to p-chains c ∈ Cp (K) suchthat c = ∂p+1 w + a for some a ∈ Cp (A) and w ∈ Cp+1 (K).
Similarly, relative p-cycles are defined by Bp (K, A) =ker(∂ˆp : Cp (K, A) → Cp−1 (K, A)) and correspond to pchains c ∈ Cp (K) with ∂c ∈ Cp−1 (A). In particular, a 1-

Fig. 2. We display a simplicial complex K in light and dark blue and a
subcomplex A ⊂ K in dark blue. The closed red edge-path c1 on the left
is a 1-cycle, which cannot be represented as the boundary of a 2-cycle dueto the hole in K. In fact [c1 ] ∈ H1 (K) forms a basis for the 1-dimensional
first homology group H1 (K). On the right, we display a relative 1-cycle c2
in red. Both c1 , c2 in fact form relative 1-cycles and [c1 ], [c2 ] ∈ H1 (K, A)
yield a basis for the 2-dimensional first relative homology group H1 (K, A).
While |K| is homotopy equivalent (deformable) to a circle, the quotient
space K/A is homotopy equivalent to two circles glued at a common
point, called a wedge of two circles - for this imagine first gluing the two
shaded connected components of A together and then shrinking the resulting
cylinder with a cut out hole until we obtain two circles. Since a wedge of
2 circles has a 2-dimensional H1 (K/A), we can also reason geometricallyto understand why dim(H1 (K, A)) = dim(H1 (K/A)) = 2.

chain corresponding to a connected oriented edge-path in K
is a relative 1-cycle if the path is cyclic (zero boundary),or when its initial and terminal points lie in A. Finally, we
have Bp (K, A) ⊆ Zp (K, A) and the p-th relative homology
is defined by Hp (K, A) = Zp (K, A)/Bp (K, A), describing
equivalence classes of relative p-cycles modulo relative pboundaries. Relative homology is of importance since itallows us to consider properties of topological quotient
spaces. Consider a simplicial complex K and a subcomplex
A ⊂ K. We can consider the quotient space K/A, within
which a point [x] corresponds to an equivalence class of apoint x ∈ K modulo A and all of A ⊂ K is identified with asingle point. An important result [16] states in particular thatthe first reduced homology of K/A with field coefficientscan be computed using the relative homology H1 (K, A),
so that we can think of relative first homology in terms of
the first homology of the quotient space. The two parts ofFig. 2 illustrate a basis for H1 (K) and H1 (K, A), where A
is the subcomplex of K containing all simplices of K whose
vertices all lie in the shaded region (which is |A|). Here,dim(H1 (K)) = 1, while dim(H1 (K, A)) = 2.
3) Filtrations: Simplicial complexes have for many
decades enjoyed popularity in pure mathematics [16] inorder to model and approximate various spaces of interest.
Only recently however has there been interest in constructing
simplicial complexes from real world data which has leadto the development of persistent homology [12], [8] which
studies the homology of an increasing sequence (a filtration)of topological spaces. Persistence is related to Morse theory,
where one studies the topology of sublevel (or superlevel)sets of a function f : X → R defined on a topologicalspace X. Each sublevel set Xr = f −1 ((−∞, r]) yields a
topological space Xr , where Xr ⊆ Xr0 whenever r 6 r0 .
As r increases, homological features can be ‘born’ anddisappear or ‘die’ as the threshold r increases. Persistence

provides a computational mechanism for understanding thesechanges. To make this precise, we work with a filtration K
of finite simplicial complexes in Rd , by which we meana sequence K : K1 ⊂ K2 ⊂ . . . ⊂ Kn = K∞ offinite simplicial complexes. Typically, each filtration index
i is associated to a real valued filtration value r so that
Ki = f −1 (−∞, r]. For example, we can assign an arbitrary
real valued scalar to each vertex of Kn . Then the function
f (σ) = maxi∈{0,...,k} f (vi ) for an arbitrary k-simplex σ =[v0 , . . . , vk ] ∈ Kn yields a valid filtration function giving
rise to a filtration of simplicial complexes when its sublevelsets are considered. When Ki = f −1 ((−∞, r]) we call r thefiltration value associated to the filtration index i.

0.3

0.3

0.15

0.15

0

0
0

0.15

0.3

0

0.15

0.3

B. Simplicial complexes from data
To study homological features of a point-cloud dataset
X = {x1 , . . . , xn } ⊂ Rd , we need to construct a sequenceof topological spaces modeling X. We will in particularfamily of union of balls spaces Xr =S consider thedx∈X {y ∈ R : kx − yk 6 r}, for r > 0. For each r,
Xr is homotopy equivalent to the Delaunay-Čech complex
DCr (X) [3], which is a simplicial complex defined for any
finite set X ⊂ Rd where each subset of d + 1 point isaffinely independent. The assumption of affine independence
is generic in that a uniform random sample satisfies thiscondition with probability one and we can also enforce
the condition by an arbitrarily small perturbation of X.
Let D(X) denote the simplicial complex corresponding to
the Delaunay triangulation of X with simplices defined by
D(X) = {[v0 , . . . , vk ] : vi ∈ X, ∩ki=0 Vvi 6= ∅}, where Vx
denotes the Voronoi cell containing x. For each k-simplex
σ = [v0 , . . . , vk ] ∈ DC(X), define f (σ) = min{r :Tkdi=1 Br (vi ) 6= 0}, where Br (x) = {y ∈ R : kx − yk 6 r}.
The Delaunay-Čech complex DCr (X), for r > 0 is the subcomplex of D(X) defined by DCr (X) = f −1 ((−∞, r]).
Since DCr (X) is homotopy equivalent to Xr , we cancompute topological information about Xr from DCr (X) at
all scales r > 0. If a function f : X → R is defined on thedata X – for example a log likelihood, probability density or
cost function, we can furthermore define an induced filtration
by extending f to a function f (σ) = maxi∈{0,...,k} f (vi ) for
the k-simplex σ = [v0 , . . . , vk ] of D(X), yielding a filtration
F (f )r of simplicial complexes by the sub-levelsets of f . As
we increase the threshold parameter r a larger and larger subcomplex of D(X) is considered. Superlevel sets can also be
studied by replacing f with −f . Fig. 3 illustrates examplesof DCr (X) at various thresholds. Besides DC(X), there
furthermore exist alternate constructions yielding simplicial
complexes, such as the Vietoris-Rips and Witness complexes[8] which are scalable to high-dimensional spaces, but which
are not necessarily homotopy equivalent to Xr at filtrationlevel r.
C. Persistent relative homology
When we apply homology to a filtration of simplicial
complexes K1 ⊂ . . . ⊂ Kn , we obtain a sequence of linear
maps fpi,j : Hp (Ki ) → Hp (Kj ) for i 6 j induced by theFig. 3. Given the pointcloud X ⊂ [0, 1]2 in the bottom left, we consider
DC(X) whose first persistence diagram is displayed in the top left. A large
persistence interval (0.04, 0.30) is marked in red and corresponds to thehole in the point-cloud. The middle row displays DCr (X) at r = 0.04
when the large hole is first enclosed by simplices, at r = 0.15 and at
r = 0.30 when the hole finally is covered by simplices. The red path
is a 1-cycle c such that [c] corresponds to marked most persistent pointin the diagram. The bottom middle and right figure displays DCr (X, Y ),where Y ⊂ X is the set of points marked in red in the bottom left andfor r = 0.03 (middle) and r = 0.04. The top right figure shows thefirst persistence diagram of DCr (X, Y ) relative to the shaded dark blue
subcomplex A(X, Y ). The diagram differs only in minor detail from thepersistence diagram of DC(X) on the left - apart from the existence of anon-finite point (0.03, ∞) indicated by the dashed line. For r ∈ (0.04, 0.3),
dim(H1 (DCr (X, Y ), A(X, Y ))) = 2, as discussed in Fig. 2.

inclusions Ki ⊂ Kj . The p-th persistent homology groupfor i 6 j is given by Hpi,j = imfpi,j , so that non-trivial
elements in Hpi,j correspond to homology classes ‘born’at or before index i and which do not ‘die’ (i.e. do notmerge or become trivial) before index j. The differencej − i is called the index persistence of such a class. For
us, Ki = f −1 ((−∞, ri ]), and rj − ri is the persistence ofthe class. In fact, all the persistent homology groups canbe computed by a decomposition of the persistence module
into interval modules [11]. The p-th persistence diagram,
captures the information about the birth and death of p-th
homology classes as the filtration value increases. It consistsof multisets of points above the diagonal. Each point (ri , rj )in the diagram corresponds to a homology class born at index
i and surviving until index j. Points that lie far above thediagonal have a large persistence and are hence considered
important features distinct from smaller scale features due
to noise. An example is presented in Fig. 3. Classes born atindex i and which do not die at the final filtration index nare called essential and are associated to points of the form
(ri , ∞) in the plane, extended formally to (R ∪ {∞})2 . Notethat the dimension of H1 (Ki ) for any ri > 0 is equal to thenumber of points above and to the left of (ri , ri ). To computea basis for the persistent homology groups, we first assume

000
001
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026
027
028
029
030
031
032
033
034
035
036
037
038
039
040
041
042
043
044
045
046
047
048
049
050
051
052
053
054

without loss of generality that the filtration KS has beenj
refined to a simplex-wise filtration, where Kj = i=1 σi , sothat Kj+1 = Kj ∪{σj+1 } and we hence add a single simplex
in each step of the filtration. Given such a simplex-wise
filtration, several algorithms (see, e.g. [4]) are available tocompute a basis of the persistent homology groups. We shalluse the library [4] and the left-to right reduction algorithm
described in [10] for this purpose. Fig. 3 illustrates anexample of a filtration and an associated homology basis.
Persistence has recently emerged as a new approach in data
analysis, since large persistence intervals in the persistence
diagrams are provably stable under noise [12] and represent
global structure information currently only extractable with
topological methods. In this work, we will use a less common
extension of persistence to relative homology. For a subcomplex A ⊆ K1 of a filtration, we consider the sequenceof linear maps on relative homology fˆpi,j : Hp (Ki , A) →Hp (Kj , A) for i 6 j induced by the inclusions Ki ⊂ Kj . A
basis for these relative homology groups can be obtained by amodification of the standard left-to right reduction algorithm.
Since H1 (Ki , A) ' H1 (Ki /A), we can in particular thinkof the first relative persistent homology groups as measuring
the persistence of the sequence of topological quotient spaces
Ki /A for i ∈ {1, . . . , n}.
An approach we introduce in his work is to defineA(X, Y ) to be the sub-complex of DC∞ (X) = D(X)of simplices whose vertices are contained in Y (shaded inFig. 3). We then augment A(X, Y ) to a filtration DCr (X, Y )by inserting the remaining simplices of D(X), with filtration
order defined by the defining function f of DC(X). In the
top left of Fig. 3, we illustrate the resulting relative first
persistence diagram for H1 (DCr (X, Y ), A(X, Y )).
III. M ETHODOLOGY
As input, we consider a set Γ = {γ1 , . . . , γm } of piecewise
linear trajectories γi : [0, 1] → Rd which can have arisen for
example as GPS traces, or as joint-configuration trajectories
of a robotic system. We assume furthermore that these
trajectories all start and end in some specified (possibly
disconnected) region R ⊂ Rd .
A. Discretization and Setup
In order to apply our simplicial complex based approach,
we require a simplicial complex L such that each trajectory
γi can be discretized as a sequence of 1-simplices in L. Fora simplicial complex such that Γ ⊂ |L| (recall |L| denotes
the union of all simplices in L) this discretization can be
performed by first mapping each vertex of the piecewise
linear γ to its nearest 0-simplex in L and to then connect
each consecutive pair of such mapped trajectory points by ashortest path of 1-simplices in L. Given R, define A ⊆ L toconsist of those simplices of L whose vertices all lie in R,
and we assume that the mapped initial and terminal positionsof the trajectories also map to corresponding vertices in A
under our discretization procedure.
Given A, L, suppose now that there exists a filtration L =K1 ⊂ K2 ⊂ . . . ⊂ Kn of simplicial complexes. In that case,

we can consider the first persistent relative homology groups
arising from the inclusions
H1 (K1 , A) → . . . → H1 (Kn , A).
Each γ ∈ Γ forms a relative cycle in Z1 (Ki , A) for alli ∈ {1, . . . , n}. For each fixed i, we can furthermore consider[γ] ∈ H1 (Ki , A) = Z1 (Ki , A)/B(Ki , A). The key insightof persistence [12] however implies that we can compute abasis for all H1 (Ki , A), i = 1, . . . , n simultaneously, and thealgorithms of persistence [4] furthermore provide the mostefficient known implementation of determining these bases
even when i is fixed. To compute these bases, we first create
an arbitrary simplex-wise filtration of L = K1 and augment
this filtration by adding each of the simplices of Kn suchthat simplices of Ki are inserted before Kj and such that thefaces of any simplex σ are inserted before σ. This results
in a refined filtration M1 ⊂ . . . ⊂ Mp = Kn , where one
simplex is inserted at each step and A = Ms for some s.
By applying the standard left-to-right matrix reduction to theboundary matrix ∂ of this filtration, one arrives at a reduction
R = ∂V such that the required bases for H1 (Ki , A) consistsof a subset of lower parts of columns of R and V below row
s. The column indices are specified by the persistence indices(see also [29] for the non-relative version of this). For each
fixed filtration index i, we thus obtain a basis b1 , . . . , bw
of H1 (Ki , A) allowing us to compute F-coordinates of theimage of our trajectories in H1 (Ki , A).
Note that this procedure amounts to a finite dimensional
binary vectorial featurization of the trajectories.
However, unlike other known featurizations, the dimensionality of the features is independent of the length of the
trajectory and instead depends on the global topology of thequotient space Ki /A.
We propose to cluster trajectories with the same homology
coordinates into a joint cluster for each filtration index i.
As i is increased, these trajectory clusters then merge hierarchically until some final filtration value. The trajectories
however do not necessarily all merge into a single cluster
at the final filtration index (because H1 (K∞ , A) might not
be trivial). Note that, compared to [29], relative homology
has the additional benefit that each trajectory already forms
a relative cycle, while in the cited work, it was necessary
to apply concatenations of trajectories to form (non-relative)
cycles in H1 (Ki ).
It follows from the fact that H1 (Ki , A) ' H1 (Ki /A)that two trajectories γ, γ 0 such that [γ] 6= [γ 0 ] ∈ H1 (Ki , A)
are not continuously deformable to one another in Ki if weallow their endpoints to also vary continuously in |A| – thisa relaxed version of the standard notion of homotopy, where
the end-points need to be fixed. Our clusters hence consistof sets of trajectories such that no trajectory from one cluster
can be continuously deformed to any trajectory of the othercluster in this sense. The converse is however not necessarily
true: it can happen that two trajectories in the same cluster
cannot be continuously deformed to one another in |Ki /A|.
Besides this information, the clustering is hierarchical innature and we can extract from the persistence intervals thefiltration indices at which two trajectories remain in separate
clusters and when they merge.
1) Filtrations of interest: The freedom in choice of filtration allows us a wide modeling capability to express
constraints. Given just a point-cloud X and terminal subset Y ⊂ X, a natural choice is the filtration given by
DCr (X, Y ) introduced earlier and shown in Fig. 3. Using
the persistence intervals, we can automatically determine
filtration values such that our trajectories lie in |DCR (X, Y )|for some minimal R > 0 and such that DCR (X, Y ) ispath connected (0-persistence). We can furthermore identify
filtration regions where the remaining Betti numbers (numberof voids, tunnels, etc.) do not vary under noise (i.e. by
focusing on large persistence intervals) using the persistence
stability theory [12]. Two particular situations to distinguishare in particular the case where X is a large and dense
set of samples from some configuration space C, in whichcase DCR (X) can be expected to be homotopy equivalent
to C and the classification of trajectories depends on the
topology of C only. In the second case, X is sparse andmight only consist of trajectory points itself – in that casethe classification of trajectories captures global intrinsic
information about the shape of the union of these trajectories
themselves – recall here that DCr (X) is topologically thesame (homotopy equivalent) to the union of ball space Xrand hence recovers topological information of Xr as theseballs grow around the samples X. We will work with bothcases, as has been done in the non-relative case in [29].
Note that in the above, only the filtration function f defining DCR (X) was used. In general, however our filtration canarise as sublevel sets of an arbitrary function. In particular,
we can consider a probability density f , log-likelihoods, costfunctions, etc. Our classification approach enables us to encode topological constraints by considering the homologicalproperties of all sublevel sets simultaneously. The notion of
topology is hence rather generic in nature as a sublevel set ofsome probability density on a simple topological space, suchas a square, can have intricate features, while the square itself
is topologically trivial. Similarly DCR (X, Y ), for sparseX ⊂ [0, 1]2 , can capture information about X itself, rather
than the square’s topology.
IV. E XPERIMENTS
We work with a dataset of GPS ship trajectories [23]from a region around Shanghai, which was recorded during
6 months in 2013. We are interested in classifying the
motions of ships between the shaded northern and southern
region in Fig. 4. Our dataset consists of 11711 trajectories
with 5102072 GPS data points which corresponded to ships
traveling between these regions and such that between eachconsecutive GPS signal point along a trajectory there wasno more than a 30 minute delay. We now demonstrate how
our approach benefits from information about the geometryof the sea region. We hence assume the knowledge of asea/land classification of this area is available and create adetailed map by uniformly sampling 100000 points X ⊂ R2
lying in the water region. Denoting the subset of X lying in

Fig. 4. A dataset of 11711 ship GPS traces of vessels traveling between
the shaded regions around Shanghai is shown in the top figure. The bottom
figure displays DCr (X, Y ) and a subset of 500 trajectories discretized as
relative 1-cycles in the simplicial complex representing the water area at
r = 0.0002.

the shaded region by Y , we constructed a triangulation of
the water region (longitude [121.897, 122.281] and latitude
range [29.8754, 30.0128]) using DCr (X, Y ) for r = 0.0002.
Next we discretized the trajectories as approximate edgepaths in DCr (X, Y ) as described in Sec. III. See Fig. 4 foran illustration. Using the filtration DCr (X, Y ), we computed
relative persistent homology coordinates for all trajectories tocluster the dataset. The filtration DC(X, Y ) was constructed
in 1.15s and, once each trajectory was represented as a relative cycle, it took less than 0.002s per trajectory to determine
its relative persistent homology coordinates. We obtained 15
classes at r = 0.0002. All classes are displayed in the topof Fig. 5, and we can see that this clearly corresponds to aclassification of the ship movements relative to the various
islands. Note the differences with a metric classification:
our classification allows for outliers in this setting as it
depends only on topological properties of the water region.
The middle row of Fig. 5 displays two example classes. In
the bottom right, we see that a single linkage clustering byFréchet distance does not result in clustering representing
the ‘hard constraints’ enforced by the environment topology

Fig. 6. OpenStreetMap [26] dataset of GPS traces intersecting a boundingbox around a highway crossing. We display the discretization of thosetrajectories traveling between the indicated gray regions as well as thecomplex DCr (X, Y ), r = 0.000109. Note that some of the scatteredpoints in the complex are due to GPS datapoints of other trajectories in thebounding box which did not travel between the shaded regions.

Fig. 5. The top figure displays all 15 found clusters at r = 0.0002, while
the middle row illustrates two classes in isolation. Note how much variationthese classes exhibit - yet they are distinguishable by the ’hard environment
constraint’ posed by the small southern island that they pass in a distinct
manner. The bottom right displays a single linkage clustering by discreteFréchet distance at distance 0.055, yielding 7 completely environment
agnostic clusters. The bottom right figure shows the classification using our
method at a higher filtration value of r = 0.01, where the smaller islands
have been covered by simplices, resulting in only 3 trajectory classes at thatfiltration level.

as these methods do not incorporate such information. The
Fréchet clustering implementation took 57.7 minutes andrequired approximately 5.8GB memory, while our approach
took a total of less than 30s and less than 300M B memory.
While the 2D filtration and persistence diagram computation
can require worse than linear time in the number of samples|X| [12], recall that |X| was independent of the number oftrajectories in our experiment. Asymptotically, the topological clustering hence requires only constant memory and time
per discretized trajectory. As a result, the topological clusters
can be computed in linear time in the number of discretized
input trajectories for a fixed simplicial representation.
A. Application to highway traffic analysis
We now focus on a traffic analysis problem, where weassume that we are not given any road/land classification databut only the trajectories themselves. We consider the set of
1685 GPS traces of cars driving across a highway crossing
next to Frankfurt airport, obtained from OpenStreetMap

Fig. 7. Subset of 45 traces of pedestrian data from [2]. On the right, we
display the resulting classification based on the probabilistic collision model
shown on the left (in color) and at filtration value of r = 0.59.

[26]. The dataset contains 60987 GPS points X and wefurthermore select the subset of 177 trajectories that startand end in the shaded regions shown in Fig. 6. Using only
X, and the subset Y of points in the shaded regions, weconstruct DCr (X, Y ), displayed in Fig. 6, for r = 0.000109
and map the trajectories to edge-paths as before. Note that inthis construction |X|, and hence the asymptotic complexityof our clustering approach (construction of DCr (X, Y ) andpersistence computation), is dependent on the size of the
trajectory dataset, but one could introduce a fixed discretization of the data as before to bound |X| by a constant.
Classifying the trajectories using our approach, we obtain
24 trajectory classes as color-coded in Fig. 1. Some of theseclasses are separately displayed in the bottom part of thatfigure. As is visible in the figure, we are able to separate
driver behaviors with our approach - note for example theclover leaf trajectory class, where a driver crossed several
bridges before continuing along the same driving direction.
B. Incorporation of probability densities
In this experiment, we illustrate the use of a probabilistic
model in conjunction with our approach. Consider the indoor
scene in Fig. 7. We would like to model the behavior ofpedestrians in this space using a dataset of 197 trajectories
with 26029 datapoints X ⊂ R2 from [2]. We extract 45
trajectories traversing between the shaded regions, displayed

in black on the left. Besides these trajectories, we consider
a Gaussian Mixture Model M , modeling the positions ofpotential collision threats in this space (shown in color). We
construct a simplicial complex filtration P (X, Y ) startingwith the sub-complex A(Y ) ⊂ D(X) at filtration zero, as
before, but inserting simplices σ = [v0 , . . . , vk ] of D(X)in order of f (σ) = maxki=0 p(vi ), where p : R2 → R
denotes the probability density of the mixture model. The
first persistence diagram for this filtration has three large
persistence intervals: one corresponding to each mode of p,with the middle Gaussian resulting in the global maximumof p. All displayed trajectories have vertices lying in thecomplex Pr (X, Y ), for r = 0.5. At r = 0.59, we obtain 4
trajectory classes as indicated in the right figure, but as rincreases these classes start to merge, starting with the redand black classes which are separated by a rather weak modeof p. Our classification scheme in this case hence allows us
to understand the movement of pedestrians relative to thesuperlevel sets of a mixture model.
V. C ONCLUSIONS
We have proposed a novel topological trajectory clustering
approach and have focused on experiments with car, ship
and pedestrian data in 2D, showing that our approach can
extract interesting motion classes from large real-world trajectory datasets. In future work, we would like to investigate
higher dimensional trajectory data and are interested instudying the interplay between topological clustering andpolicy learning in particular. We are intending to makethe topological tools developed for this work available (seewww.csc.kth.se/˜fpokorny). A list of related papers
and software will also be maintained at the ICRA 2016
Workshop on Emerging Topological Techniques in Robotics
website: www.csc.kth.se/icra2016topology/.
VI. ACKNOWLEDGEMENTS
This research is supported in part by the Knut and Alice Wallenberg Foundation, the EU grant FLEXBOT (FP7-ERC-279933),
the NSF CISE Expeditions Award CCF-1139158, DOE Award
SN10040 DE-SC0012463, and DARPA XData Award FA8750-122-0331, and gifts from Amazon Web Services, Google, IBM, SAP,The Thomas and Stacey Siebel Foundation, Adatao, Adobe, Apple
Inc., Blue Goji, Bosch, Cisco, Cray, Cloudera, Ericsson, Facebook,
Fujitsu, Guavus, HP, Huawei, Intel, Microsoft, Pivotal, Samsung,
Schlumberger, Splunk, State Farm, Virdata and VMware.

R EFERENCES[1] B. D. Argall, S. Chernova, M. Veloso, and B. Browning, “A surveyof robot learning from demonstration,” Robotics and autonomous
systems, vol. 57, no. 5, pp. 469–483, 2009.
[2] B. Fisher, “Edinburgh informatics forum pedestrian database,” http://homepages.inf.ed.ac.uk/rbf/FORUMTRACKING/, 2015.
[3] U. Bauer and H. Edelsbrunner, “The morse theory of Čech andDelaunay filtrations,” in Proc. of the Thirtieth Annual Symp. on Comp.
Geometry, ser. SOCG’14. New York, NY, USA: ACM, 2014, pp.
484:484–484:490.
[4] U. Bauer, M. Kerber, and J. Reininghaus, “PHAT (Persistent Homology Algorithm Toolbox),” http://code.google.com/p/phat/.
[5] A. Billard, S. Calinon, R. Dillmann, and S. Schaal, “Robot programming by demonstration,” in Springer handbook of robotics. Springer,
2008, pp. 1371–1394.

[6] L. Brun, A. Saggese, and M. Vento, “A clustering algorithm oftrajectories for behaviour understanding based on string kernels,” inEighth Int. Conf. on Signal Image Technology and Internet Based
Systems (SITIS). IEEE, 2012, pp. 267–274.
[7] K. Buchin, M. Buchin, J. Gudmundsson, M. Löffler, and J. Luo,
“Detecting commuting patterns by clustering subtrajectories,” Int.
Journal of Computational Geometry & Applications, vol. 21, no. 03,
pp. 253–282, 2011.
[8] G. Carlsson, “Topology and data,” Bull. Amer. Math. Soc. (N.S.),
vol. 46, no. 2, pp. 255–308, 2009.
[9] CMU Graphics Lab, “CMU graphics lab motion capture database,”
http://mocap.cs.cmu.edu/, 2015.
[10] D. Cohen-Steiner, H. Edelsbrunner, and D. Morozov, “Vines andvineyards by updating persistence in linear time,” in Proc. of thetwenty-second annual symposium on Computational geometry. ACM,
2006, pp. 119–126.
[11] V. De Silva, D. Morozov, and M. Vejdemo-Johansson, “Dualities
in persistent (co) homology,” Inverse Problems, vol. 27, no. 12, p.
124003, 2011.
[12] H. Edelsbrunner and J. Harer, “Persistent homology-a survey,” Contemporary mathematics, vol. 453, pp. 257–282, 2008.
[13] H. Edelsbrunner and J. L. Harer, Computational topology: an introduction. AMS Bookstore, 2010.
[14] A. Fod, M. J. Matarić, and O. C. Jenkins, “Automated derivation ofprimitives for movement classification,” Autonomous robots, vol. 12,
no. 1, pp. 39–54, 2002.
[15] S. Gaffney and P. Smyth, “Trajectory clustering with mixtures ofregression models,” in Proc. of the fifth ACM SIGKDD Int. Conf. on
Knowledge Discovery and Data Mining. ACM, 1999, pp. 63–72.
[16] A. Hatcher, Algebraic topology. Cambridge: Cambridge University
Press, 2002.
[17] ICML, “Workshop: Topological methods for machine learning,” 2014.
[18] O. C. Jenkins and M. Mataric, “Automated derivation of behavior
vocabularies for autonomous humanoid motion,” in Proceedings of
the Second International Joint Conference on Autonomous Agents andMultiagent Systems. ACM, 2003, pp. 225–232.
[19] B. Kehoe, S. Patil, P. Abbeel, and K. Goldberg, “A survey ofresearch on cloud robotics and automation,” Automation Science andEngineering, IEEE Transactions on, vol. 12, no. 2, pp. 398–409, 2015.
[20] K. Kihwan, L. Dongryeol, and E. Irfan, “Gaussian process regression
flow for analysis of motion trajectories,” in IEEE ICCV.
IEEE
Computer Society, November 2011.
[21] R. A. Knepper, S. S. Srinivasa, and M. T. Mason, “Toward a deeper
understanding of motion alternatives via an equivalence relation on
local paths,” Int. Journal of Robotics Research, vol. 31, no. 2, pp.
167–186, 2012.
[22] J.-G. Lee, J. Han, and K.-Y. Whang, “Trajectory clustering: a partitionand-group framework,” in Proceedings of the 2007 ACM SIGMOD int.
conf. on Management of data. ACM, 2007, pp. 593–604.
[23] MarineTraffic, http://marinetraffic.com/, 2014.
[24] N. Nguyen, D. Phung, S. Venkatesh, and H. Bui, “Learning anddetecting activities from movement trajectories using the hierarchicalhidden markov model,” in CVPR 2005, vol. 2. IEEE, 2005, pp.
955–960.
[25] NIPS, “Workshop: Algebraic topology and machine learning,” 2012.
[26] OpenStreetMap, http://www.openstreetmap.org, 2015.
[27] P. Pastor, H. Hoffmann, T. Asfour, and S. Schaal, “Learning and generalization of motor skills by learning from demonstration,” in Robotics
and Automation, 2009. ICRA’09. IEEE International Conference on.
IEEE, 2009, pp. 763–768.
[28] C. Piciarelli, C. Micheloni, and G. L. Foresti, “Trajectory-based
anomalous event detection,” Circuits and Systems for Video Technology, IEEE Trans., vol. 18, no. 11, pp. 1544–1554, 2008.
[29] F. T. Pokorny, M. Hawasly, and S. Ramamoorthy, “Multiscale topological trajectory classification with persistent homology,” in Proceedingsof Robotics: Science and Systems, Berkeley, USA, July 2014.
[30] ——, “Topological trajectory classification with filtrations of simplicial complexes and persistent homology,” International Journal of
Robotics Research, 2015.
[31] C. Saunders, D. R. Hardoon, J. Shawe-Taylor, and G. Widmer, “Using
string kernels to identify famous performers from their playing style,”in Machine Learning: ECML 2004. Springer, 2004, pp. 384–395.
[32] H. Wang, A. Klaser, C. Schmid, and C.-L. Liu, “Action recognition
by dense trajectories,” in CVPR. IEEE, 2011, pp. 3169–3176.
[33] Y. Zheng and X. Zhou, Computing with spatial trajectories. Springer,
2011.