Efficient Greedy Discrete Subtrajectory Clustering

Comparing whole trajectories gives little information about shared structures. Consider two trajectories $T_1$ and $T_2$ that represent individuals who drive the same route to work, but live in different locations. For almost all distance metrics, the distance between $T_1$ and $T_2$ is large because their endpoints are (relatively) far apart. However, $T_1$ and $T_2$ share similar subtrajectories. Subtrajectory clustering allows us to model a wide range of behaviours: from detecting commuting patterns, to flocking behaviour, to congestion areas [12]. In the special case where trajectories correspond to traffic on a known road network, trajectories can be mapped to the road network. In many cases, the underlying road network might not be known (or it may be too dense to represent). In such cases, it is interesting to construct a sparse representation of the underlying road network based on the trajectories. The latter algorithmic problem is called map construction, and it has also been studied extensively [2, 6, 3, 14, 17, 25, 26]. We note that subtrajectory clustering algorithms may be used for map construction: if each cluster is assigned one centre, then the collection of centres can be used to construct an underlying road network [9, 10].

BBGLL [12] introduce the following problem. Rather than clustering $𝒯$, their goal is to compute a single cluster. Given $𝒯$, a $\mathrm{\Delta }$-cluster is any pair $(P,𝒫)$ where:

• $\mathrm{■}$

  $P$ is a subtrajectory of a trajectory in $𝒯$,

• $\mathrm{■}$

  $𝒫$ is a set of subtrajectories, where for all $𝒯[a,b],𝒯[c,d]\in 𝒫$, $[a,b]\cap [c,d]=\mathrm{\varnothing }$,

• $\mathrm{■}$

  and, for all $P^{\prime }\in 𝒫$, the Fréchet distance between $P$ and $P^{\prime }$ is at most $\mathrm{\Delta }$.

See Figure 1. They define two functions that each output a single $\mathrm{\Delta }$-cluster $(P,𝒫)$:

• $\mathrm{■}$

  SC($\max$, $\mathrm{\ell }$, $\mathrm{\Delta }$, $𝒯$) requires that $|P|=\mathrm{\ell }$ and maximises $|𝒫|$.

• $\mathrm{■}$

  SC($m$, $\max$, $\mathrm{\Delta }$, $𝒯$) requires that $|𝒫|=m$ and maximises $|P|$.

If $|P|=\mathrm{\ell }$ and $|𝒫|=m$, their algorithms use $O(n^2+nm\mathrm{\ell })$ time and $O(n\mathrm{\ell })$ space.

AFMNPT [1] propose a clustering framework. They define a clustering $C$ as any set of clusters. They propose a scoring function with constants $(c_1,c_2,c_3)$ that respectively weigh: the number of clusters, the radius of each cluster, and the number of uncovered vertices. Their implement a heuristic algorithm that uses this scoring function.

BBDFJSSSW [9] do map construction. They compute a set of $\mathrm{\Delta }$-clusters $C$ for various $\mathrm{\Delta }$. For each $(P,𝒫)\in C$, they use $P$ as a road on the map they are constructing. Thus, it may be argued that [9] does subtrajectory clustering. As a preprocessing step, they implement SC($\max$, $\mathrm{\ell }$, $\mathrm{\Delta }$, $𝒯$). Choosing various $\mathrm{\ell }$ and $\mathrm{\Delta }$, they compute a collection $S$ of “candidate” clusters. They greedily select a cluster $(P,𝒫)\in S$, add it to the clustering, and remove all subtrajectories in $𝒫$ from all remaining clusters in $S$. For a cluster $(Q,𝒬)\in S$, this may split a $Q^{\prime }\in 𝒬$ into two subtrajectories. BBGHSSSSSW [10] give improved implementations.

ABCD [4] do not require that for all $𝒯[a,b],𝒯[c,d]\in 𝒫$, $[a,b]\cap [c,d]=\mathrm{\varnothing }$. They do require that the final clustering $C$ has no uncovered vertices, and that for each cluster $(P,𝒫)$, $|P|\le \mathrm{\ell }$ for some given $\mathrm{\ell }$. For any $\mathrm{\Delta }$, denote by $k_{\mathrm{\Delta }}$ the minimum-sized clustering that meets their requirements. Their randomised polynomial-time algorithm computes a clustering of $O(\mathrm{\Delta })$-clusters of size $\tilde{O}(k_{\mathrm{\Delta }}\mathrm{\ell })$. BCD [8] give polynomial-time clustering algorithm of $O(\mathrm{\Delta })$-clusters of size $\tilde{O}(k_{\mathrm{\Delta }})$. Conradi and Driemel [16] implement an $\tilde{O}({\mathrm{\ell }}^2{n}^4+k_{\mathrm{\Delta }}\mathrm{\ell }{n}^4)$-time algorithm that computes a set of $O(\mathrm{\Delta })$-clusters of size $O(k_{\mathrm{\Delta }}\log n)$.

Our contribution is threefold. First, in Section 5, we describe an efficient implementation of SC($\max$, $\mathrm{\ell }$, $\mathrm{\Delta }$, $𝒯$) and SC($m$, $\max$, $\mathrm{\Delta }$, $𝒯$). We compare our implementation to all existing single-core implementations.

In Section 5 we propose that SC($\max$, $\mathrm{\ell }$, $\mathrm{\Delta }$, $𝒯$) produces high-quality clusters across all clustering metrics. We create a greedy clustering algorithm that fixes $(\mathrm{\ell },\mathrm{\Delta })$ and iteratively adds the solution $(P,𝒫)=$ SC($\max$, $\mathrm{\ell }$, $\mathrm{\Delta }$, $𝒯$) to the clustering (removing all vertices in $𝒫$ from $𝒯$). An analogous approach for SC($m$, $\max$, $\mathrm{\Delta }$, $𝒯$) gives two clustering algorithms.

In Section 6 we note that the functions SC($\max$, $\mathrm{\ell }$, $\mathrm{\Delta }$, $𝒯$) and SC($m$, $\max$, $\mathrm{\Delta }$, $𝒯$) search a highly restricted solution space, having fixed either the minimum length of the centre of a cluster or the minimum cardinality. We observe that these functions always output a cluster on the Pareto front of some bivariate function $\theta (m,\mathrm{\ell })$. Since points on this Pareto front give high-quality clusters, we argue that it may be beneficial to approximate all vertices of this Pareto front directly – sacrificing some cluster quality for cluster flexibility. We create a new algorithm PSC($\mathrm{\Delta }$, $𝒯$) that computes in $O(n^2{\log }^{4}n)$ time and $O(n)$ space a set of clusters that is a $2$-approximation of the Pareto font of $\theta (m,\mathrm{\ell })$. We show that using PSC($\mathrm{\Delta }$, $𝒯$) as a subroutine for greedy clustering gives even better results.

In Section 7 we do experimental analysis. We compare our clustering algorithms to AFMNPT [1] and BBGHSSSSSW [10]. We adapt the algorithm in [10] so that it terminates immediately after computing the subtrajectory clustering. We do not compare to Conradi and Driemel [16] as their setting is so different that it produces incomparable clusters. We compare the algorithms based on their overall run time and space usage. We evaluate the clustering quality along various metrics such as the number of clusters, the average cardinality of each cluster, and the Fréchet distance between subtrajectories in a cluster. Finally, we consider the scoring function by AFMNPT [1].

We first define discrete walks for pair of trajectories $P=(p_1,\mathrm{\dots },p_x)$ and $T=(t_1,\mathrm{\dots },t_y)$ in ${ℝ}^2$. We say that an ordered sequence $F$ of points in $[x]\times [y]$ is a discrete walk if consecutive pair $(i,j),(s,l)\in F$ have $s\in \{i,i-1\}$ and $l\in \{j,j-1\}$. The cost of a discrete walk $F$ is the maximum distance $d(P(i),Q(j))$ over $(i,j)\in F$. The discrete Fréchet distance is the minimum-cost discrete walk from $(x,y)$ to $(1,1)$:

Given $𝒯$ and a value $\mathrm{\Delta }$, the Free-space matrix $M_{\mathrm{\Delta }}(𝒯,𝒯)$ is an $n\times n$ $01$-matrix where for all $(i,j)\in [n]\times [n]$ the matrix has a zero at position $i$, $j$ if $d(𝒯(i),𝒯(j))\le \mathrm{\Delta }$. We use matrix notation where $(i,j)$ denotes row $i$ and column $j$. In the plane, $(0,0)$ denotes the top left corner. We also use $M_{\mathrm{\Delta }}(𝒯,𝒯)$ to denote the digraph where the vertices are all $(i,j)\in [n]\times [n]$ and there exists an edge from $(i,j)$ to $(s,k)$ whenever (see Figure 2):

• $\mathrm{■}$

  $s\in \{i,i-1\}$ and $k\in \{j,j-1\}$, and

• $\mathrm{■}$

  $M_{\mathrm{\Delta }}(𝒯,𝒯)[i,j]=M_{\mathrm{\Delta }}(𝒯,𝒯)[s,k]=0$.

Alt and Godeau [5] show that for any two subtrajectories $𝒯[a,b]$ and $𝒯[c,d]$:

${𝒟}_F(T[a,b],T[c,d])\le \mathrm{\Delta }$ if and only if there is a directed path from $(b,d)$ to $(a,c)$.

Note that we follow [1, 9, 10, 12] and consider only disjoint clusterings, meaning $\forall (P,{𝒫}^{\prime })\in C$, $\forall 𝒯[a,b],𝒯[c,d]\in 𝒫$, the intervals $[a,b]$ and $[c,d]$ are disjoint. For brevity, we say that a $\mathrm{\Delta }$-cluster is any cluster $(P,𝒫)$ where $\forall P^{\prime }\in 𝒫:{𝒟}_F(P,P^{\prime })\le \mathrm{\Delta }$.

AFMNPT [1] propose a cost function to determine the quality of a clustering $C$. We distinguish between $k$-centre and $k$-means clustering:

We also consider unweighted metrics: the maximal and average Fréchet distance between subtrajectories in a cluster, and the maximal and average cardinality of the clusters.

All competitor algorithms do not terminate on the larger real-world data sets, even when the time limit is 24 hours. Thus, we primarily use these sets to compare algorithmic efficiency. Our qualitative analysis uses the subsampled data sets and synthetic data. We generate data by creating a random domain $X$. Then, for some parameter $c$, we randomly select $c\%$ of the domain and invoke an approximate TSP solver. We add the result as a trajectory. This results in a collection of trajectories with many similar subtrajectories. Synthetic-A has a random domain of 200 vertices. Synthetic-B has a random domain of 100 vertices. For both of these, we consider $c\in \{50,90,95\}$.

We sketch how we optimise the BBGLL algorithm in the full version Considers two integers $a,b$ with $\mathrm{\ell }=b-a$ (Figure 3). It alternates between StepFirst which increments $a$, StepSecond which increments $b$, and Query that given $(\mathrm{\Delta },\mathrm{\ell }=b-a,m)$ aims to find a $\mathrm{\Delta }$-cluster $(P,𝒫)$ with $P=𝒯[a,b]$ and $|𝒫|=m$. We apply two ideas:

1. 1.

   Windowing. By [12], it suffices to keep only $M_{\mathrm{\Delta }}(𝒯[a,b],𝒯)$ in memory. Furthermore, one may store only the $z^{\prime }\in O(n\mathrm{\ell })$ zeroes in this sub-matrix. This is implemented in [9, 10]. By cleverly abandoning memory, our StepFirst deletes a row in constant time. The StepSecond function is then the bottleneck, as it adds a new row.

2. 2.

   Row generation. In [12], the StepSecond and Query functions iterate over all $n$ cells in the bottommost row of $M_{\mathrm{\Delta }}(𝒯,𝒯)$. By applying a range searching data structure, we instead only iterate over all zeroes in this row.

We note that map-construct-rtree [10] uses r-trees for row generation (we refer to the journal version [11]). However, we believe that the size of their r-tree implementation dominates the time and space used by their algorithm in our experiments. In the full version, we describe how to use a light-weight range tree to implement this principle.

We compare the implementations across our data sets for a variety of choices of $\mathrm{\Delta }$ and $m$. We use the implementations of SC($m$, $\max$, $\mathrm{\Delta }$, $𝒯$) so that our experiments have to fix fewer parameters. Figures 4 and 5 shows the time and space usage across a subset of our experiments. The data clearly shows that our implementation is significantly more efficient than the previous single-core implementations. It is often more than a factor 1000 more efficient in terms of runtime and memory usage. We do not provide a more extensive analysis of these results since this is not our main contribution.

Let $B(𝒯)$ be a balanced binary tree over $𝒯$. Let $S(\mathrm{\Delta },𝒯)$ start out as empty. Our algorithm PSC$(\mathrm{\Delta },𝒯)$ (defined by Algorithm 3) does the following:

• $\mathrm{■}$

  For each inner node in $B$, denote by $𝒯[a,b]$ the corresponding subtrajectory of $𝒯$.

• $\mathrm{■}$

  For every prefix or suffix $P$ of $𝒯[a,b]$, let $(P,𝒫)$ be a $\mathrm{\Delta }$-cluster maximising $|𝒫|$.

• $\mathrm{■}$

  We add $(P,𝒫)$ to $S(\mathrm{\Delta },𝒯)$.

Given a trajectory $𝒯[a,b]$ and our current set $S(\mathrm{\Delta },𝒯)$, we simultaneously compute for all prefixes $P$ of $𝒯[a,b]$ a maximum-cardinality cluster $(P,𝒫)$ and add this cluster to $S(\mathrm{\Delta },𝒯)$. We define an analogue of Definition 1:

Given $𝒯[a,b]$ we initialize our algorithm by storing integers $a$ and $c\leftarrow a$. We increment $c$ and maintain for all $y\in [n]$ the value ${\text{col}}_a(c,y)$. This uses $O(n)$ space, as does storing row $c$. Given this data structure, we alternate between a StepSecond and Query subroutine:

StepSecond: We increment $c$ by $1$. For each $y\in [n]$, we check if $d(𝒯(c),𝒯(y))>\mathrm{\Delta }$. If it is then $(c,y)$ is an isolated vertex in the graph $M_{\mathrm{\Delta }}(𝒯,𝒯)$. Thus, there exists no directed path from $(c,y)$ to any vertex in row $a$ and we set ${\text{col}}_a(c,y)=-\mathrm{\infty }$. Otherwise, ${\text{col}}_a(c,y)=\max \{{\text{col}}_a(\alpha ,\beta )\mid \alpha \in \{c,c-1\},\beta \in \{y,y-1\}\}$ which we compute in $O(1)$ time.

Query: For a prefix $P=𝒯[a,c]$ of $𝒯[a,b]$ we compute a maximum-cardinality $\mathrm{\Delta }$-cluster $(P,𝒫)$ using the values ${\text{col}}_a(c,y)$ (see Algorithm 4).

Note that Algorithm 4 is essentially the same as the algorithm in [12] that we described in Section 5. Thus, it computes a maximum-cardinality cluster. The key difference is that we, through ${\text{col}}_a(b,j_2)$, avoid walking through the free-space matrix in $O(\mathrm{\ell })$ time to compute $j_1$. We denote by Query’($T[c,b]$, $\mathrm{\Delta }$, $𝒯$, $S(\mathrm{\Delta },𝒯)$) the function which computes for a suffix $P=𝒯[c,b]$ of $𝒯[a,b]$ a maximum-cardinality cluster $(P,𝒫)$ in an analogous manner.

Consider any $\mathrm{\Delta }$-cluster $(P,𝒫)$. Let $P=𝒯[a,b]$. Let $𝒯[i,k]$ be the minimal subtrajectory in our tree that has $𝒯[a,b]$ as a subtrajectory. Let the children of $𝒯[i,k]$ in our tree be $𝒯[i,j]$ and $𝒯[j+1,k]$. Then, by minimality of $𝒯[i,k]$. Hence, $𝒯[a,j]$ is a suffix of $𝒯[i,j]$, and $𝒯[j+1,b]$ is a prefix of $𝒯[j+1,k]$. There must exist a trajectory $P^{\prime }\in \{𝒯[a,j],𝒯[j+1,b])\}$ which contains at least half as many vertices as $T[a,b]$. Since $P^{\prime }$ is either a prefix or suffix of a subtrajectory in our tree, we must have added some maximum-cardinality cluster $(P^{\prime },{𝒫}^{\prime })$ to $S(\mathrm{\Delta },𝒯)$. Moreover, since $P^{\prime }$ is a subtrajectory of $P$, $|{𝒫}^{\prime }|\ge |𝒫|$ (since any discrete walk from row $b$ to row $a$ is also a discrete walk in the Free-space matrix from row $d$ to row $c$ for $[c,d]\subseteq [a,b]$). $◀$

Consider the balanced binary tree $B(𝒯)$. Each level of depth of $B(𝒯)$ corresponds to a set of subtrajectories that together cover $𝒯$. For every such subtrajectory $𝒯[a,b]$, both StepSecond and Query get invoked for each row $c\in [a,b]$. Since the tree has depth $\log n$, it follows that we invoke these functions $O(\log n)$ times per row $c$ of $M_{\mathrm{\Delta }}(𝒯,𝒯)$ (i.e. $O(n\log n)$ times in total). Naively, both the StepSecond and Query functions iterate over all columns in their given row $c$. However, just as in Section 5, we note that our functions skip over all $(j_2,c)$ for which $d(𝒯(j_2),𝒯(c))>\mathrm{\Delta }$. It follows that we may use our row generation trick to spend $O(Z_c{\log }^{3}n)$ time instead ($Z_c$ denotes the number of zeroes in row $c$). Since each row $c$ gets processed $O(\log n)$ times, this upper bounds our running time by $O(z{\log }^{4}n)$.

For space, we note that invoking StepSecond and Query for a row $c$ need only rows $c$ and $c-1$ in memory. So the total memory usage is $O(n)$. Explicitly storing $S(\mathrm{\Delta },𝒯)$ takes $O(n^2\log n)$ space. However, since we only iterate over $S(\mathrm{\Delta },𝒯)$, the total space usage remains $O(n)$. The range tree therefore dominates the space with its $O(n\log n)$ space usage. $◀$

For $k$-means clustering, we obtain three algorithms by running Algorithm 2 instead. Note that our clustering algorithms invoked Algorithm 1 across several runs for exponentially scaling $\mathrm{\Delta }$. Our clustering algorithms do not invoke Algorithm 2 in this manner. Instead, Algorithm 2 considers exponentially scaling $\mathrm{\Delta }$’s internally. Our experiments show that this heavily impacts the clustering’s performance on large input.

We again obtain ten algorithm configurations. In addition, we run the map construction algorithm from [10] and the $k$-means subtrajectory clustering algorithms from [1]:

• $\mathrm{■}$

  Envelope runs the algorithm from [1] with the scoring parameters $(c_1,c_2,c_3)$.

We design different experiments to measure performance and the scoring function. For performance, we run fourteen experiments across the five real-world data sets. Seven experiments for $k$-centre clustering, and seven for $k$-means clustering. Each experiment ran three times, using a different vector $(c_1,c_2,c_3)$. Each experiment runs $11$ programs for $k$-centre clustering, or $12$ for $k$-means clustering. If a program fails to terminate after 24 hours, we do not extract any output or memory usage.

We restrict scoring measurements to the three smaller or subsampled real-world data sets. This is necessary to make Map-construct and Envelope terminate. We add six synthetic data sets. We only compare under $k$-means clustering since only our algorithms can optimise for $k$-centre clustering. The experiments were conducted on a machine with a 4.2GHz AMD Ryzen 9 7950X3D and 128GB memory. Over 90 runs hit the timeout. Running the total set of experiments sequentially takes over 100 days on this machine.

The full results per experiment are organised in tables in the full version. Here, we restrict our attention to only three plots.

Recall that $c_1$ weighs the number of clusters in the clustering and that $c_3$ weighs the fraction of uncovered points. We follow the precedent set in [1] and choose $c_1=1$, and $c_3$ equal to the number of trajectories in $𝒯$. Finally, $c_2$ weighs the Fréchet distance across clusters. For each data set, we first find a choice for $c_2$ such the algorithms lead to a non-trivial clustering. We then also run the algorithms using the scoring vectors $(c_1,10\cdot c_2,c_3)$ and $(c_1,0.1\cdot c_2,c_3)$, for a total of three scoring vectors per data set.

To compare solution quality, we measure the scoring function and other qualities such as: the total number of clusters in the clustering, the maximum Fréchet distance from a subtrajectory to its centre, and the average Fréchet distance across clusters. Our competitor algorithms frequently time out on the real-world data sets. Hence, the qualitative comparison also makes use of our synthetic data set. We use three different configurations of our synthetic set to create input of different size. We choose to compare only under $k$-means clustering for a more fair comparison to other algorithms. Indeed, the algorithm by [1] was designed for the $k$-means metric ([10] is scoring function oblivious). Consequently, comparing under $k$-centre would only give us an unfair advantage.